CONVERGENCE OF MIXED METHODS IN CONTINUUM MECHANICS AND FINITE ELEMENT ANALYSIS FAROOQUE AQUIL MIRZA B.Sc, University of Karachi, 1966 B.Eng., McGill University, 1971 M.Eng., University of British Columbia, 1973 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of C i v i l Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA January, 1977 (^ cT) Farooque Aquil Mirza, 1977. 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 i t 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/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 C i v i l Engineering The University of British Columbia Vancouver, B.C., V6T 1W5 Canada January, 1977 CONVERGENCE OF MIXED METHODS IN CONTINUUM MECHANICS AND FINITE ELEMENT ANALYSIS ABSTRACT The energy convergence of mixed methods of approximate analysis for problems involving linear self-adjoint operators is investigated. A new energy product and the associated energy norm are defined for such indefinite systems and then used in establishing the strain energy con-vergence and estimation of error for problems in continuum mechanics. In the process, the completeness requirements are laid out for approximate solutions. Also established is the mean convergence of the basic variable e.g. displacements and stresses. After accomplishing a new mathematical framework for the mixed methods i n continuum, the theory is then extended to the f i n i t e element method. The completeness requirements, convergence c r i t e r i a and the effect of continuity requirements on convergence are established. The f l e x i b i l i t y offered by the mixed methods in incorporating the boundary con ditions is also demonstrated. For stress singular problems, the strain energy convergence is established and an energy release method for deter-mining the crack intensity factor K^. is presented. A detailed eigenvalue-eigenvector analysis of the mixed f i n i t e element matrix is carried out for various combinations of interpolations for the plane stress linear elasticity and the linear part of the Navier-Stokes equations. Also discussed is i t s relation to the completeness requirements. Finally, numerical results are obtained from applying the mixed f i n i t e element method to several examples. These include beam bending, i i i a plane stress square plate with parabolically varying end loads, a plane stress cantilever and plane strain stress concentration around a circular hole. A plane stress example of a square plate with symmetric edge cracks is also solved to study the strain energy convergence. Lastly, two rec-tangular plates, one with symmetric edge cracks and the other with a central crack are considered to determine the crack intensity factor K^.. In most of the examples, the strain energy convergence rates are predicted and compared with the numerical results, and excellent agreement is observed. iv TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS iv LIST OF TABLES v i i LIST OF FIGURES x NOTATION x i i i ACKNOWLEDGEMENTS xiv CHAPTER 1. INTRODUCTION 1 1.1 Background 1 1.2 Purpose and Scope 3 1.3 Limitations 4 2. MATHEMATICAL PRELIMINARIES 5 2.1 Basic Concepts 5 2.2 The Limit Process in Hilbert Spaces 7 2.3 Orthogonality and Orthonormal Basis 8 2.4 Subspaces and Projections 9 2.5 Application to Linear Elasticity 12 3. MIXED METHODS 15 3.1 Mixed Variational Principle 16 3.2 Projection Operators 20 3.3 Mixed Galerkin Method 22 3.4 Concepts and Definitions 24 3.5 Convergence of the Mixed Method 31 3.6 Completeness 39 V CHAPTER Page 3.7 Estimation of Error in the Energy Product 39 3.8 Extension 41 4. FORMULATION AND CONVERGENCE OF THE MIXED FINITE ELEMENT METHOD 47 4.1 Generation of a Finite Element Approximation . . . 47 4.2 Generation of an Element Matrix 51 4.3 Assemblage of Element Matrices 53 4.4 Convergence Criteria 55 4.5 Completeness Criterion 58 4.6 Stress Singularities 67 4.6.A Strain Energy Convergence for Problems with Stress Singularity 68 4.6.B Determination of Stress Intensity Factor 73 5. APPLICATION OF BOUNDARY CONDITIONS 78 5.1 Homogeneous Boundary Conditions 80 5.2 Homogeneous Mixed Boundary Conditions 81 5.3 Nonhomogeneous Boundary Conditions 83 5.4 Boundary Residual Concept . . . . . 91 6. EIGENVALUE ANALYSIS OF THE ELEMENT MATRIX 93 6.1 Linear Ela s t i c i t y Problem 93 6.2 Linear Part of the Navier-Stokes Equations . . . . 103 7. APPLICATIONS OF THE MIXED FINITE ELEMENT METHOD . . . . 118 7.1 Beam Problem 118 7.1. A Two Second Order Equations 119 7. I.B Four First Order Equations 125 7.2 Plane Linear Ela s t i c i t y 128 7.2. A Plane Stress: Square Plate with Parabolically Varying End Loads 129 v i CHAPTER Page 7.2.B Cantilever (Plane Stress) 130 7.2. C Stress Concentration around a Circular Hole (Plane Strain) 133 7.3 Stress Singularities 134 7.3. A Strain Energy Convergence 134 7.3.B Evaluation of the Stress Intensity Factor K^. 136 8. CONCLUSIONS 138 BIBLIOGRAPHY 142 APPENDIX A: CALCULATION OF THE FINITE ELEMENT MATRIX; LINEAR DISPLACEMENTS AND STRESSES OVER A TRIANGULAR ELEMENT 206 APPENDIX B: EIGENVALUE DISTRIBUTION FOR A REAL SYMMETRIC MATRIX 212 APPENDIX C: EQUIVALENCE OF ENERGY PRODUCTS FOR FOUR FIRST ORDER BEAM EQUATIONS AND PLANE LINEAR ELASTICITY EQUATIONS 219 APPENDIX D: ELASTICITY SOLUTION FOR A CANTILEVER 223 v i i LIST OF TABLES Table(s) Page I Eigenvalues and eigenvectors of element matrix for t r i -angular elements using different combinations of interpolations for U,V,T'S; linear elasticity plane stress 146 II Eigenvalues and eigenvectors of element matrix for rectangular elements using different combinations of interpolations for u,v, and T'S; linear el a s t i c i t y plane stress 147 III Eigenvalues and eigenvectors of element matrix for t r i -angular elements using different combinations of interpolations for u,v,p and x's (x ,x ,x ); two-xx yy xy dimensional, incompressible creeping flow (linear part of the Navier-Stokes equations) 148 IV Eigenvalues and eigenvectors of element matrix for rec-tangular elements using different combinations of interpolations for u,v,p and x's (x ,x ,x ); two-r r xx yy xy dimensional, incompressible creeping flow (linear part of the Navier-Stokes equations) 149 V Numerical results for two second order beam equations; v-linear, M-linear . . . 150 VI Numerical results for two second order beam equations; v-quadratic, M-quadratic 152 VII Numerical results for two second order beam equations; v-cubic, M-cubic 154 v i i i Table(s) Page VIII Numerical results for four f i r s t order beam equations; forced boundary conditions on v and M; v,6,M and V a l l linear 156 IX Strain energy estimations for four f i r s t order beam equations with forced boundary conditions on v and 9 ; M,V,6, and v a l l linear 158 X Numerical results for four f i r s t order beam equations; forced boundary conditions on v and 0; shear strain energy term V included; M,V,9,v a l l linear 158 XI . Numerical results for parabolically loaded plane stress problem 160 XII Numerical results for the cantilever (plane stress) with boundary conditions B.C.I (Figure 19). Mixed f i n i t e element; displacements and stresses linear 161 XIII Numerical results for the cantilever (plane stress) with boundary conditions B.C.2 (Figure 19). Mixed f i n i t e element; displacements and stresses linear 161 XIV Comparison amongst C.S.T., L.S.T., Q.S.T., and mixed fi n i t e element. Cantilever with boundary conditions B.C.2 (Figure 19) 162 XV Numerical results for the plane stress problem of square plate with edge cracks, Figure 29 163 XVI Stress intensity factors from the fi n i t e element analysis of the rectangular plate with symmetric edge cracks, Figure 33(a) 164 XVII Stress intensity factors from the f i n i t e element analysis of the rectangular plate with a central crack, Figure 33(b) 164 ix Table(s) Page XVIII Comparison of stress intensity factors obtained from energy release rate using different elements and procedures 165 X LIST OF FIGURES Figure(s) Page 1 Mode I crack (opening mode) 166 2 Typical contour for evaluation of J-integral . . . . . . 166 3 Accommodation of crack extension Aa by advancing nodes on the path T 0 ' 167 4 Node numbers and degrees of freedom for a triangular element 168 5 Node numbers and degrees of freedom for a rectangular element 168 6 Mechanism for the fourth zero eigenvalue 169 7 Forces acting on an infinitesimal beam element 169 8 Forced boundary conditions on v and M for the beam problem 170 9 Degrees of freedom for the beam element—two second order equations 171 10 Two second order beam equations; displacement linear and moment linear 172 11 Two second order beam equations; displacement quadratic, moment quadratic 176 12 Two second order beam equations; displacement cubic and moment cubic 177 13 Four f i r s t order beam equations, forced boundary conditions on v and M; v,M,9 and V a l l linear 178 14 Four f i r s t order beam equations, forced boundary conditions on v and 0; v,0,M and V a l l linear . . 180 15 Four f i r s t order beam equations, forced boundary conditions on v and 0; shear strain energy included 181 xi Figure(s) Page 16 Parabolically loaded plane stress problem (N=4) 183 17 Parabolically loaded plane stress problem using mixed element; displacements and stresses linear • . 184 18 Parabolically loaded plane stress problem using displace-ment element (cubic displacements), Ref. [5] 186 19 Linear el a s t i c i t y cantilever problem and forced boundary conditions used 187 20 Convergence plots for the cantilever with boundary condi-tions B.C.I, using mixed f i n i t e element; stresses and displacements a l l linear 188 21 (a): Grids used for the constant stress triangular elements . 189 (b) : Grids used for the linear stress triangles . . . . 190 (c) : Grids used for the quadratic stress triangles . . . 191 22 Typical mesh used for mixed f i n i t e element (N=4) . . . . 192 23 Plots for the cantilever with boundary conditions B.C.2 using mixed f i n i t e element; stresses and displacements a l l linear 192 24 Convergence of tip deflection for C.S.T., L.S.T., Q.S.T. and mixed f i n i t e element with stresses and displace-ments linear for boundary conditions B.C.2 194 25 Model for i n f i n i t e plate by a square plate with a circu-lar hole at centre, isotropic and orthotropic 195 26 Finite element grid for the square plate with a circular hole in the middle. Isotropic and orthotropic cases 196 27 Comparison of theoretical and f i n i t e element (mixed) results for i n f i n i t e plate with a circular hole in the middle 197 x i i Figure(s) Page 28 Comparison of theoretical and constant stress f i n i t e element, Zienkiewicz [42], results with the same material properties as used in the mixed method; TO—1.0 197 29 Plane stress problem considered for investigation of energy convergence in the case of stress singularities 198 30 Plots of strain energy versus the mesh size for the plane stress problem with symmetric edge cracks—Figure 29 199 31 Strain energy convergence for the plane stress problem; square plate with symmetric edge cracks—Figure 29 • . 200 32 Normal stress distribution along the middle of a square plate with symmetric edge cracks—Figure 29(b), (N=12) 201 33 Rectangular plates with cracks used for determining the crack intensity factor K. 202 34 Finite element mesh for determining the crack intensity factor Kj., used for both symmetric edge and central cracks 203 35 Normal stress distribution along the edge OA of the rec-tangular plates with cracks—Figure 33 204 36 Boundary tractions and conditions for cantilever linear elasticity solution (Appendix D) 205 x i i i NOTATION The specific use and meaning of symbols is defined in the text where they are introduced. The summation convention holds for subscripted variables with repeated lower case indices; i t does not apply to repeated upper case indices. The range of summation is indicated where the variables are f i r s t introduced. The lower case letter Z which is also frequently used for element length or diameter appears as 1 in the text and the tables and as £ in the figures. The Greek symbol e implies "belongs to" unless otherwise specified. xiv ACKNOWLEDGEMENTS The author wishes to express his gratitude to his advisor, Dr. M. D. Olson, for his invaluable advice and guidance during the research and preparation of this thesis. He also expresses thanks to Dr. D. L. Anderson and Dr. N. D. Nathan for their valuable suggestions during the research work. The financial support of the Defence Research Board and the National Research Council of Canada, University of British Columbia Fellowships (1973-75) and Killam Predoctoral Fellowship (1975-76) is gratefully acknowledged. Finally, many thanks to Ruth Mirza for tremendous inspiration, assistance and diligent typing of this thesis. 1 CHAPTER 1 INTRODUCTION 1.1 Background The f i n i t e element method has proved to be an extremely powerful tool for the analysis of engineering problems for which no closed form solu-tions exist. A good introduction to the subject, which is s t i l l undergoing continuing development, can be found in the books by Zienkiewicz [41], Gallagher [9], Cook [4], Strang and Fix [33], Huebner [14], and Oden [19], etc. For more than a decade now, the use of either displacement or equilibrium f i n i t e element methods have dominated in the analysis of problems in continuum and structural mechanics. These methods precipitate from well-established extremum principles; the displacement method from the principle of minimum potential energy or Rayleigh-Ritz method and the equilibrium method from the principle of minimum complementary energy. Both energy principles mentioned here, involve positive definite operators with well developed mathematical properties, such as lower and upper bounds, complete-ness requirements, convergence in the energy sense, etc., and provided tremendous assistance in establishing the convergence of the displacement and equilibrium f i n i t e element methods. Hutton and Anderson [16] used the Galerkin method in the f i n i t e element analysis and established the convergence properties for the problems which may not have any variational principle, i.e. for non-self-adjoint operators. This encouraged many researchers to use the Galerkin method in the f i n i t e element analysis of incompressible viscous flow problems with velocities and pressure as basic dependent variables, Taylor and Hood [34], Olson and Tuann [26], etc. But the use of pressure as a dependent variable amounts to a mixed formulation which has caused some di f f i c u l t i e s regarding the degree of the polynomial assumed for pressure within the element in relation to the polynomials used for approximating the 2 velocity f i e l d . This is also observed in the derivation of the variational theorem for incompressible and nearly incompressible materials by Herrmann [13] where the pressure, taken as the mean of the normal stresses, is consider-ed as a dependent variable along with the displacements. The earlier appli-cations of the mixed method in the f i n i t e element analysis date back to the mid-1960's when the use of mixed f i n i t e element models for plate bending were proposed, independently, by Herrmann [12] and Hellan [10]. These involved the simultaneous approximation of two dependent variables, the bending moment and the transverse deflection of thin elastic plates and were based on station-ary rather than extremum variational principles. Dunham and Pister [6] used the Hellinger-Reissner variational principle to develop mixed f i n i t e element models for plane elasticity and plate bending problems while Wilnderlich [40] exploited the idea of mixed models in a f i n i t e element analysis of non-linear shell behaviour. Somewhat similar to mixed f i n i t e elements was the develop-ment of the hybrid elements by Pian and Tong [28] and Tong [36]. However the assumed approximations for stresses and displacements were not considered over the entire domain as in the mixed f i n i t e elements mentioned above. In a l l these studies, higher accuracies as well as rapid convergence were obtained for certain quantities, e.g., stresses, than from the corresponding displace-ment elements. However, the mixed methods involve indefinite systems and despite their wide spread use their mathematical properties are not as well understood as those of the displacement and equilibrium methods. This situation has caused considerable d i f f i c u l t i e s in establishing completeness requirements and convergence proofs for the mixed methods. Consequently the theoretical basis for these methods is far behind that for the displacement and equilibrium methods. See for example Tong and Pian [37] and Oliveira [25]. More recently, Oden [20] has discussed some generalizations of the 3 theory of mixed methods and Reddy and Oden [30] established the convergence of dependent variables, i.e., stresses and displacements, by decomposing the linear, positive definite operator into two self-adjoint linear operators (A=T*T) and then applied the theory of projections in order to formulate the mixed method. However, their numerical examples were limited to one dimensional problems and further, the convergence of the strain energy was not explored. 1.2 Purpose and Scope The purpose of this thesis is to extend the work of previous in-vestigators in the f i e l d of mixed fi n i t e element method and to define the mathematical, framework in which the procedure can be used most advantageously. A new energy product and the associated energy norm are defined which are then used in establishing the energy convergence of mixed methods in continuum mechanics. The completeness requirements for the approximating solutions are also established. The theory is then extended to i t s applica-tion in the f i n i t e element analysis. The continuity requirement for displace-ments and stresses and their influence on error in the energy product is also discussed. The forced and natural boundary conditions seem to be inter-changeable depending on the boundary integrals and how these are treated during the derivation of the mixed variational principle. In the case of the mixed Galerkin method, which provides exactly the same results as the mixed variational principle for self-adjoint linear boundary-value problems, the forced and natural boundary conditions w i l l depend on how the boundary residuals are accounted for, i.e. either through the displacement boundary residual or the stress boundary residual. In the mixed f i n i t e element formulation, some approximations for the dependent variables can lead to mechanisms in general and self-equilibrat-ing systems for incompressible cases. In order to explain this, the eigen-values and the composition of the eigenvectors for the linear elasticity plane 4 stress and linear part of the Navier-Stokes equations are studied for various combinations of polynomial approximations of the dependent variables. These are stress and displacement for plane stress, and stress, pressure and velocity for the incompressible, viscous flow. The latter is very similar to the plane strain e l a s t i c i t y problem for incompressible material. A number of numerical examples for one and two dimensional problems are presented. In the one dimensional case, the fourth order beam equation is f i r s t decomposed into two second order equations and then into four f i r s t order equations. The different combinations of the forced boundary conditions are also demonstrated for the beam problem. The plane stress mixed f i n i t e element is formulated, using linear stresses and displacements, and used in the analysis of a square plate with parabolically varying end loads, a canti-lever with parabolic end load and a rectangular plate with symmetric edge cracks. Also, the plane strain problem of stresses around a circular hole, both isotropic and orthotropic cases, is analysed. The latter orthotropic case only requires slight adjustments to the element matrix of plane stress case. Finally, the stress intensity factor for both symmetric edge cracks and a central crack is determined from the energy release rates. 1.3 Limitations In the development, analysis and applications of the mixed methods carried out in this thesis, the linear boundary-value problems of self-adjoint operators are considered. These cover the largest class of problems in continuum and structural mechanics. It is hoped that the theory developed here w i l l provide better grounds for extension to problems involving non-self-adjoint non-linear operators. 5 CHAPTER 2 MATHEMATICAL PRELIMINARIES In this chapter basic concepts, definitions and theorems are pre-sented to recapitulate some of the mathematical fundamentals which shall be used later in the development of the theory of mixed methods. For further cl a r i f i c a t i o n and proofs of the theorems, the texts by Lorch [17], Mikhlin [18], Balakrishnan [43] and Hellwig [11] are hereby referred to and shall not be repeated again. 2.1 Basic Concepts The f i r s t step in predicting the behaviour of a physical system by mathematical analysis is to idealize the system so as to obtain a mathematical model. In many cases i t is a differential equation with some boundary condi-tions. The purpose of the fi n i t e element method is to provide an approximate solution to the differential equation and boundary conditions that w i l l con-verge to the right solution in some sense, e.g. energy, mean square, etc. The solution of a given equation is obtained by finding that function which, when acted upon by the given operator, yields a known function. Attention w i l l be restricted to those functions that are square integrable over a given domain ft, i.e. functions u^ such that the Lebesque integral j \ u.u.dft < °°. (2.1) 'ft i i The functions considered w i l l , in general, be vector valued and u^ is the i* " * 1 component of a column vector u, i.e. T U = < U i U o . . . U > . 1 n The class of square integrable functions over ft, denoted by L2(ft), constitutes a vector space over the f i e l d of real numbers. 6 In order to compare different approximations of a given equation, i t is necessary to introduce a norm. This can be accomplished by f i r s t defining an inner product. (Only real vector spaces are considered.) Definition 2.1.1 A real vector space H is called an inner product space (also pre-Hilbert space) i f there is defined a real-valued function of pairs of vectors u and v in H denoted by (u,v) which satisfies the following conditions: (i) (u1+u2,y) = (u l sy) + (u 2,y) ( i i ) (ciu,v) = ct(u,v); a = a constant ( i i i ) (u,y) = (y,u) (iv) (u,u) > 0, where the equality holds i f , and only i f , u = 0. Definition 2.1.2 An inner product space in which every Cauchy sequence is a con-vergent sequence is said to be complete. A complete inner product space is called a Hilbert space. If the L2 norm of u is defined as || uB = /(u~u7 (2.2) then the following theorem holds: Theorem 2.1.1 For any vectors u, y and w E H (i) ||u|| > 0; and ||u|| = 0, i f , and only i f , u = 0; ( i i ) I ctu|| = | ct | I u|| ; a = constant; ( i i i ) |(u,y)| < ||u|| ||y|| (Schwarz inequality); (iv) I u+y|| < ||u|| + I y| (triangle inequality). The difference between two approximating functions can be characterized by the norm of their difference, i.e. ||u-y||. 7 2.2 The Limit Process in Hilbert Spaces Definition 2.2.1 Let there be given a Hilbert space H and let {u^} be a sequence of i t s elements. This sequence converges, or tends to u i f u is an element of H and lim II u -u|| = 0. (2.3) " -n -" n-x» The element u i s called the limit of the sequence ^Hn^» a n c* * s written as u •> u, lim u = u. -n -n n-*x> The following are examples of the above definition: (a) In the space L 2(ft) the convergence u " H i means that lim f |u -u|2dft = 0 ' il -n - ' i.e. that u^ converges to u i n the mean. (b) Convergence L2(ft)=W2(ft) denotes convergence in the mean of the functions of the sequence (u } a n a a l l their derivatives up to order r inclusive to a limiting function and a l l i t s corresponding derivatives. Theorem 2.2.1 If u and v e L o(fl) converge to u and v, respectively, in the L o -n -n ^ •> norm, then (u ,v ) -* (u,v). (2.4) -n -n - - • Corollary 1. If u -m then (u ,v)->-(u,v) > -n - -n - - -Corollary 2. If u -Ma then ||uJHuJ. If {u nJ is a specified sequence in a Hilbert space H which converges to an element u, then by definition 2.2.1 limlu -u|=0. This means that for any specified e>0 i t is possible to find a number no( e) such that for n>r)o(e), l l " r , - i j l < e -8 Let k>no("|") and n>n0(-|), t h e n H-k"-'l < f a n d I-n"-H < 2* Now, the norm of the difference of u, and u can be estimated by -k -n the triangular inequality H-k"-J = I ( V- ) ' ( V- ) I I < " v - l + " v - l i < e-The last equation by virtue of the arbitrariness of e implies that lim 1 u. -u || = 0. (2.5) . "-k -n k-x>° Thus i f un"*"u> then equation (2.5) is necessarily satisfied. Since, in the space L 2(ft), equation (2.5) necessitates the existence of a limiting element u, the converse also holds. 2.3 Orthogonality and Orthonormal Basis An important notion i n any inner product space is that of ortho-gonality. Two vector functions u and v are said to be orthogonal i f their inner product is zero, i.e. (u,v) = 0. Definition 2.3.1 The orthogonal complement of a set 0 in a Hilbert space is the set of a l l elements orthogonal to every element in 0 . It is denoted by 0 X . Definition 2.3.2 An orthonormal set is one in which any two elements are orthogonal to each other, and each element of the set is of unit norm. Theorem 2.3.1 Every non-trivial Hilbert space has an orthonormal basis. (For proof, see Balakrishnan [A3]). 9 Definition 2.3.3 A set is said to be dense in a Hilbert space H, i f i t s closure is equal to H. Definition 2.3.4 A Hilbert space is said to be separable i f i t has a countable dense set. Definition 2.3.5 A sequence of functions u^u^, . . . ,u^ is said to be complete in L 2(fi) i f for a function u with a f i n i t e L 2 norm and any e>0 i t is possible to find a natural number N and constants cti,ct2, . . . ,ct^ such that N || u- (aiui+ct2u2+ . . . +a^u^)| = || u- £ ct.u.|| < e. (2.6) i=l Theorem 2.3.2 If an orthonormal sequence of functions u^uj, . . . ,u N is complete, then N the Fourier series E a.u. of some function u with a f i n i t e L? norm converges i = l i i to this function in L 2 norm. The coefficients a^ are given by a. = (u,u.). (2.7) x 1 In this case there occurs the so-called Parseval equation; OO 00 ||u||2 = X a? = I (u,u.) 2. (2.8) i=l 1 i=l 1 2.4 Subspaces and Projections Consider the space L 2(fi) of scalar or v e c t o r functions which are defined in a certain f i n i t e domain Q, and have a f i n i t e L 2 norm. Select some linear set of functions belonging to this class and add a l l of i t s limiting elements to this set (i.e., functions which are the limits in the mean of a sequence of functions belonging to the given linear set). Such linear sets 10 of functions are called closed subspaces of the basic space L 2(fi) (e.g. Let fi be equivalent to the segment (0,2TT) of the x-axis. The set of functions whose Fourier series contains only sin(nx),n=l,2, . . . constitute a subspace. Obviously the functions of this subspace are characterized by the fact that they are orthogonal to the functions cos(nx),n=l,2, . . .). Sometimes the linear set leads to a subspace which coincides with L 2(fi). Thus adding a l l the limiting functions to the set of a l l polynomials (which is obviously linear) would constitute the class L 2(fi). Let there be given a separable Hilbert space H and one of it s sub-spaces Hj. Then from definition 2.3.4, there exists in Hi a complete f i n i t e or enumerable orthonormalized system ui,u 2, . . •>un' Take an arbitrary function u, which has a f i n i t e norm but does not necessarily belong to H^ , and construct its Fourier series in terms of the functions ui,u2 u . n OO E a u : a = (u,u ). (2.9) , n n n n n=l The sequence ( u n) m a y be composed of a f i n i t e number of terms. In such a case the sum in (2.9) w i l l contain only a f i n i t e number of terms. From theorem 2.3.2, the series converges in the L 2 norm. Let i t s sum be oo denoted by u, that i s , u = E a u . In fact, u is the limit as N-**> , n n n=l of the functions N E a u , n n n=l which belong to the subspace H l s since the sequence (u^le Hi. Therefore, u E H]. The function u is called the orthogonal projection, or simply the projection of the function u onto the subspace Hi. The difference u=u-u is orthogonal to the subspace Hi. 11 Therefore any function u e H can be represented as the sum of two terms, i.e. u=u+u, of which the f i r s t is the projection of the given function u onto the given subspace Hi and the second one is orthogonal to this subspace. The projection of u onto H} does not change when one system of functions {u }, which is complete and orthonormalized, is replaced by some other complete, n orthonormal system {w^ } in Hj. A projection possesses an extremal property. If u e L 2(fi) and u is the projection of the function u onto a subspace , then the norm of the difference u=u-v where v is an arbitrary function from Hi, becomes a minimum for v=u. In fact || u-vl| 2 = ||u+(u-v)||2 = ||u||2 + ||u-v|2 (2.10) since (u-v) e Hi and thus (u,u-v)=0. It becomes obvious now that || u-v||2 becomes a minimum at v=u; this minimum equals ||u||2. Consider a l l possible functions u e L 2(fi) and their projections u onto a given subspace Hi. The differences u=u-u constitute some new subspace H2. Every function from H 2 is orthogonal to every function from Hi and i t is said that the subspaces Hi and H 2 are orthogonal. Further, any function u e L 2(fi) expands into a sum u=u+u, where u e Hi and u e H2. This fact is usually formulated by saying that L 2(fi) is an orthogonal sum of the subspaces Hi and H2, and each of the subspaces Hi and H 2 is the orthogonal complement of'the other subspace. The properties mentioned about L 2(fi) in this section apply without change to any arbitrary Hilbert space. If a given Hilbert space is the orthogonal sum of the subspaces Hi and H 2 then i t can be written as H = Hi © H2 or Hi = H0 H2; H 2 = H0 Hj. 12 2.5 Application to Linear Elasticity A typical differential equation in linear e l a s t i c i t y is of the form Aw = f in tt (2.11) where w and f are members of some Hilbert space H(ft) and f is known. A solution w sought for the differential equation (2.11) w i l l also be required to satisfy certain boundary conditions. A class of functions that satisfy exactly a l l the boundary conditions of the problem and possess sufficient continuity so that Aw i s defined is known as the f i e l d of definition of the operator A and is denoted by D^ . In general, A is considered defined for a dense set of H(fi). For a large majority of problems in linear e l a s t i c i t y , the operator A has the following properties: (i) linear; A(au+gv) = aAu + 6Av; where a, B are constants and u, v e D ; (i i ) symmetric; (Au,v) = (u,Av) for a l l u, v e D^ ; ( i i i ) positive definite; (Au,u) > 0 for a l l u e D , where the equality holds i f , and only i f , u = 0; (iv) positive bounded below; (Au,u) > r 2(u,u) for a l l u e D^ , where r is a positive constant. Let u and v be two functions in D,. One convenient measure of - . A closeness of these functions is the square root of the energy of their difference, i.e. the energy norm of u-v. This is now defined. Definition 2.5.1 For a positive definite operator A, the energy product of functions u and v e D^ over tt is given by [u,v] = (Au,v) = / vTAudft. (2.12) 13 The energy norm ju| i s defined as |u| = /[u,v]. (2.13) The energy product and the energy norm defined above s a t i s f y the properties of an inner product and a norm presented i n d e f i n i t i o n 2.1.1 and theorem 2.1.1, re s p e c t i v e l y . Thus a measure of closeness of the functions u,v e D. i s - - A |u-v| = /[u-y,u-y] = /(A(u-v),u-v). (2.14) The space D may be incomplete with respect to the energy norm, i . e . not a l l Cauchy sequences i n D converge to a function i n D i n the energy A A norm. I f so, a function u i s defined to be a member of the space by the following l i m i t i n g process: l^n'-l 0 a s n °°» (2.15) whereu i s a t y p i c a l member of a sequence {u } e D.. Thus, the space D, -n -n A A i s completed by including a l l the l i m i t i n g elements. The completed space so obtained i s a H i l b e r t space and i s denoted by H to emphasize i t s dependence upon the operator A. The d e f i n i t i o n of the energy product i n equation (2.12) can be extended to a l l functions i n H.: A X [u,v], = l i m (Au ,v ) = l i m /_ v Au dfi; u ,v e D.. (2.16) - - A — n -n ; f i - n — n -n -n A n-x» n-x" The completeness d e f i n i t i o n 2.3.5 and theorem 2.3.2 can now be rewritten as: D e f i n i t i o n 2.5.2 A sequence of functions u i , U 2 , . . .,u^ i s said to be complete i n H^ i f f o r a function u with a f i n i t e energy norm and any e>0 i t i s possible to fi n d a natural number N and constants cti ,o i2> • • • s u c h that N Ju-(a 1ui+a 2u 2+ . . . +a Nu N)| = Ju- E a i u i | < e* (2.17) 14 Theorem 2.5.1 If an orthonormal sequence of functions u l 5 u 2 , . . .,u^ is complete, N then the Fourier series I a - t u ^ °f some function u E H with a f i n i t e energy norm converges to this function in the energy norm. The coefficients a^ are given by a± = [u,u iJ. (2.18) In this case there occurs the so-called Parseval equation |u|2 = Z [u,u ] 2 . (2.19) 1=1 1 15 CHAPTER 3 MIXED METHODS In this chapter, a brief introduction to the theory of mixed methods is presented f i r s t . Then the definitions and theorems given in Chapter 2 are used to establish the convergence of these mixed methods. Since, these methods involve indefinite operators, their fields of definition have to be restricted in order to define their energy products and energy norms. The procedure for doing this is illustrated for a typical indefinite operator, and the resulting energy norm is then used to prove convergence in energy and to establish an error estimate. In the process the completeness requirements are also laid out. A remarkably large class of problems in mathematical physics involves equations of the form -Au-f = 0 (3.1a) Bu-gj = 0 on 3i?! and B*Tu-g2 = 0 on dR2, (3.1b) n where u=u(x) is a function defined on a bounded region R of E ; dR is the smooth boundary of R. The operator A is assumed to have the following properties: (i) A i s factorable, i.e. A=T*T, and ( i i ) A is positive definite. B and B* are linear boundary operators. Thus equation (3.1a) can be written as -T*T-f = 0. (3.2) Here T is a linear operator whose domain D^ is in a Hilbert space U and i t s range in another Hilbert space V. The operator T* is the formal adjoint of T; i t s domain D_. is in V and i t s range is in U. If the boundary conditions in (3.1b) are homogeneous then, as the operator T* is the formal adjoint of T, the following Green's formula holds: (Tu,v)^ = (T*v,u) y {for every u e U and v e v} (3.3) where ( , ) and ( , )(/. represent inner products in the spaces U and V, respectively Through equations of the form (3.2), by using direct integral methods, a collection of variational principles can be developed, Oden [21], 16 which generalize the classical ones existing in the theory of linear elasticity, In many cases, bounds on the solution or some other quantity of interest can be obtained even i f the problem is not expl i c i t l y solvable. 3.1 Mixed Variational Principle A mixed variational principle associated with equation (3.2) can be developed by splitting this equation into canonical form, equivalent to the following pair of equations: Tu = v in R; and B(u) = g^ on dRi (3.4) T*v = -f(u) in R; and B*(v) = g 2 on dR2 where dRi+dR2=^R, and the linear boundary operator B (and B*) depends upon T. For the non-homogeneous boundary condition of (3.4), T* is the formal adjoint of T i f i t satisfies the generalized Green's formula (v,Tu)^= (u,T*v) y+ ( v , B u ) m (3.5) where (v,Bu)^.^ denotes an inner product in the space V associated with the boundary terms. The only difference between (3.3) and (3.5) is the addition of the boundary term to the right, hand side. The boundary operator B* is the adjoint of B in the sense that (v'BuW = (U'B*VW <3-5a> In certain cases i t is also possible to write equations (3.4) in a generalized Hamiltonian form. Assume that there exists a Gateaux differentiable, (Balakrishnan [43])} bilinear functi onal H(u,v), called the generalized Hamiltonian, whose total Gateaux differential is of the form 6H(u,v;n,w) = (-f(u).n)^ + (v,w)^ (3.6) where n is the variation of u and w that of v. In equation (3.6), H(u,v) is assumed to have the property that i t s partial Gateaux derivatives with 17 respect to u and v, 6H/6u and 6H/<5v, respectively, coincide with the right hand sides of the canonical equations (3.4), i.e. (3.7) Then the generalized Hamiltonian forms are analogous to those of analytical dynamics, Tu = T*v = 6H <5v 6H (3.8) with B(u) = gi on dRi and B*(v) = g 2 on 8i?2. The direct integral method can now be employed to derive a func-tional for (3.4) or (3.8). Let W denote the tensor product space, W = UxV. Then elements A of If are the ordered pairs T A = <u v> ; u e U, v e V. The equations (3.4) can be put into a compact form (Oden [21]) P(A) - r = 0 where the matrix operator P and T are: P = 0 T* T 0 -f (u)l and T = (3.9) (3.10) (3.11) Similarly the g e n e r a l boundary conditions can be written symbolically as (3.12) 5(A) - ldR = 0 on dR where B -0 B* -B 0 • T ' -dR §2| 1-81. (3.13) Denote the inner product of elements A e w=uxv as (A!,A2) = (u 1,u 2) f /+ (v!,v 2)^, (3.14) 18 T T and A =<u^ v^> , A2=<u2 v2> are the ordered pairs. A functional can now be constructed by integrating the inner product of the residual of equations (3.10) and (3.12) and a variation of an arbitrary element dl\ in W, i.e. J(u,v) = J(A) = fR {[P(A)-r] + [B(A)-T9i?]}6Adi< = l/2(Tu,v) y + 1/2(1^^)^ + F(u) - l/2(v,v) u + l/2(B*v,u) y 3 R 2 - l/2(Bu,v) 7 3 i ? 2 - l/ 2 ( g 2 , u ) [ / 9 i ? 2 + 1 / 2 ( 8 1 , v ) ^ (3.15) Here • F(u) = (f,u) y . (3.16) In computing the boundary terms i n (3.15), i t is important to realize that both B and B* depend upon T. Moreover B(u(x))=0 i f x e dR2 and B*(v(x))=0 i f x e dR± for the boundary conditions given in (3.4). In the formulation of the variational statement (3.15) inclusion of the boundary terms is analogous to the boundary residual concept presented by Finlayson and Scriven [8] and to the principle of virtual work. Theorem 3.1.1 The functional J(u,v) of (3.15) assumes a stationary value at the point (u,v) which satisfies the canonical pair (3.4), where the operator P is defined for some dense set of the space W. Proof : - - w Let the varied solution be u=u+cxn and v=v+6w, where u,n and v,w e D^ . Substitute i t into the expression (3.15) for J(u,v), 19 JCu+an,v+Bw) = 1/2CT(u+an),v+6w) + 1/2(T*(v+Sw),u+an)y + F(u+an) - 1/2(v+Bw,v+Bw)^ + 1/2(B*(v+Bw) , u + a n ) - 1/2(B(u+an) ,v+Bw)y " ^ 2 , ^ ) ^ 4 - ( g l , v + 3 w ) u 9 i ? i . (3.17) Then for a stationary point of J(u+an,v+Bw); 6J(u,v;n ,w) = lim 3J (u+an ,v) + Um 3J(u,v+gw) = 0 a+0 3a B+0 3 6 = [ l / 2 ( T u , w ) K + l/2(T*v ,n) y + (f(u),n) + 1/2 ( B * v , n ) u 3 " 1 / 2 ^ > " \ , R l ~ + [l / 2(Tn,w) v 3 R i + 1 / 2 ( ^ , 0 ) ^ " <v,w)y + l / 2 ( B ^ , n ) ^ 2 " l / 2 ( B u , w ) ^ i + ( g i . w ) ^ ] - [ ( T u , w ) ^ - (v.w), - ( B u , w ) ^ i + (gi,w) v 3 R i] + [ ( T ^ . r O ^ + ( B * v , n ) u 3 i ? 2 - ( g 2 > n ) u 3 i ? 2 ] = t(T*v+f(u),n)y + ( B ^ - g z . n ) ^ ] + [(Tu-v , w ) y - (B5-gl,w) ] = 0 (3.18) Since the variations n and w are arbitrary, then from Lagrange's lemma T*v + f(u) =0; B*v = g 2 on dR2 Tu - v = 0; Bu = gj on dRx. Hence the equations obtained are the same as equations (3.4), and therefore the solution at the stationary point (u,v) does satisfy the canonical pair (3.4). Boundary Conditions Rewriting equation (3.5) by spli t t i n g the boundary inner product (Tu,v) y = (u,T*v) y + ( B ^ u ) ^ + (Bu,v) K 3 K i (3.19) the variational statement in (3.15) can be rewritten as J(u,v) = (Tu,v) v - l/2(y,v)v + F(u) - (Bu-gi,v) - ( 8 2 , " ) ^ (3.20a) 20 or J(u,v) = {T*vlu)u - l/2(v,v) 7 + F(u) + (B*v-g 2.u)^ ~ CBu-8l ,v>^3jRi + CBOv.u)^. (3.20b) The equivalent forms of J(u,v) in equations (3.20a) and (3.20b) provide some f l e x i b i l i t y as to how boundary conditions can be incorporated. In (3.20a) the boundary integrals are extracted from the second of the canonical pair (3.4) (the equilibrium equation), and in (3.20b), from the f i r s t of (3.4) (the constitutive-compatibility combined equation). This idea shall be explained later in detail (Chapter 5 on boundary conditions). It is worthwhile noting that under the assumption Tu=v being exactly satisfied, the functional (Ju,v) in (3.20a) reduces to I(u) = l/2(Tu,Tu) + F(u) - (g2>u) 3 i ? 2 (3-21) This is the functional for the principle of minimum potential energy, the completeness, energy convergence and bounds for which are well established, Mikhlin [18]. 3.2 Projection Operators The orthogonal projection, or projection of a function into subspaces was discussed in Section 2.4. The method of orthogonal projections can also be used for estimating errors in approximations of linear boundary value problems, Mikhlin [18]. The idea of projections of linear operators has been exploited by Reddy and Oden [30] in establishing convergence rates for the basic variables involved in the mixed methods. In general, there exists a number of possible projections for a given operator, e.g., T*T. The four important projections are cited here: 21 (i) Primal Projection (equivalent to the conventional Ritz-Galerkin approximation) ( i i ) Dual Projection (similar to the Equilibrium Finite Element method) ( i i i ) Primal-Dual Projection -\ > (both together lead to the mixed formulation) (iv) Dual-Primal Projection J Consider two linear vector spaces, U and V, defined over the same f i e l d , and let {i> }, (m=l,2, . . . M) denote a set of M linearly independent m elements in U and {ip }, (n=l,2, . . . N) a set of N linearly independent n elements in V. The sets {cb } and {ib } define an M-dimensional subspace $ e u and m n M an N-dimensional subspace ¥ e V, respectively. The Gram matrices associated with the subspaces and 4* are M N G i j = ( V * j V H i j ° ( W r ( 3 ' 2 2 ) which are not singular since {cf> } and {tj; } are linearly independent and the m n biorthogonal bases can be computed directly (j)1 = G ? U • 5 ' <J>1 = HTH. . (3.23) i j J i j 3 From equations (3.22) and (3.23), i t can be seen that the b i -orthoganality conditions are O i ' ^ U = &3±; (^ i,^ J) 1 / = 63±. (6^ =1, when i=j and zero otherwise). (3.24) Note that there is no relation between the spaces U and V, and the biorthogonal bases i n $ w and are completely independent, unless some additional infor-M N mation i s provided. Definition 3.2.1 The orthogonal projection operators ]T:L7->-$m and P:!/"-*^ can now be defined in the following sense: i f u is an arbitrary element in U and v is an arbitrary element in V, the projection u of u into subspace $ M and the projection v of v into subspace V are of the form 22 M N H(u) = u = Z a 1^ and P (v) = v = £ b 1^. (3.25) i=l i=l 1 where a 1 = (ujtj)1)^ and = (v,^ "*") . (3.26) Now let U be a Hilbert space consisting of functions u defined on a compact, convex subset R of E n with a smooth boundary <$R, and let T be a bounded linear operator mapping U into another Hilbert space V, and T*:V •> U i t s formal adjoint satisfying equation (3.5) while B* and B are their boun-dary operators which satisfy equation (3.5a). Consider the cases in which S> e D , the domain of T, and 4* e D„,.. In general, T($ ) is not a subspace M l N 1 * M of YN, and T*(4'N) is not a subspace of $ . Operators T and T* can be approximated by projecting T($^) into 4^ and T*^ ) into $ . This projection process leads to rectangular matrices of the following type: PT($ M): PT(ci_L) = m^J (Dual-Primal) (3.27) IIT*(YN): HT*^±) = njcbj . (Primal-Dual) (3.28) Here PiV-vf^ and n:L7->-$^ are projection operators defined by bases cb1 and i p 1 , and m i = ( T < ) )i ' * j V ; n j = ( T * v * i V ( 3 , 2 9 ) Similar projections can also be obtained for B and B*. 3.3 Mixed Galerkin Method Consider the problem of equation (3.4) where f is a function of spatial coordinates x only. That i s , T*v + f = 0 (3.30) TU._ v = 0. (3.31) Choose linearly independent sequences of coordinate function <(>i >^ 2 > • • • ^ e L7° for approximating u and ty\,ty2> • • • ^ e V^ ^ o r V o v e r fc^e s a m e 23 domain R which are continuously differentiable i n the closed domain R=R+dR as many times as required for the specified problem, and satisfy a l l the boundary conditions of the problem. Then the functions M "M = i=l (3.32) N v M = Z b.t|». N . 1 i r i 1=1 where a^ and b^ are arbitrary constants, also satisfy a l l the boundary conditions of the problem. In the mixed Galerkin method, the coefficients a. and b. are deter-mined from the requirement that the residual of equation (3.30), is made orthogonal to the coordinate functions for u^ and the residual of equation (3. i s made orthogonal to the coordinate functions for over R. This then leads to the following system of algebraic equations: N I b ( T*iK , c b . ) = -(f ,Q> ) ; 1=1,2,3, . . . M • (3.33) i =1 1 M N E a.(T(j) ,^ ) - Z b.Op. ,ifO = 0; j=l J j 1 k=l J J 1 1=1,2,3 . . . N (3.34) which consists of (M+N) equations with (M+N) unknowns (a^, i=l,2, . . . M; b , j=l,2, . . . N). For e q u a t i o n s (3.30) and (3.31), the operator T and i t s adjoint T* are defined for sets which are dense in some separable Hilbert spaces U® and V®, respectively, and the sequences {<f>^} e and {^} e D^ ,A are complete. The derivation of equations (3.33) and (3.34) is the same as the estimation of T and T* by projection operators i n the previous 24 section, equations (3.27) to (3.29), except here the bases are not b i -orthogonal. Furthermore, e q u a t i o n s (3.33) and (3.34) would be identi-cal to those derived from the functional J(u,v) of equation (3.15), since T* is the formal adjoint of T. Before the completeness and the convergence of mixed methods are presented, certain definitions and concepts have to be introduced. 3.4 Concepts and Definitions The differential equations and boundary conditions in equation (3.4) can be put into matrix forms somewhat similar to equations (3.10) and (3.12) as 0 T* T -1 in R and 0. B* B 0 on DR or AA = P in R and BA = g on c where 0 T*~ "o B* A = and B = _T -1 B 0 which are linear operators, and T where A forms ordered pairs <u v> just like in equation (3.9) (3.35) (3.36) (3.36a) Definition 3.4.1 The class of functions that satisfy a l l the boundary conditions of the problem (3.35) and possess sufficient continuity to make the evalua-tion of AA possible is known as the f i e l d of definition of A and is denoted 25 by W^. For example i f A=T*T were a fourth order differential operator, then T and T* would be second order differential operators and the functions u and v in W would have continuous second derivatives in R and would satisfy the boundary conditions in (3.35) on dR. Definition 3.4.2 The operator A is called symmetric i f for any elements A and A from the f i e l d of definition of this operator W , the following identity is valid: (AA,A) = (A >AA) (3.37) Theorem 3.4.1 Operator A as defined in (3.36a) i s a symmetric operator. Proof: For A and A e W^, form the product < A M V = < T * v ' a > t 7 + (Tu-v,v)^ A = (T*v,u) y + (Tu,v)v - (v,v)^. Using equation (3.5) A Since A and A satisfy the same boundary conditions, then (AA,A) = (T*v,u) + (Tu-v,v) = (A,AA) v — - WA v ' 'u V - -VI, A A Hence, as long as T* is the formal adjoint of T, the operator A is symmetric, 26 In order to define a norm in the product space IV , an inner product has to be introduced. In the following i t is assumed that the approximations u^ and v^ given in equations (3.32) are to be used for u and v, respectively. Then when the Mixed Galerkin Method is applied to equations (3.35), the second equation yields: ( v V v = ^ v ^ V j = 1 ' 2 , • • " N ( 3- 3 8 ) M where This equation may be best understood by thinking of Tu^ as a known function f. Then i f the are orthonormal, equation (3.38) reduces to b_.= (f , i(Oy,, i.e. the Fourier coefficients. Further, i f the sequence of functions is complete and Tu^ e V, then v^ can be made arbitrarily close to Tu^, Lorch [17]. However, in formulating an energy product for the mixed method, another question arises and that i s : is Tu^=0 when v^O ? This may be answered by noting that for vN=0 (which implies that the coefficients b^ .=0) , equation (3.38) yields (Tu^ijO^ = 0, for each j . (3.39) Then, i f the sequence of functions is complete and i f Tu^ is restricted to be in V, Tu^=0 necessarily, Lorch [17], This is now used to define a w restricted f i e l d of definition D. for the operator A. A -Definition 3.4.3 w The restricted f i e l d of definition D^ for the operator A is defined as the product space U * V where the restricted spaces U and V are subspaces of U and V, respectively. A sequence of functions {((J }^ E U 27 M are used for the approximation of u, i.e. u^ . = E a^-^- Then the space V may be defined somewhat arbitrarily provided i t includes Tu^. That i s , there is a minimum requirement that TUL^. e V. Thereafter, any sequence of N functions {ib.A used for axproximating v, i.e. v„ = E b .ib. , are assumed to N N j=l J J be complete and to belong to V, i.e. £ V. Conversely, any sequence { i k , } e V restricts the choice of {d>w} e U. Such restricted sequences {d>w} then N M M constitute the restricted subspace U. Furthermore equations (3.38) are always to be enforced and therefore < v V ^ - ( T v V ^ ( 3 - 4 0 ) w where v^ and u^ now belong to the restricted space D^ . Now the energy product can be defined for this restricted space. Definition 3.4.4 The energy product is defined as [Aj,A 2] = (AAj, A 2) = / <u2 v2> 0 T* T -1 (3.41) (T*vi,u 2)~+ ( T U l , v 2 ) ~ - (v 1,v 2)~ where A^ = <U], v^> and A2=<u2 v2> e D . (Note: If Tu=v is satisfied exactly i t can be shown that [A,A]=(Tu,Tu), i.e. the above energy product is twice the strain energy). The energy product so defined has to satisfy the properties of an inner product in definition 2.1.1. The properties ( i ) , ( i i ) and ( i i i ) 2 8 follow from linearity and symmetry of the operator A whereas the property (iv), which follows from positivity of the inner product, is not as obvious. This is proved here. The energy product in (3.41) can be rewritten as [A,A] = (T*v,u)~ + (Tu,v)~ - (v,v)~; {A e = U*V). (3.41a) Assuming homogeneous boundary conditions, therefore from (3.3) (T*v,u)~ = (Tu,v)~ and substitution into (3.41a) yields [A,A] = 2(Tu,v)~ - (v,v)~. (3.41b) w Since A e D and equation (3.40) is satisfied, replacing (Tu,v)~ by (v,v)~ A V V in (3.41b) leads to [A,AJ = (v,v)~. ( 3 . 4 2 ) The right hand side of (3.42) is always greater than zero and equals zero i f , and only i f , v=0. Further, for any A in the restricted space D^ , i f vEO then also Tu=0 and i t follows from the homogeneous boundary conditions that u=0. Hence, the positive definiteness of the mixed operator A, when w every A is chosen from the restricted space D , can be ascertained as [A,A] >. 0, where the equality holds i f , and only i f , A = 0. (3.43) Further i t is bounded below as [A,A] > v2j|v||2 (3.43a) where v is a positive constant. This then proves that the property (iv) 29 of an inner product i n definition 2.1.1 also holds. The energy norm can now be defined as |A| = /[ATAT; (A e D"} and the following theorem holds. Theorem 3.4.2 If the energy norm is defined as |AJ = /[A,A] and the conditions of the energy product in definition 3.4.3 are satisfied, i.e. A is chosen from w the restricted space D , then (i) |A| > 0, the equality holds i f , and only i f A = 0; (i i ) |aA| = |a| JA|, a = constant; ( i i i ) [ [AlA 2]| < |A x j |A 2| (Schwarz inequality); (iv) [Ax + A2| < |Ai( + |A2j (triangle inequality). Proof ( i ) : This property of the energy norm follows from the positive definiteness of the energy product. Proof ( i i ) : From linearity of the operator A and substituting Aj-A^aA in equation (3.38), where a is constant [aA,aA] = a 2[A,A]. By taking the square root of both sides |aA| = /[aA,aA] = |a| |A|. Proof ( i i i ) : Consider an arbitrary positive real number X and construct the 30 non-zero function A j + X A 2 . The energy product of this function, by using linearity and symmetry of the operator A yields [ A i + A A 2 , A i + A A 2 ] = U i , A 2 ] + 2X [ A l 5 A 2 ] + A 2 [ A 2 , A 2 ] . By virtue of property (i) [A2+AA2,Aj+ A A 2 ] >0. (3.44) Therefore [ A L A J + 2X [Ai,hz] + X2 [ A 2 , A 2 ] >Q. (3.45) The l e f t hand side of (3.45) is quadratic in X, and a l l i t s values are non-negative. Hence i t follows that i t s discriminant i s less than or equal to zero [ A ^ A o J 2 - [ A i . A i ] [ A 2 , A 2 ] « 0 or [ A l 5 A 2 ] 2 < [ A i , A i ] [ A 2 , A 2 J . Taking the square root of both sides yields the required inequality | [ A l 5 A 2 ] | < IA!J |A 2|. Proof ( i v ) : If A=l i n proof ( i i i ) , the following results: [Ai+A^Ai+As] = [ A L A J ] + 2 [ A l 5 A 2 ] + [ A 2 , A 2 ] or | A l + A 2 | 2 < | A l | 2 + 2 | [ A 1 , A 2 ] | + | A 2 | 2 . Substituting the Schwarz inequality | [ A 1 , A 2 ] | <. \h1\ J A 2 J ; | A 1 + A 2 | 2 « | A i | 2 + 2|A!.| |A 2| + | A 2 | 2 or | A l + A 2 | 2 < ( | A i | + | A 2 | ) 2 . 31 Taking the square root of both sides gives the required inequality |Al+A2| < |Ai| + |A2|. Similarly i t can be shown that lAi-A 2| « jAj | + |A2|. 3.5 Convergence of the Mixed Method The mixed variational principle of equation (3.15) for homogeneous boundary conditions, i.e., gi=g2=0 in equations (3.4), in lieu of definition 3.4.3 of the energy product, c a n be rewritten in the following form: F(A) = (AA.A) - 2(p T,A) (3.46) where A e From theorem 3.1.1, let AQ=(uo,vg) satisfy the canonical pair (3.4) for which the functional F(AQ) assumes a stationary value, and has a f i n i t e value in general. Then the matrix form (3.36) becomes AA0 = p_. (3.47) Substituting for p in (3.46), F(A) = (AA,A) - 2(AA0,A) = [A,A] -2[A 0,A]. (3.48) Adding and subtracting [AQ.AQ] from the right hand side of (3.48) F(A) = [A,A] - 2[A0,A] + [A 0,A 0] - [A 0,A 0] which can be formally shown to be F(A) = [A-A0,A-A0] - [A 0,A 0]. (3.49) Now i f A=AQ the exact solution, then F(A 0) = -[A 0,A 0]. (3.50) Let d, some real number be the exact stationary value of the functional. Therefore - d = F(A 0) = -[A 0,A 0]. (3.51) 32 In general, the functional F(A) does not lead to a minimum problem. However, as a consequence of the conditions imposed on the energy product^ i t ' s positive definiteness provides the sufficient condition for the minimum. Therefore the following minimum functional theorem is presented. Theorem 3.5.1 Let A be a symmetric and indefinite operator as defined i n equation (3.35) and further that this equation has a solution. Then of a l l the values which are given to the quadratic functional F(A) = (AA.A) - 2(pT,A) (3.52) by a l l possible functions from the restricted f i e l d of definition of the operator A, the actual minimum occurs only for the solution of equations (3.35). Converse: W T If there exists in D, a function A=<u v> which also satisfies A the conditions of the definition 3.4.3 and gives the minimal value to the functional (3.52), then this function is also the solution of equation (3.35). Proof: W T T Assume that A,Ai e D and A=<u v> and Ai=<u n> . Set A-A0=Ai where AQ is the solution of equations (3.35). Thus A=AQ+AJ. Therefore from linearity and symmetry of the operator, F(A) = [A,A] - 2(pT,A) = (A(AQ+AI)JAQ+AJ) - 2(pT,A0+A1) = [A 0,A 0] - 2(p T,A 0) + 2(AA0-pT,A1) + [ A L A J ] = F(A 0) + 2(AA0-p,A1) + [ A L A J . 33 But AAo-p=0 by hypothesis. Hence F(A) = F(A 0) + [ A ^ A J . Now from positive definiteness of the energy product [Ai.AjJ > 0, where the equality holds i f , and only i f , A^O. Thus F(A) > F(A 0) with equality only valid i f A^EO or A=A0. Hence the functional attains i t s minimum value when A=A0 and from (3.49) min F(A) = F(A 0) = -[A 0,A 0] = -|A0|2. The converse follows from theorem 3.1.1. The functional T F(A) assumes a stationary value at A=<u v> which satisfies the equation (3.35) and further, this stationary value is a minimum when the conditions of the definition 3.4.3 are satisfied. Therefore the function AQ E D which gives the functional in (3.52) the minimal value is also the solution of equation (3.35). Theorem 3.5.2 The approximate solution An=<u v >^ E of equation (3.35) r -m m n A constitutes a minimizing sequence for the functional (3.52) provided that equation (3.35) has a solution with f i n i t e energy and that the conditions of definition 3.4.3 are satisfied. Proof: From theorem 3.5.1 let d be the minimum value of F(A), i.e. d = min F(A) = -|A0|2. n *•* — T Let the approximation A =<u v > be given by -m m n _ m u = E 6, u. m i=l 1 1 (3.53a) n v = E IJJ . v. n j - l ^ 34 n n W — — Therefore F(A )>d, since A e D, and u and v are l i n e a r combinations of -m -m A m n the sequences {cp^ } and { ^ l , r e s p e c t i v e l y , which enter into D^. The quantity d i s the exact lower bound of the function F(A). I f z i s an a r b i t r a r y small p o s i t i v e number, then, by d e f i n i t i o n of the exact lower bound, there e x i s t s a function A e such that d<F(A)<d+j. The task now reduces to choosing natural numbers m and n and constants Qi,uo, . . . , G and vi,vo, . . . , v such that j.» z > ' m ' n they s a t i s f y the in e q u a l i t y F ( A " ) - F(A) < f . (3.53b) -m - l Using equation (3.49) this i n e q u a l i t y reduces to ^-A 0>^-A 0] - [A-A0,A-A0] < |. (3.53c) Consider the l e f t hand side t ^ - A 0 . A ^ - A o 3 " [A-A0,A-A0] = | ^ - A 0 j 2 - |A-A0|2 = t | A ^ - A 0 | + |A-A0|} {|A*-A0| - IA-Aol>- (3.54) By the t r i a n g l e i n e q u a l i t y , property (iv) of theorem 2.4.2 I^ -Ao j « |A^-A| + |A-A0|. Therefore |A*-Ao| - | A - A 0 | < |A*-A|. Subs t i t u t i o n into (3.54) y i e l d s f A ^ - A o » A ^ - A 0 ] - [A-A0,A-A0] < { | / £ - A 0 | + | A - A „ | } | / £ - A | . Select m,n and c o e f f i c i e n t s ui,u?, . . . u ; v-i.vo, . . . v such that The value of k w i l l be chosen l a t e r . Also from the t r i a n g l e i n e q u a l i t y < 1-0 + 141 • 35 Therefore |?J-A 0| + |A-A0| « + | A | + 2|A0|. Now |A"| < |A| + f and [^-A0,A^-A0] - [A-A 0,A-A 0] < {2|A| + 21A0 | + -} p (3.54a) By taking k such that £{2|A| + 2|A0| + f ) <\ inequality (3.54a) reduces to (3.53), i.e. Therefore F(A N) - F(A) < -|. -m - z Hence, i t follows that d « F(A N) « F(A) + \\ < d + e. (3.55) -m — / Let A^ be the solution obtained by minimizing the functional F(A) in (3.46). Then d <s F(A N) ^ F(A N) (3.56a) -m -m or d « F(A N) < d + e. (3.56b) —m Thus by letting £->-0, the functional F(A N) converges to the exact value d, -m i.e. the sequence A =<u v > is minimizing. -m m n Rewriting inequality (3.56b) as 0 < F(A N) — d E , -m and using equations (3.49), (3.51) and positive definiteness of the energy product, this inequality reduces to 36 or Again as E->0, therefore A -A0 < E . -m - u jAn-A0 -v 0. •-m - u (3.57) But from equation (3.42) Therefore i v v °« 0 (3.58) Hence convergence in energy also implies the mean square convergence of the approxi-mating stress v to the exact value vn. m In order to show the mean convergence of the displacement, consider a free vibration problem for which the equation (3.35) takes the form 0 u = w 2 " l o" u (3.59) _T -1 V 0 0 V where co is the frequency of vibration. Alternatively, AA = ABA. (3.60) Here A is a symmetric operator as before, to2=A and 2 = Q 0 ' T Let Ag and AQ=<UQ vg> be an eigenvalue and i t s eigenfunction, respectively. Then AA0 = A0BA0, (3.61) and substituting into the energy product (3.41) yields (3.62) But Therefore [A 0,A 0] = A 0(BA 0,A 0) D^. A (BA0,A0)Dpv? = (u 0,u 0)~. A [A 0,A 0] 0 (uo,u0). An = (3.63) (3.64) 37 Since the energy product [A 0,AQ] represents strain energy, equation (3.64) is similar to the Rayleigh's quotient. Further, from the positive definiteness of [An,An] and observing that (ug,ug) is always a positive number, Ag is also a positive real number. W In general for any A e D , the equation (3.64) can be rewritten as ™" A [ A , A ] A = T r-^ ' (3.65) T Let A Q = < U Q vg> be the exact solution and consider the energy product of the difference A n - A n , where A n and A n satisfy the same boundary conditions. Then -m - u -m - u J from symmetry of the operator A Using equations (3.62), (3.64) and (3.65) ^m"-°'-m"-° ] = X(um'XXm)U + Xo( u0> uo) y " 2A 0(u 0,u m)~. Therefore from equation (3.42) ( v n - v 0 , v n - v 0 ) y = x ( u m > u m ) y + ^ o ( u 0 . u 0 ) ~ - 2A 0(u 0,u n)~. Because of the definition 3.4.3, the problem becomes that of a minimum. Therefore A > AQ and replacing A by A0 yields the following inequality: A ° ( V U ° ' V U O ) L ^ ( V V ° ' V V o ) y or a>0||um-u0|| < ||vn-v0||. Since tog is not zero, from (3.58) also II u -u0|| -> 0. (3.66) m u" Hence convergence in energy also implies the mean square convergence of the approximation to the exact value U Q of the displacement. The following theorem can now be stated for convergence of stress and displacement when the approximate solution converges in the energy sense. 38 Theorem 3.5.3 n - - T W If the approximations A =.<u v > e DA for the displacement and stress r -m m n A converge In the energy sense, i.e. |A^—AQ |->-0 as m and n ->• °°, then i t follows that the displacement u and stress v converge to the exact values un and vn , m n u u respectively, in the mean square sense. So far the energy convergence and the convergence of displacement and stress have been discussed in general terms. Now i n order to establish the completeness i n energy, i t is sufficient to reiterate the requirements in w definition 3.4.3. In defining the restricted space D , i t was required that Tu^ e V. Also the sequence of functions {ty } used for approximating v, i.e. n - ' • v = E li.v. must belong to V. Thus, when the Mixed Galerkin Method is applied n j=rj 3 to equations (3.35) using the approximations in (3.53a), the second equation yields n m E 0p . , i | O v . = E (ii. ,T<f>.)u. ; i=l,2,.. . ., n. (3.67) j=l 1 J J j=l 1 J . J Assuming e V are orthonormal for convenience, i.e. r 1 i f i = j (*, J = { 1 J U i f i j f j , equation (3.67) then reduces to m v. = E (I|J. ,T<f>.)u. = (Tu , i K ) . (3.67a) I . , I i i m i 3=1 J These are the Fourier coefficients for Tu with respect to the orthonormal se-in quence {ifj^ } e V. The equation (3.67a) can also be obtained by minimizing the L 2 norm II Tu -v II, Mikhlin [18], given by ^ 1 m n" IITu -v ||2 = L (Tu -v )2dR, (3.68) " m n" JR m n with respect to the unknown constants v.. This represents the mean square J convergence of v^ to Tu^ and plays an important role in defining completeness in energy of the mixed methods. Unless, for any Tu e V the set of functions m 39 {\pn} is complete for convergence in the mean square sense, the energy product in definition 3.4.4 is not positive definite. Therefore, the mean square con-n m vergence of the approximation of v, i.e. v = Z il.v. to Tu (u = E <b .u. : <b, n j=i 3 3 m m i = 1 x x i e u) acts as a prerequisite for completeness in energy. 3.6 Completeness The sequences of functions <)>i><)>2, • • • a n c^ ^ 1 J ^ 2 > • • • a r e said W to be complete in D i f for a function A with a f i n i t e energy norm and any £\>0, £ 2 > 0 , i t is possible to find positive integers M and N and constants cti , c t 2 , . . . ct^; &i,&2> ' * * s u c ^ t n a t the following inequalities are satisfied: (i) | | T U m - ( B I 1 | ; 1 + B 2 ^ 2 + • • • +eNVl'<e2 M T N - - T - _ N where A=<u v> , L,=<\x£ v > and u = I a.A., v = E B.i|>.. -M M N M xx N ±-i xx The f i r s t inequality (which comes from the second of equations (3.35), the constitutive equation) implies that the approximation of the stress N function v, i.e. v = £ 8.il>. for f i n i t e N should contain a l l the stress modes n i = 1 i r i that correspond to strain modes present in the strain Tu^ obtained from the N displacement approximation u,, and only then A w w i l l converge to A in energy. M -M -3.7 Estimation of Error in the Energy Product n W It was shown i n Section 3.5 that the mixed approximation A^ e of equation (3.53a) does converge in the energy sense. Therefore i t is desir-able to seek some estimate of error in the energy product (or strain energy). From equation (3.41) "0 T*| [A^-A0,An-A0] = / <u -u 0 v -v0> -m-u-m-u JR m u ™ u T -1 u -u 0 i m u 1 v n-v 0 J 'dR (u m-u 0,T*(v n-v 0))~ + (v n-v 0,T(u m-u 0))~ - ( V V 0 ' V V o ) ^ °- 6 9 ) 40 T where AQ=<UO vg> is the exact solution. Consider homogeneous boundary conditions. Therefore from (3.3) (u m-uo,T*v n-v 0)J = (v n-v 0,T U ] i i-uo)-and reduces equation (3.69) to I ^ - A 0 , A ° - A 0 ] -'2(v n-vo,Tu m-u 0); - ( V ^ V Q ~. (3.70) But from (3.42) (v n-v 0,T(u m-u 0))~ = (v n-v 0,v n-v 0)~. Substituting this into (3.70) yields [ A £ - A O . A £ - A 0 ] = (v n-v 0,v n-v 0);. Therefore, when the conditions in definition 3.4.3 of the energy product are satisfied, the error in the energy product is given by the mean square error in the stress v; [ A ^ - A 0 , A ^ - A 0 ] = ||vn-v0||2 (3.71) and the error in the energy norm as |A^-A0| = Fn-v0||. (3.72) Perhaps i t should be noted that i f the second of the matrix equations(3.35) i s satisfied exactly, the following are obtained for the energy product and energy norm: [ A ^ - A 0 , A ^ - A 0 ] = (A(u m-u 0),u m-u 0) = | | T ^ - T U Q || 2, and |VU°I = / ( A ( V U o ) ' V U o ) = l T V T u ° i ' The last two equations can be recognized as the errors in the energy product and the energy norm, respectively, as one would obtain in the Ritz method where A is a positive definite operator as in equation (2.11). 41 3.8 Extension It can occur that for some functions f e w there does not exist a func-tion A in the restricted f i e l d of definition that w i l l satisfy equations (3.35) Further, i t is important that an approximation be so chosen that only the kinematic boundary conditions are required to be satisfied. Such situations arise in the majority of problems in continuum mechanics. The governing differential equations of the type in (3.35) are derived under the assumption that the load f is continuous. This then requires that the f i e l d of defini-tion of the operator A be the totality of functions u and v defined over ™~ th R=R+dR that possess continuous m derivatives and satisfy a l l the boundary conditions of the problem on u and v. Here m is the order of derivatives in T. W Thus i f f is continuous then there exists a solution in D. but i f f A W is discontinuous, no solution can be found in D^ . This d i f f i c u l t y can be w overcome by considering limits of functions that l i e in D . Then i t is possible A to formulate the functional F(A) in such a manner that a generalized solution of equations (3.35) is obtained. Just as a discontinuous load may be con-sidered as a limit of a sequence of continuous loads, so functions with th discontinuous m derivatives are introduced that are the limits of sequences th of functions with continuous m derivatives. Thus i t can be said that amongst the new set of functions lies the solution (or generalized solution i f i t is not in D ) of equations (3.35) for any f e H . These ideas are now developed. A A Using equation (3.5), the f i r s t term on the right hand side in the energy product of equation (3.41) is replaced by (T*V!,u2)J= ( v i , T u 2 ) ~ - ( V l,Bu 2)~ 9 i ?. The boundary term is then deleted and the modified energy product in the symmetric form, denoted by [A 1,A 2] , is given by *"* A 42 [ A l , A 2 ] A = (v i,Tu 2)~+ ( T u l 5 v 2 ) ~ - ( v 1 > V 2 ) ~ . (3.73) T T w In which A1=<u1 v p and A2=<u2 v2> e D . The conditions of d e f i n i t i o n 3.4.3 A are s t i l l required to be s a t i s f i e d so that the positive definiteness holds. The energy product i n (3.73) also s a t i s f i e s the properties of an inner product ( d e f i n i t i o n 2.1.1). The energy norm i s now defined as |A|A = / [ ^ A ] " A (3.74) and also s a t i s f i e s the axioms of a norm presented i n theorem 3.4.2. The space w may be incomplete with respect to the energy norm, i.e. not a l l Cauchy ~ W W W sequences i n D converge to a function i n D . If this i s so, D i s completed A A A by defining A to be a member of the space i f |An - Al ->• 0, as m •> °° and n -> °°. (3.75) 1 -m -n - - T W -Where A =<u v > e D. since u and v are li n e a r combinations of the -m m n A m n sequences {<j> } and {\b }, respectively, which are i n Df. The com-, m n A W pleted cross product space so obtained i s denoted by H where the A subscript emphasizes the dependence on the operator A. The energy product i n (3.73) i s only defined for D^ 7 but may now be extended for a l l functions i n llW; A [ A 1 , A 2 ] A = lim U v l n , T u 2 m ) ~ + ( T u l m , v 2 n ) ~ - (v v )-} (3.76) m-**> V V V n-KX> where Aj=<u v, >T and A2=<uri v >T e D^ . Thus the energy product and im in 2m 2n A. r W the energy norm have meaning for a function A e H . A The f i e l d of d e f i n i t i o n of the functional F(A) i n equation (3.46) w w can now be extended from D. to H, and theorem 3.5.1 becomes: A A 43 Theorem 3.8.1 For a symmetric and indefinite operator A as defined in equation (3.35), of a l l the values which are given to the functional F(A) = [A,A] A - 2(gT,A) (3.77) by a l l possible functions A in H^ which also satisfy the conditions in definition 3.4.3 of the energy product, the minimum occurs only for the solution of equations (3.35). Proof: According to theorem 3.5.1 i f equation (3.35) has a solution A in w w D , this solution uniquely determines the minimum value of F(A) in D . It w w i l l be shown that the minimum value of F(A) in the wider class H is not altered and that the functions u and v only give this value. w — w Let d denote the minimum value of F(A) in D, and d that in.H.. A A W W -Then as H, includes D. , d d. A A - T W Assume d < d.. Then there exists a function A=<u v> e H, such that A F(A) < d, i.e. F(A) = [A,A] - 2(pT,A) = |A|2 - 2(f,u) < d. - W n T W But as A e H, then there exist sequences A =<u v > e D, such that - A -m m n A |An-A|->-0. Tnerefore from theorem 3.5.3 II u -u||-K). Thus |An|->-JAJ and (f,u )-»• '—m —' m " l-m' '-' m (f,u). Therefore for sufficiently large m and n, F(A) and F(A^) differ by an arbitrary small amount and i t follows that F(A^) < d. This i s impossible n W -as A e D.. Hence the contradiction shows that d=d. -m A To show that the minimum value of the functional is unique assume -that A e H also minimizes the functional. From the proof of theorem 3.1.1, i f the boundary conditions are homogeneous, equation (3.18) reduces to [(T*v-f ,n)~ + (Tu-v,w)~] = 0 (3.78) 4 4 w where u,n and v,w e D , and n and w are arbitrary variations of u and v, A respectively. Equation (3.78) can also be written as [A,A] = ( £ T ,A ) (3.79) where A=<n w>^ . If A is set equal to A, i.e. u=n, v=w, the equation (3.79) becomes [A,A] = (p T,A). (3.80) Equation (3.79) is also valid for any function in H^ . In fact i f A e i f A A then there exists A n e such that lA-A^-K), | u-u II+0 and [A ,A n]= ( p T , A n ) -m A l - -ml " m - -m - -m Proceeding to the limit gives equation (3.79) which is also valid for arbitrary functions in H^ . Thus A [A,A] = (p T , A ) . (3.81) By repeating the proof of theorem 3.1.1 with F(A) expressed in (3.77) the relation [A,A] = (p T ,A) is obtained where A=<n w>T is now an W " -arbitrary function in H . Putting A=A and A separately yields [A,A] = (p T ,A) (3.82) and [A,A] = (p T , A ) . (3.83) Subtracting (3.82) from (3.81) and (3.83) from (3.80) gives [A-A,A] = 0; [A-A,A] = 0. Finally subtracting one from the other of the last two equations yields [A-A,A-A] = 0. This implies that A=A. Therefore the minimum value of the functional and W W the function from H, which gives this minimum are the same as in D, and from A A theorem 3.5.1 such a function is also the solution of the equations (3.35). If the minimum of the functional F(A) in (3.77) is given by a w function A that is not in D , then this function is known as a generalized A. 45 s o l u t i o n o f e q u a t i o n ( 3 . 3 5 ) . To i l l u s t r a t e t h e s e c o n c e p t s , c o n s i d e r t h e e x a m p l e o f u n i a x i a l t e n s i o n o f a b a r w i t h u n i f o r m c r o s s - s e c t i o n . T h e d i f f e r e n t i a l e q u a t i o n s a r e e x p r e s s e d i n t h e m a t r i x f o r m a s f 00] o 0 - D 1 EA } • ( 3 . 8 4 ) a n d t h e b o u n d a r y c o n d i t i o n s a r e u ( 0 ) = 0 , P(L) = 0 ( 3 . 8 5 ) w h e r e n = ^ J P i s t h e a x i a l l o a d ; u , t h e e l o n g a t i o n ; f ( x ) , t h e d i s t r i b u t e d a x i a l l o a d a n d E A , t h e a x i a l s t i f f n e s s . F o r t h i s e x a m p l e , t h e o p e r a t o r s A , T a n d T * a n d A , g a r e g i v e n b y A = 0 - D 1 D E A , T = -T * =D ; A=<u P > T a n d p = < f ( x ) 0 > T . Now t h e e n e r g y p r o d u c t ' i n ( 3 . 4 1 ) i s g i v e n b y PiP-[ A i , A 2 ] = / { - u 2 D P i + P 2 D u ! gX~ } d x -I n t e g r a t i n g b y p a r t s t h e f i r s t t e r m o n t h e r i g h t h a n d s i d e y i e l d s , L rL P 1 P 2 [ A l f A 2 ] = - u 2 P i + / { P j D u 2 + P 2 D U l - -—-}dx, - - 0 0 J i A a n d d e l e t i n g t h e b o u n d a r y t e r m g i v e s t h e m o d i f i e d , s y m m e t r i c e n e r g y p r o d u c t o f ( 3 . 7 3 ) ; ( 3 . 8 6 ) [ A i , A 2 ] A = | L ( P i D u 2 + P 2 D U l P n P l r 2 - } dx . A J 0 " i — ^ - - i — i E A T h e f u n c t i o n a l F ( A ) i n ( 3 . 7 7 ) , f o r t h e e x a m p l e c o n s i d e r e d , b e c o m e s ( 3 . 8 7 ) F ( A ) = / { 2 P D u - | r } d x - 2 / f ( x ) u d x . 0 E A 0 ( 3 . 8 8 ) 46 If the solution A for equation (3.84) with f(x) continuous and the boundary W A' w conditions (3.85) was sought in D , then both boundary conditions on u W - - - T and P have to be satisfied by A e D,. However, the functions A=<u P> in - A H^ are defined such that A JA-An| = / ( L 2(P-P )D(u-u ) - (P-P ) 2dx 0; as m and n •> - (3.89) I - -m1 ' o n m r-A n n T W in which A =<u P > e D. and therefore u and P have continuous f i r s t -m m n A m n derivatives and satisfy a l l the boundary conditions of the problem. Such a definition means that functions u need only have generalized or piecewise f i r s t derivatives which in this case implies that u i t s e l f be continuous and - - - - T W P, piecewise continuous for A=<u P> e H . While the boundary condition to A be satisfied is that on u only, e.g. kinematic boundary condition u(0)=0, whereas P on the boundary turns out to be a natural boundary condition at the extremum of F(A) and P need not satisfy P(L)=0. In lieu of theorem 3.8.1, theorems 3.5.2 and 3.5.3 and the complete-ness definition of section 3.6 s t i l l hold. The only changes needed are the modified definition of the energy norm as in equation (3.74) and the extension of the f i e l d of definition of the functional from D*f to RW. A A In concluding this chapter, i t should be pointed out that in the development of the theory and the consequent proofs of the theorems, the problem considered involved s i n g l e dependent variables, u for displacement and v for stress. However, the theory is not limited to one dimensional problems and the extension to problems involving multi-dependent variables in displacements and stresses i s only a simple matter. 47 CHAPTER 4 FORMULATION AND CONVERGENCE OF THE MIXED FINITE ELEMENT METHOD It is well established that the formulation of a f i n i t e element method is a purely topological construction of a function sought as an approximate solution of a well-posed problem in mathematical physics, Oden [22,23]. The approximate solution i s expressed in terms of known coordinate functions and unknown parameters which can be determined by applying various different techniques, i.e. Ritz method, Galerkin method, least square method, etc. The technique employed here is the mixed method discussed in Chapter 3. A complete formulation of the mixed fi n i t e element method, i t s convergence and completeness c r i t e r i a are presented in this chapter. The strain energy convergence of the mixed f i n i t e element method for problems with stress singularities i s also established and a method for determining the stress intensity factor K^. is laid out. 4.1 Generation of a Finite Element Approximation In the f i n i t e element method, the f i r s t step is to replace the domain of the problem ft=ft+S by ft* such that ft* may be exactly subdivided into a number E of non-overlapping subdomains of simple geometrical forms called elements. Here, ft is the closed domain, ft, the open domain and S, the boundary. e The domain of a typical element w i l l be denoted by ft and adjacent elements are specified to have a common boundary. Thus ftm ft ftn = 0, m ^ n; m,n=l,2, . . . E where 0 is the empty set, and - F - p ft* = (j ft . e=l Further lim ft* = ft. E-*» 48 The elements are chosen, i f possible, such that ft* coincides with ft, but i f not, in such a manner that for a f i n i t e number E the error in the last equation i s acceptable. If so, for notation purposes ft and ft w i l l be used to represent ft* and ft*. The second step in the method involves the assumption of approxi-mate solutions for u and v in each of the elements and can be expressed in the form m u = E u, d>, k-1 k k n -e _ e.e v = E v ty k=l k k e=l,2, . . . E (4.1) wh ere <(>!" and tyT_ are the coordinate functions defined only in ft^" and u", v~ —e —e are the values u and v , or one of i t s derivatives, respectively, at certain nodal points generally situated on the boundary of the element. For example, e —e i f the u , k=l,2, . . . m correspond to the values of u at the nodes with i . m , coordinates x. then l (^(x™) = 1, i f k = m 1 k=l,2, . . . M. (4.2) =0, i f k $ m e e A similar relation also holds for ty^. Such a definition ensures that <J>^ e —e and ty^ are linearly independent throughout ft . A polynomial, which always contains linearly independent terms, with an appropriate number of coefficients (i.e., equal to the number of degrees of freedom in u or v) can always be put into the form (4.1) and leads to linearly independent (j^ and ty^. It i s only necessary to assume coordinate functions defined over individual elements to obtain a solution in the fi n i t e element procedure. However, i t can be rigorously shown that approximate solutions of the form (4.1) in fact do lead to approximations of the form (3.32) in terms of global degrees of freedom. It i s convenient to introduce two functions 49 e e — $k and defined over the entire domain fi such that —e for e fi = 0 otherwise; (4.3) and Y,!(x,) = ^ ( x 4 ) for x e fie (4.4) = 0 otherwise; where the x^ fe represent a point in the domain. Then the assumed approxi-mations for u and v throughout the whole domain fi can be written as E m e=l k=l k k (4.5) E n v = E E v ^ . e=l k=l fc k There has not been introduced any type of continuity in the generation of approximate solutions u and v in (4.5). Clearly u has (mxE) unknowns in u^ and v has (nxE) unknowns in v^. On interelement boundaries where nodes of adjacent elements coincide, i t is customary to specify that these nodal values should be the same. Assume that there are M and N in-dependent degrees of freedom for u and v over fi denoted by u_^ and v_^ , respectively. Then the element degrees of freedom are related to the global degrees of freedom by the relationship M uf = E F % , k . n ik l i=l (4.5a) N e e and v, = E G.. v. k . , jk j J=l J where e e F,, = 1 , i f nodal value u, coincides with u. ik ' k l = 0, otherwise. 50 Similarly e e ' G.. =1, i f nodal value v, coincides with v. jk k 3 = 0, otherwise. e e Therefore and G ^ are rectangular matrices in general and compare with the compatibility matrix one obtains in the usual displacement method relating the element degrees of freedom to global degrees of freedom. Then E m M j j ^ k u = £ £ £ u. F6, $f and v = E Z Z v. G6, f f e=l k=l j=l E n N Ze=l k=l 1-1 J J K k Defining E m cb. = z z F e $f J e=l k=l J k k E n and i|>. = Z Z G?. J e-l k=l J k k allows equation (4.5) to be written as M u = Z u. d>. i=l 1 1 (4.6) N v = Z v ,\b. j=l J J in which 4^ and xp are linearly independent throughout Q. Thus the fi n i t e element approximations of (4.6) have the same form as (3.32). The approximations for u and v within each element are chosen according to a certain continuity requirement and completeness so that the convergence of the fi n i t e element method in some sense can be guaranteed. This has been excluded in the discussion above and w i l l be discussed later. 51 4.2 Generation of an Element Matrix The derivation of an element matrix is described here for the problem of equation (3.30) and (3.31) over ft T*v - f = 0 (4.7) Tu - v = 0. (4.8) The element matrix can be derived either by using the variational principle F(A) of theorem 3.8.1 or by using the mixed Galerkin method of section 3.3 in which both residuals, that of the differential equations within the open e e domain ft and boundary S are included. Since the linear operators T and T* are self-adjoint, the matrices obtained by both methods w i l l be identical. Consider the variational formulation F(A) = [A,A] A - 2(pT,A) (4.9) T T where A=<u v> , p=<f 0> and [A,A] = J [2vTu - v 2]dft (4.10) (pT,A) = JQ fudft. (4.11) —Q. m e e —e ^ e e Let u and v be approximated by u = .E-<i.u. and v =.E,il>.v. within the element. 1=1 x i 1=1 l l Substitution of these into (4.10) yields, m n n n [A e,A el = / [2 E E <j;?T<j>eu!V:- E E i ^ V f v 6 ] d f t 6 - - A ' e . . . T r i l i i . N . T ± 3 l j ft i=l 3=1 J J 1=1 j=l m n n n or [A ,A ] = 2 E E u V f f ijjeT.<j>edfte- E E v V f f i ^ V f d f t 6 (4.12) X=l 3=1 J ft J 1=1 J=l J f t J and (4.11) gives (p T,A e) = I ueJ f<Oedfte. ( 4 - 1 3 ) " " i=I 1 V 1 e e e Now F(A ) in terms of u^ and v_. takes the form m n n n F(A e) = 21 E u e v e / tyeT^dnG- E E v V f / ^doT i=l j=l 1 2 ft6 J 1 i=l j=l 1 2 fte 1 2 - 2 E u ? J f<|>?dfte i=l 1 fte 1 52 3F(A e) 3F(Ae") and for a stationary point i t s derivatives — ~ c r - and „ ~ ' must equal zero. 9u? ov? i 1 This results in the following equations: Z f Tcb 7* dfi v - / e f ^ d f t 6 = 0; i=l,2, . . . m (4.15) m Z / 41 T i f d f i u - Z / ^?^ edfi ev? = 0; 1=1,2, j = i n J J j = i n J J n. (4.16) Equations (4.15) and (4.16) together yield (m+n) equations for (m+n) unknowns and can be put into a matrix form by defining a., and b.. as i l i j The matrix form is i . = / T O e d f t ; 1=1,2, . . u v i j j=l,2, . . m, b 6. = / yp%eda e; i=j=l,2. xj J^e r i r j . . n (4.17) p^ = / f ^ d f i e ; i=l,2, . . . n 1 • m. 0 (mxm) eT a (nxm) (mxm) (nxn) e u l e u 2 m e v l e v 2 e Pl e P 2 m 0 0 (4.18) 53 In general m need not be equal to n and matrix a can be rectangular. Furthermore, from equations (4.17), matrix b is symmetric and positive definite as long as the functions \b are square summable over fi . It can be observed that the matrix of coefficients of equation (4.18), i.e. 1 0 (mxm) eT a (nxm) a (mxn) nxn (4.19) (m+n) is symmetric but indefinite. The latter property follows from the non-extremum character of the mixed operators. 4.3 Assemblage of Element Matrices In the displacement approach, the stiffness matrix relates the nodal forces to nodal displacements. Hence during assemblage of element stiffness matrices, the addition of columns corresponds to equating the nodal displacements of adjacent elements to achieve certain compatibility 54 and add i t i o n of rows corresponds to equilibrium of the nodal forces. In the present case, the matrix of (4.19) rel a t e s both stresses and displacements (v and u, respectively) to the generalized applied loads. Although, the mechanics of the process of adding columns and rows i s s t i l l the same, the additions of both corresponding rows and columns of the element matrices are due to equating the nodal variables i n stresses and displacements with respect to c e r t a i n continuity requirements. This, i n f a c t , transforms the mixed problem of equation (4.14) from l o c a l (element) degrees of freedom i n u and v to global degrees of freedom. This i s analogous to the transformation £ £ of approximate solutions f o r each element of (4.1) i n u^ and to approxi-mate solutions i n global degrees of freedom u^ . and v. of (4.6). S i m i l a r l y f o r the mixed Galerkin method, i n which the r e s i d u a l of equation (4.7), a f t e r s u b s t i t u t i o n of an approximate s o l u t i o n i s made orthogonal to cf>? and that of equation (4.8) i s made orthogonal to \p^, the addition of corres-ponding rows and columns extends the orthogonalizing process from the element domain ft to the o v e r a l l domain ft. The submatrices a and b can be assembled separately or the complete matrix (4.10) can be assembled as a whole. The l a t t e r i s used i n t h i s thesis f o r analyzing problems. This requires interchanging of rows and columns of 6 e the matrix so that the degrees of freedom i n matrix equation (4.18) u i , u 2 , e e e e e e e e e e e e ' ' ' Um' v l ' v 2 > * " - v n a P P e a r a s u l > v l > u 2 > v 2 > • • • u m ' v m ' v m + i ' ' ' • v n f o r n m. This i s demonstrated i n Appendix A. By sui t a b l y arranging the glob a l degrees of freedom u^ and v.. the master matrix can be obtained i n an optimum banded form. I t was mentioned i n the previous section that the matrix of (4.19) i s i n d e f i n i t e . Therefore the method of Gaussian eli m i n a t i o n with p a r t i a l p i v o t i n g i s used to solve the glo b a l equations. 55 4.4. Convergence Criteria In Chapter 3 and especially section 3.8, the requirements that would ensure energy convergence in mixed methods for a continuum were presented. The same results are now used to provide sufficient convergence c r i t e r i a for a f i n i t e element approximation. Theorem 3.8.1 implies that the convergence of the mixed approxi-mation in the energy norm of equation (3.74) is ensured when applied to the equations (4.7) and (4.8) i f the coordinate functions <j> e U. and ifi e Vk m A n A W used for approximating u and v, respectively, are complete in H . This is explained for a f i n i t e element approximation by considering the following example (similar to (3.84)) representing the one dimensional linear e l a s t i c i t y problem with unit stiffness EA=1; u,v are displacement and stress, respec-tively, and f(x) is the applied load in the x-direction. -d l 0 dx dx -1 (4.20) Therefore dx' dx" (4.21) The homogeneous boundary conditions are u(0)=u(L)=0. It is convenient to w ~ construct a cross product space H =U *V\ in which A A A | A | A = / [ A , A ] A -J /J (2vf-v 2)dx (4.22) Then the functions A=<u v>T that are in rlf must satisfy the condition d(A,An) = | A-An| -v 0, as m -+ °°, n ->• °°; A n E D^. - -m' I- -mlA ' -m A (4.23) The space H A therefore contains the exact solution according to theorem 3.8.1. Thus the convergence would be ensured i f the coordinate functions <j> and il; m n W are complete in U^. 56 From section 4.1, a f i n i t e element approximation for u and v is obtained by dividing the interval (0,L) into E elements of lengths 1B and within each element, assuming solutions of the form m e=l,2, . . . E. (4.24) -e e.e k=l fc R The equations obtained within each element (section 4.2) are n e , e E / <j>* i|/Jdxv® = J1 fcb.dx; i=l,2, . . . m (4.25) j=l 0 1 1 J 0 1 m e , n c E / ipe<pe dxu e - E f i j j e^ edxv e = 0; i=l,2, . . . n (4.26) j=l 0 1 3 J j=l 0 1 3 3 where the prime denotes a derivative with respect to x. The equations (4.25) and (4.26) can be assembled to yield the following global equations: N . , , E / <fr.iji.dxv. * / f<)>_,dx; i=l,2, . . . M (4.27) j = 1 0 1 J J 0 x M L N L E / xp 4 dxu. - E / i^.^.dxv. = 0; i=l,2, . . . N. (4.28) j_2 0 1 J J j=1 0 1 3 3 Here <j>^ and ^ are the same as in equation (4.6) and M,N are the numbers of global degrees of freedom in u and v, respectively. Thus the global approxi-mate solutions for u and v M u = E <j> u i=l (4.29) N v = E \b .v. i-1 1 1 W converge in the energy norm i f the <b and ip are complete in H . The condi-tions that must be satisfied by the assumed solution within each element in order to ensure that the coordinate functions defining the global solution 57 W be complete in must therefore be established. For the example considered, the equations (4.22) and (4.23) W suggest that functions <J> in H must have generalized f i r s t derivatives i A w and therefore are continuous over (0,L). However, xp. in H over (0,L) 3 _ may only be piecewise continuous. Further, these functions must also w satisfy the forced boundary conditions as a l l functions in H are required to do so, and since u on the boundary takes on forced or kinematic boundary conditions, therefore <J>^ are to satisfy these conditions. In order to establish a criterion which leads to satisfying the completeness requirements of section 3.6, their definition in terms of fi n i t e element approximations is necessary. This is done next. Completeness for Finite Element Approximations The finite element approximations for u and v as linear combinations of sequences of functions d> and \p , m and n = 1,2, . . . , in HW are complete m n A * i f for a A with a f i n i t e energy norm and any e^O, £2>0> there exists a subdivision of the domain of the problem corresponding to M and N degrees of freedom in u and v, respectively, such that N (i) ||Tu^ - I V.4/J < e 2 i=l ( ^ ) | A " $ < E l M N where A=<~u v>T, A^=<uM v N> T and ^ . E^.cb. , v^.^v.^.. Again, the completeness requirement (i) implies that the stress approximation, N v =.E,v.il>., should contain a l l the stress modes that correspond to strain N J=l J J modes present in the strain Tu^ for f i n i t e degrees of freedom N in stress v. If the operator T involves derivatives of order m^ , then the energy 58 T T product in equation (3.73) for A =<u. v.> and A.=<u. v > - i i i -J J J [A.,A.]. = (v. ,Tu. ) + (v. ,Tu. ) - (v. ,v. ) (4.30) - i ' - j J A v jn' im'v i n ' jm'v v in' jn V \t.->w is symmetric and involves derivatives of order m^, in u and zero in v. To establish certain completeness criterion for f i n i t e element approximations for u and v within each element, interest lies in the least requirement that would allow the energy norm in (4.22) to be evaluated, i.e. requirement ( i i ) , and also enable the mean convergence of approximation for v to Tu, i.e. requirement (i) of completeness. Therefore continuity of (m^-1) derivatives is required for u. However v can be piecewise continuous. Completeness criterion for the fi n i t e element approximations within the element i s presented in the following section for linear e l a s t i c i t y plane stress. The explanation of the notation used for certain quantities can be found in the beginning of Chapter 5. 4.5 Completeness Criterion A general criterion for completeness is stated and j u s t i f i e d here and is analogous to the one presented by Oliveira [24,25]. Let H^ be the set of compatible and equilibrated elastic fields where u^ have continuous and bounded second order derivatives and v. have continuous and bounded f i r s t l order derivatives within each element. It is to be demonstrated that the ne e e T e completeness for a f i n i t e element approximation A =<u T > where u is the ^ r -m -m -n -m vector of displacement components and T , the stress components, within the T W element with respect to any A=<u x> e Hg w i l l be obtained i f : — — (a) the general analytical expression for u^ and T within each element are polynomials with the number of arbitrary parameters equal to the number of unit modes corresponding to the element, 59 (b) the terms of degree higher than the f i r s t in u^ and the constants in e x,. can vanish regardless of the values taken by the constant terms i j e and the coefficients which affect the linear terms in u_^ and the constant terms in x. . . i j (c) the constant terms and the coefficients which affect the linear terms e e in u^ and the constant terms i n are completely arbitrary. e e (d) x,. should at least possess a l l the modes that are present in T,,u,. The conditions (a) to (c) result in displacements u_^ or their f i r s t derivatives and stresses T to take up any arbitrary value throughout e e the element. The condition (d) assures convergence of x.. to T..u. in the i j i j 3 mean square sense, i.e. completeness requirement ( i ) . w The displacements u. belonging to H 0 can be represented inside fi by the following Taylor's expansion u = u (P) + u (P)(x -x^) + 1/2! u i k ( P ) ( x - - x ^ < V X k } ( 4 ' 3 1 ) where P and P are points in 0, and P depends upon the coordinates of the point where the u. are to be determined. Similarly x.. can be represented by the following Taylor's expansion T i j • T u ( p ) + Ti j >k c 5> (V< )- ( 4- 3 2 ) Now consider the displacement f i e l d with components u*\ = u (P) + u (P) (x - x ) ( 4 . 3 3 ) t i i i >J j J within fi and c a l l i t a tangent f i e l d to u^ at P. Similarly consider the stress f i e l d x6.. = x ..(P) ( 4 . 3 4 ) c i j c i j within fi and c a l l i t a constant f i e l d at P. From hypotheses, a l l second derivatives of u. and f i r s t derivatives of x,. are supposed to be bounded i i j inside Q . Then equations ( 4 . 3 1 ) and ( A . 3 3 ) yield 60 • U i " U t i l K 2 T L l ( I e ) 2 ( 4 ' 3 5 ) e where Lj is the upper bound for a l l the second derivatives of u^, I is the maximum diameter of element e and d is the total number of second order derivatives. Similarly from equations (4.32) and (4.34) K j - T c l j l < f L" (V (4-36> e where L 2 is the upper bound for a l l f f i r s t derivatives of x . By considering similar expansions for the f i r s t derivatives of u_^ , i t is possible to derive the inequality l ui , j " u t i , j l < ( 4 - 3 7 ) The distance between A and A over the entire domain Q, can be evaluated from the energy norm of equation (3.74) E . _ d(A, A) = / Z [A-A.A-A] - - n - - - - A e=l e e T where E is the total number of elements and A=<u x > . By using equation (5.13) d(A,A) = / Z / [2 ( x e - x e) TT(u e-u e) - ( x e - x e ) T C ( x e - x e ) T ] d f i e - - / , ' „e - -c - - -t - -c - - -c y e=l n or d(A,A) < / z [2|f ( x e - x e) TT(u e-uf)dQ e| + 1/ ( x e - x e ) T C ( x e - x e ) d f t e I ] . - - -\ / i ~e - -c - - - t ' 1 J - -c - - -c ' Y e = l tt tt As the operator T involves f i r s t derivatives and C is a compliance matrix which involves no derivatives, by virtue of equations (4.36) and (4.37) E ? C 7 ^2 . ,.?/, ^,2X0*1 d(A,A) < / Z [{Lf(l ) z + L^Q ) 2}ft c] " " ' e=l for _ e sufficiently small. L 3 and are positive numbers. The inequality above can be further simplified; < 61 d ( A , A ) < / ( L 3 + L " ) I 2 E fi ( 4 . 3 8 ) e=l xi n e e where L^yL^ and 1^ denote the maximum values of L3 ,Li4 and I in the whole set of elements. Now i f a=(Ln+L,"), then d(A,A) < / a l z f i . (4.39) - - n Therefore the energy distance between A and A tends to zero when the size of elements decreases indefinitely or the number of elements E goes to in f i n i t y . Consider now a type of element generating a sequence of families W N T of fields, H^, whose completeness Is to be investigated. Let A^=<u^ j_N> be one of these fields whose generalized parameters or nodal degrees of W n 6 freedom take the same values as that of A from Hn., and A corresponds to UA -m N e A., within fi . If the general criterion i s satisfied then A is one of the -M 0 -family of fields which can occur within the f i n i t e element. Let such a e e f i e l d correspond to values u , and x ... of the parameters, i.e. L l K , C X J K. e m ke e u t i = E A u t i k ; in fi , (4.40) k=l n x-6. . = E i\> ex e. ; in fie, (4.41) cij k = 1 i j cijk* or A 6 = <uf f°> T ; in fi6. - - t -c e On the other hand for A , -m e m ke e e u . = E A. u ., ; in fi , (4.42) mi , T r i mik k=l ' e n ke e e x ., = E U). . x ... ; in fi , (4.43) nij k = 1 r i j mjk .ne e e T . ^e or A = <u x > ; in fi . -m -m -n 62 From equations (4.40) and (4.42) m . I ~ e e 1 I v ,ke. e e . k=l which then leads to a difference of derivatives as m . .e . e i i „ ,ke , e u . . - u e . .1 = I E <f> e. ( u . . - u ., ) I ; i n if. (4.44) t i , j m i j 1 ' k = 1 T i , j ^ t x k mik y 1 Similarly | x 6 . - x e . | = | E 4v(x e. ., - xe... )| ; in f f . (4.45) 1 c i j n i j 1 ' k = 1 r i j cijk nijk Ice Ice Ice Since the coordinate functions 4>. ,<J>. . and \p. . can be expressed in dimension-x-less form, i.e. as a function of • j - L , these functions remain bounded as the e size of the element decreases. Assume that the moduli of the magnitudes of Ice _ce Ice e a l l functions <j>. ,<j>. . and \b. . remain below Lc. a positive number. Then rx ' Tx,j yxj 5 ' * m |if. - u e. | < Lf E I uf.. - u6., I (4.46) 1 tx mx1 , , txk mxk1 k=l T e m |<. . - u e. .| < ^ E |u^.: - u E . J (4.47) 1 tx,j mx,j' I , 1 txk mxk1 e ie x n l x 6 . . - xe..| < L 5 E |xe... - x6.., I (4.48) 1 cxj nxj 1 3 k = 1 1 cxjk nxjk' within 0, . e e e Since A and A take the same values at nodes, i.e. u , and x.., : - ik i j k ' equations (4.35) and (4.36) permit one to write l U t i k - U m i k l < i r ^ V 2 W . 4 9 ) K i j k - T n i j k l < fLtCV' (4-50) Introducing (4.49) and (4.50) into (4.47) and (4.48) yields 63 I T 6 , . - x e..| < f L f L ? n ( l e ) = L?2 e. (4.52) Equations (4.51) and (4.52) s t i l l hold even i f u . and x .. are discontinuous, mi,j n i j as long as they are bounded. From s i m i l a r i t y between (4.36), (4.37) and (4.52), (4.51), r e s p e c t i v e l y , the following i n e q u a l i t y which i s analogous to (4.39) can be written d i r e c t l y d(A,A^) < /3l£n. (4.53) Combining equations (4.39) and (4.53) and using the t r i a n g l e i n -equality of the energy norm d C ^ , A ) < d(AjJ,A) + d(A,A) < + / B ) or |AJJ - A | < ( v £ " + / B ) ^ T ( I ) . (4.54) N Equation (4.54) means that the energy distance between A^ and A tends to zero as the s i z e of the largest element decreases i n d e f i n i t e l y . Note, condition (d) of completeness c r i t e r i o n i s s a t i s f i e d f o r the case where u. . and x.. are piecewise continuous and taken as constants within each e e element, i . e . the best u. can represent i s constant s t r a i n and T , . constant i i j s t r e s s . Although, from the completeness c r i t e r i o n , continuity of x i s not required; improved accuracy i n the energy convergence i s observed i n the res u l t s (Chapter 7) by making x continuous across element boundaries. Also, i f the stresses x ^ are not continuous, and completeness requirements s a t i s f i e d , the mixed f i n i t e element method y i e l d s the same r e s u l t s as the displacement method would using i d e n t i c a l i n t e r p o l a t i o n s f o r the displacements, 64 Theorem 4.5.1 The mixed f i n i t e element method, with approximations for the displace and stresses satisfying completeness requirements, yields exactly the same results as the displacement method would using the same approximations for displacements i f the stresses were not continuous across the interelement boundaries• Proof: Suppose with the f i n i t e element solution of equations (4.5) for u and v the functional F(A ) of (4.9) has been arrived at. Using (4.17) and f=0,F(A ) can be expressed in the matrix form T T T T / , rr\ ,.e. e e e e e e e e T. (4.55) F ( A ) = u a y + y a u - y b y Assume u and v to be continuous and let u 6 = ru (4.56) y e = sy (4.57) e e where r and s relate the total element degrees of freedom Q and v to the global degrees of freedom u and v so as to restore interelement continuity N for u and v, respectively, equations ( 4 . 5 A ) . Therefore F(A^) for the entire domain can be expressed as N T T T T T T T F(A^) = u r asy + y s a r u - v s bsv (4.58) e e where a and b matrices are formed by collecting the a and b matrices for each element in uncoupled form, respectively, and substituting (4.56) and (4.57) to couple them. M and N are the global degrees of freedom in u and v, respectively. If only the displacement u is continuous, then T T _e, .eT T ~eT, %e c n . M = H E a v 2 a I ^ - Y k ^ (4.59) where Ni~eN because at common nodes the stresses are not equated. 65 But within each element eT e , e e a u = b v and since b is an (nxn) symmetric positive definite matrix, i t can be inverted to give e , e _ i eT e v = b 1a u also ^ e = b ^ a V 2 = b _ 1 a T r u . (4.60) Substitution of (4.60) into (4.59) yields w „ N l N T T . i T T T , _ i T T T , _ i T F ( A * ) = u r ab 1a ru + u r ab 1a ru - u r ab 1bb xa ru. - M - - — - — - - — - — - - — — - — Since b _ 1b = I; therefore FCAJI 1) = u Tr Tab- 1a Tru. (4.61) - M - - — - -If i t can be shown that the elimination of v in the matrix equation (4.18), i . e . a 6b e 1ae'"'u = k 6 u 6 = p 6 (4.62) -e and that k is the same as the element stiffness matrix one would obtain from the displacement method using the same approximation for u, then the theorem would hold. In a system of fi n i t e degrees of freedom, the requirement (i) of completeness requires that the approximation for stress should contain a l l the modes present i n Tu , the strain f i e l d from u . Then no matter how much m m better the stress field v is chosen to be, for discontinuous stress fields n across the element boundaries, the Fourier type convergence of the strain f i e l d from v cannot be better than Tu . Further, from equation (3.68) n m ||u' - E v.ip.W2 = ||u* l | 2 " " v? = 0. (4.63) 1=1 1=1 du Here Tu =-r—=u' for the example considered in equation (4.20). m dx m 66 Then, for orthonormal \> , equation (4.28) yields v = (u' A ) ; 1=1,2, . . . n. (4.64) l m i Also matrix b is just an identity matrix I. Therefore e, e_ i eT e T eT e eT , n a b xa = a la = a a . (4.65) From (4.63) n if _*|| 2 = Z v 2. (4.66) 1=1 1 m But from (4.1), u' = E i f ' u . , therefore m . i 1 i i=l m m | | u j 2 = E E / C <f,e'dx u 6u e *" i=l j=l ^ 1 J n m m n. and from (4.64) E v^ = E E )(f | ? ' ) u u ; k = l , 2 ( . i=l 1 i=l j=l 1 k k 3 1 3 For the problem considered and e ,,et . x a ± j = ( * ± - . * j > • Therefore the matrix form of equation (4.66) is eT, e e eT e eT - _ _ = _ a a u e or e e eT k = a a . (4.67) The functions tfj.'sused, were assumed to be orthonormal for convenience, l However, no generality i s lost and for any linearly independent ty^s equation (4.67) can be written as i T k e = a e b e _ 1 a e . (4.68) Therefore equation (4.61) takes the form 67 F(A N l) = i/r^kru = i^Ku -m T T, (4.69) where k is the assembled uncoupled global stiffness matrix and K is the coupled form as one would obtain from the displacement approach. Therefore, the energy convergence cannot be any better in the mixed method when the stresses are discontinuous across the interelement boundaries. This then completes the proof. 4.6 Stress Singularities Stress singularity at the crack tip in plane ela s t i c i t y problems, i t ' s influence on strain energy convergence and estimation of it s crack intensity factor of mode type I (the opening mode, Figure 1) are considered in this section. The displacements and stresses near the tip of a sharp crack (mode type I) are given by v K /2r 8G/TT~ (2K-1)COS| - Cos|^-(2K+l)Sin| - Sin|^ + 0(r); (4.70) -T X X I T yy v2iTr X . xy I2'Sin2~> rj.Siny-, , n 3 6 .osy + X X yy xy + 0(/r~) (4.71) Here, K_ is the stress intensity factor; K takes the value (3-4v) for plane strain and (3-v)/(l+v) for plane stress state; v is Poisson's ratio and G, the shear modulus. As for the constant stresses S , S and S in equation xx yy xy (4.71), the stress free condition on the crack surface near the tip leads to S =S =0, whereas the component S in the direction of the crack line remains yy xy xx unknown. 68 The strain energy convergence of the mixed f i n i t e element method for stress singular plane elasticity is established f i r s t , and then a procedure for determining the stress intensity factor K^. is presented which is somewhat similar to the one presented by Parks [27] for displacement method. 4.6.A Strain Energy Convergence for Problems with Stress Singularity In plane elasticity with stress singularities, the convergence rate for the displacement f i n i t e element method is often controlled by the nature of the solution near the singularities. Unless the singularities are properly handled, the regular so-called high accuracy element w i l l not be able to improve the convergence rate. Tong and Pian [38] showed that the convergence in strain energy for displacement and hybrid type elements is only of order 1 or 0 ( j ) where 1 is the maximum size of the elements J e e e near the crack tip. They established i t by arguing that in stress singular problems, the typical polynomial type approximations for the displacements lead to strain fields which exclude the singular stresses near the crack tip. Hence the missing terms govern the error in strain and consequently lead to slower convergence of the strain energy. A similar argument is followed for establishing the strain energy convergence for mixed methods applied to problems with stress singularities. In linear e l a s t i c i t y , twice the potential energy for zero body forces can be written as (from equation (3.77)) F(A) = [A,A] A- 2/s T ±u 1ds (4.72) T where [A»A]A= / v (2[Tu]Tx - i TCT)dV (4.73) represents twice the strain energy and C, T*, T, u and T are the extensions 69 to three dimensions from plane stress el a s t i c i t y , equations (5.10) and T^'s are the surface stresses specified on part of the boundary. Let An T =<u0 xn> be the exact solution of the problem over a domain V. For simplicity, assume that the solution is sufficiently smooth except that u 0 = r ag(x) (A.74a) x 0 = r a _ 1 h ( x ) (4.74b) near the point of the singularity R in V. In equations (4.74), r is the radial distance from R, g and h are smooth functions, x are the spatial co-ordinates, a is not an integer and Pl + 1 > a + j > 1 and p 2 + 2 > a + -| £ 1 (4.75) where n is the spatial dimension of the domain V, p 1 and p 2 are the degrees of polynomials used as approximations for u and T, (satisfying completeness) respectively. The f i n i t e element approximations can be put into the form (section 4.1) m u = Z $.(x)u. i=l (4.76) n T = E f.(X)T.. i=l 1 " x Let the solution obtained from the mixed f i n i t e element method be m _ u = Z <b (x)u 1=1 (4.77) n T = E 4>. (x)x. " i = l 1 - - 1 and the solution by forcing the nodal degrees of freedom to take the exact values of A 0 at the nodes, be 70 m ii = E $ (x)u. i=l 1 " _ 1 r = E m.(x)x . . i=l 1 1 (A.78) Using equation (3.49), inequality (3.56a) yields [A-A 0,A-A]< [A-A0,A-A0] (4.79) _ x T T where A = <u T> , A = <u x> and A 0 = <UQ T Q > . Then from theorem 3.5.3 / v(u i-u. o) 2dV < X![A-A0,A-A0] (4.80) / V ( T . J - T . j 0 ) 2 d V < A2[A-A0,A-A0] (4.81) where \± and A2 are positive constants and A} depends on the lowest non-zero vibration frequency. The inequalities (4.79) to (4.81) show that the mixed fi n i t e element solution is bounded by the rate at which |A — A Q| approaches zero. If the interpolations are chosen such that the equations (4.78) can exactly represent any polynomial of degree jpi for displacements and p 2 for stresses; further that way from the singularity, the Pi+1 derivatives of ug and p 2+l derivatives of tg are bounded then i t can be shown (section 4.5) I u. .-u? .I ^ A J r a 1 in V i 1 i,3 1,3' 1 1 A.h pl |u. .-u? .1 « — ^ - z in V-Vi; i,3 i,3 „Pl+l-a I T . , -T? . | « B. . r a 1 in V i 1 13 13 1 13 B. .hP 2 + 1 T . . - T ? . | < - i J — in V-Vi. 13 13 r P 2 + 2 - a (4.82) (4.83) 71 The positive constants and B are, respectively, the combinations of the upper bounds of a l l Pi+1 derivatives in the Taylor expansions of u~ and a l l P2+1 derivatives for . is the domain covered by the elements adjacent ,0 L, U L . L X V a L - X V V ^ O -L. L L L U C -L. CL y X U X. C A U < a i l O X U U O W J- w .0 i -to R and h is the maximum size of the elements. From equation (4.73), the energy product for A-AQ is given by [A-A 0,A-A 0] A = / v {2[T(u-u 0)] T(x-x 0) - (x-x 0) TC(x-x 0)}dV (4.84) or [A-A0,A-A0] < 2Jv\ [T(u-u 0)] T||(x-x 0)|dV + / y|(x-x 0) T|C|(x-x 0)|dV. (4.85) As the operator T involves f i r s t derivatives and C is a positive definite compliance matrix, by virtue of (4.82) and (4.83) hPl+P2+l h 2 ( P 2 + l ) [A-A0,A-A0] « / V i { c 1 r 2 ( o - 1 ) + c 2 r 2 ( a - 1 ) } d V + / ^ ^ P i + p ^ a + C\2(p2+2-a)}dV (4.86) The constants C j and c 2 are positive; c^ depends on A_^ and B and c 2 depends on B.. and the compliance matrix C. Assuming that the inequalities in (4.75) are satisfied, integration of the right hand side of (4.86) yields ~ O f _ 1 W u P l + P 2 + l u 2 ( P 2 + x ) [A-A0,A-A0] < k 1 ( c 1 + c 2 ) r Z ^ a ± } +k2c1-JL-— +k 3c 2 J) > . (4.87) - - U - - U 1 1 ^ M A X ^ J. Pl+p2+3-2a-n 0 ^ 2(p2+2-a)-n min min The constants k±,k2 and k 3 are also positive and depend on the geometry and the arrangement of the f i n i t e element mesh; r and r . are the maximum and max mxn the minimum radial distances from R to the boundaries of , respectively. From inequality (4.87), the contribution to error in the energy product from the elements immediately adjacent to the point of singularity is of order 72 i"2(a l ) + n and the contribution from the rest of the domain is of order either max h P l + P 2 + l h2(P2+D or — — , — — — r — , whichever is smaller. Pl+P2+3-2a-n 2(p2+2-a)-n min min Further, the main contribution from the latter part is from the elements that are close to the singularity. In the f i n i t e element analysis, r , r , and ° max min h are usually of the same order, i.e. r ^ r . ^ h. Thus inequality (4.87) can max mm be rewritten as [A-A0,A-A0] < UiO^+kz) + c 2 ( k 1 + k 3 ) } h 2 ( a _ 1 ) + n . (4.88) It should be noted that the constants c and c , through their dependence on A. and B.., also depend on the behaviour of u and near the point of x xj -singularity. Now from inequality (4.79) [A-A0,A-A0] < {ciCki+kz) + c 2 ( k 1 + k 3 ) } h 2 ( a ~ 1 ) + n . (4.89) Inequality (4.89) shows that the order of convergence of the energy product is controlled by the order of singularity, rather than by the order of the polynomials used for interpolation of u and x provided the inequalities (4.75) hold and completeness criterion (section 4.5) is satisfied. For plane e l a s t i c i t y , from equations (4.70) and (4.71), n=2 and a=l/2. Therefore 2(a-l)+n=l and for any Pi»l and P2>0, the inequalities (4.75) are satisfied. Inequality (4.89) now becomes [A-A0,A-A0] < {c^kj+kz) + c 2(k 1+k 3)}h. (4.90) 73 Therefore i t can be concluded that the strain energy w i l l converge at least as h and perhaps faster i f there is a cancellation of errors in the energy product of (4.84) due to singular terms which are lef t out when regular polynomials are used for approximating displacements and stresses. 4.6.B Determination of Stress Intensity Factor K^. The crack intensity factor K_ (mode type I crack, Figure 1) is related to the potential energy release rate G_, i.e. the rate of change in the potential energy due to crack extension. For plane stress and plane strain problems with unit thickness; this relationship i s given by the following equation, Rice [29]. where TT is the potential energy, a the crack length and G and K are the same as defined for equations (4.70) and (4.71). It has been concluded by many investigators, Watwood [39], Anderson, et al. [1], Parks [27], etc. that the most accurate f i n i t e element scheme for determining the stress intensity fact is by means of the evaluation of the energy release rate without requiring any special elements to model the stress singularities. This method i s now applied to determine K_ for plane strain linear e l a s t i c i t y using the mixed f i n i t e elements. For plane strain, the equation (4.91) reduces to °i - - H - - v (4-92> where E is the Young's modulus. The procedure used in the application of the mixed f i n i t e element 74 method is somewhat similar to the one presented by Parks [27]. For two dimensional planar configurations, an alternative computation of G_ can be based on the path independent J-integral, Rice [29]; J = /_ Wdy - T.|^ds. (4.93) Here W is the strain energy density (W=l/2 x..e.. for linear e l a s t i c i t y ) ; T i s an arbitrary contour enclosing the crack tip, Figure 2; T is the traction vector associated with n, the outward unit normal and u is the displacement vector, and ds is an element of arc-length along T. Near the crack tip, the parameters comprising the integrand of the J-integral as determined from a f i n i t e element solution are likely to be in greatest error. However, the path independence allows the* contour to be chosen farther away from the crack tip, hopefully resulting in improved accuracy. Rice [29] also showed that the J-integral in fact is equal to the potential energy release rate, euation (4.91). J - G i = < ¥ ^ - TJ-Suppose that the f i n i t e element analysis has been performed on a given planar linear elastic body of unit thickness containing a crack. In the discrete form, the potential energy in the mixed method can be expressed as T I T T i = T Au - r t B T - f u (4.95) M - — 2- — - -where the matrices A and B are separately assembled element sub-matrices e e a and b of equation (4.18), f is the generalized load vector, and u and T are the solution vectors. The equation (4.95) is analogous to the expression for twice the potential energy in a continuum, equation (5.12). The potential energy release rate can now be obtained by differentiating TT^ with respect 75 to the crack length a: 3TT T T M - 8 - TA u i j . TA * i 9 u . T3A 1 T3B df 1 7 " + 9a - ¥a^ ~ 2^ 3 ? ~ i a ~ ^ ' (4.96) But from the f i n i t e element analysis, A U = B T and Ax=f; therefore the f i r s t two terms on the right hand side are zero and (4.96) reduces to 3ir. _M = 1^ T T 3 3a "2s- - ^ 0 A " u T A -B T — _af1 3 ^ (4.97) Furthermore, i f the loading on the body is accomplished by surface tractions applied on the boundary other than the crack face, then the load vector f 3 f is independent of the infinitesimal crack advance, i.e. -r^O. Therefore, da from the equations (4.94) and (4.97) 3TT M 3a ( — ) k 2 = - ^ u T > ^ [ 7 ] , (4.98) where S is the master f i n i t e element matrix; "0 A I S = T A -B (4.99) 3S and — represents the change in the master f i n i t e element matrix per unit crack advance. Aa, the infinitesimal extension of the crack tip can be approximated by rigidly translating a l l nodes on and within a contour T 0 about the crack tip in the x-direction, Figure 3. A l l other nodes remain in their i n i t i a l positions. Thus the master f i n i t e element matrix, which depends on individual element geometries and elastic material properties, remains unchanged in the regions interior to T 0 and exterior t o T j , and the 3S only contributions to -r3" come from the band of elements between the contours 3a T 0 and T1. Thus 76 F° 3S° 3TTM T T 3_ru, .T .T * -x ru, — = <u T > — [-] = <Q t > E - — [ 7 ] (4.100) da - - da T - - . . da f 1=1 -where Q and 1? are the nodal variables for the nodes on TQ and Tj; E° is the number of elements between the contours TQ and and S° their element matrices. The change in the element matrices can be calculated directly as SS? 3S? 3x. = T^" (4.101) 3a 3x. 3a 3 where the nodal coordinates x^ are thought of as functions of the crack length a. The derivatives 3x_./3a are then unity or zero, depending on whether or not x^ . is the x-coordinate of a node located on TQ, respectively. Alternatively, 3S^/3a may be approximated by a simple forward f i n i t e difference scheme 3S? AS? "v —• - — -• = "T [S? -S? ]. (4.102) 3a Aa Aa -x -•• "a+Aa """a Here S° are the element matrices for the elements between the contours T 0 and Ti, calculated for the i n i t i a l crack length a, and S? when x coordinate a+Aa of each of the nodes lying on TQ have been incremented by Aa. The equations (4.98) and (4.100) suggest that to calculate the stress intensity factor R^., the master f i n i t e element matrix equation need only be solved once, i.e. for the i n i t i a l crack length a. After obtaining this solution, pre- and post-multiplying the differentiated element matrices of equation (4.102) with the solution vectors for the corresponding nodal variables and then summing these over a l l the elements between TQ and Tj yields the rate of change of potential energy in the discrete sense ATr„/Aa. Alternatively, the differentiated matrices can be assembled and then M pre- and post-multiplied by the solution vectors u and t for the nodal dis-placements and stresses, respectively, for the nodes on TQ and Ti , as in 77 equation (4.100) a+Aa -so ][„]. (4.103) Finally, the stress intensity factor Kj for plane strain state can be computed from the requirement that i t be internal to the body and enclose the crack tip. It can also be shrunk to a single node at the crack tip so that the sum-mation in (4.100) extends over the elements adjacent to the crack tip only A glance at Figure 3 and the procedure outlined above for determining the potential energy release rate indicate that It is an area-analogue of the path independent J-integral. (4.104) The contour TQ to be translated is thus far arbitrary except for 78 CHAPTER 5 APPLICATION OF BOUNDARY CONDITIONS So far only the homogeneous boundary conditions have been considered in the development of the theory of mixed methods. The treat-ment of homogeneous, mixed homogeneous and non-homogeneous boundary condi-tions , how these can be incorporated in the mixed f i n i t e element method, and equivalence to boundary residual concept are presented in this chapter. For i l l u s t r a t i o n purposes the plane stress linear e l a s t i c i t y problem with unit thickness i s considered. The governing differential equations are - T . . . = f. (5.1a) T.. = 2ue.. + Xe 6.. (5.lb) 13 13 kk 13 e = l/2(u, , + u. ) (5.1c) i j i»j 3,1 where x „ and e are symmetric second order tensors, and i=j=l,2. Equations (5.1a) are the equilibrium equations relating stress gradients to the body forces f^, (5.1b) are the constitutive equations where X and u are Lam£s constants and (5.1c) are the kinematic equations relating strains to displacement gradients. Assuming that the equations (5.1c) are satisfied identically equations (5.1) reduce to the following set of equations - T = f ; i=j=l,2 (5.2) ••-J >J -1-1 / 2 ( U i , 3 + U 3 , i ) " W k l " 0 5 i-J-K-J-1.2. (5.3) Where C^j^j is the fourth order compliance tensor. For E, the Young's modulus and v, the Poisson's ratio defined as 79 _ u(3A+u), _ J (X+y) ' V 2(X+u) the equations (5.2) and (5.3), for a plane stress problem, can be written in the following form: x - x = fx xx,x xy,y - x - x = fy xy,x yy,y (5.4) (5.5) u, - - (x - vx ) = 0 x E xx yy v, - - (-vx + x ) = 0 y E xx yy 2(l+v)_ u, + v, - — x = 0. y x E xy (5.6) (5.7) (5.8) The equations (5.4) to (5.8) can now be put into matrix form 0 0 0 0 5 3x 0 0 T-3y 3y 3x 3x u 3y 0 _3y_ _3x ^ - i 0 E E 0 0 -2(l+v) v XX yy xy which takes the equivalent form of (3.35) as "0 T*l IT -C u f X 0 0 0 0 (5.9) (5.10) or Where AA = p, u = <u v> ; or = u and u 2 X = <X X X > , xx yy xy v X* = - x = - i - o 3x 3y 0 3_ 3_ "3y ~3x (5.10a) (5.10b) (5.10c) (5.10d) 80 T T T A = <u x > and f = <f f > , x y and the compliance matrix C is a symmetric positive definite matrix 1 -v 0 ,C=f -v 1 0 0 0 2(l+v) The energy product of definition 3.4.3, for unit thickness, becomes (AA,A) = [A,A] = / [uTT*x + xTTu - xVrjdfi. (5.10e) (5.11) This, on integration by parts of the f i r s t term on the right hand side, yields [A,A] = - £ x u-nds - $ x u-sds + / [ 2 T T T U - T T C T ] d f i (5.12) where T and x are stresses normal and tangential to the boundary and n nn ns J -and s are unit outward normal and tangential vectors, respectively. If the boundary conditions are homogeneous, the boundary integrals i n (5.12) can be dropped resulting in the energy product of equation (3.73). fA,A] A = / [2x Tu - x Cxjdfi. (5.12a) The mixed variational principle of (3.77) for plane stress linear el a s t i c i t y can now be written as F(A) = J n [2xTTu-xTCx]dft - 2 / f Tudfi. (5.13) T W and A=<u x > E H.. - - - A 5.1 Homogeneous Boundary Conditions The homogeneous boundary conditions for plane stress can be expressed as on S m x . . n. =0 u. = 0 l (5.14) on S u 81 where ru's are the components of unit outward normal, S_ and are the portions of the boundary S where the stresses and the displacements are prescribed to be zero, respectively. The element matrices can be generated directly from (5.13) (Appendix A) and assembled by using the procedure discussed in chapter 4. Since the stresses are incorporated into the functional F(A) as natural boundary conditions only the kinematic boundary conditions on u are to be satisfied. This can be achieved by forcing the corresponding displacement nodal variables to be zero on the boundary S^. The process is identical to forcing the homogeneous kinematic boundary conditions in the displacement method. 5.2 Homogeneous Mixed Boundary Conditions The boundary conditions for the plane stress problem of equations (5.4) to (5.8) in this case are u. = 0 on S x u T . . n . = 0 on S_ (5.15) T . , n , + au. = 0 on S„ xj j x M where S^, S_ and S^ are the portions of the boundary S on which displacements, stresses and mixed conditions are specified, respectively, and a is a constant. Consider the energy product of (5.11). Since the variable u and x are in the f i e l d of definition of operator A, they must satisfy a l l the boundary conditions in equations (5.15). The energy product of two elements w A and A from D, is - - A = IQ ^-T-*~ + I T - - " d i s -integration by parts of the f i r s t two terms on the right hand side yields 82 [A,A] = - $ x u>nds + $ f u»nds - <$ x u-sds + <P x u-sds - - 1 nn- - i J s nn- - ; g ns- - J g ns- -+ /{j [HTT*! + ITI2 ~ ITCxJdQ. (5.16) Since A and A satisfy the same boundary conditions, the boundary integrals above cancel each other. Further, since C is symmetric, therefore T T f Cx = T Cf and [A.A] = IQ t ^ T I * l + l T T a - xTCJf]dfi or [A,A] = [A,A] which proves that the operator A is symmetric. Therefore the energy product W in the space H is given by JTL. [A,A]. = / [2xTTu - xTCx]dfi + / ctuTuds (5.17) M since the contribution to the line integral arises only from the S^ part of the boundary. In this case, the mixed boundary conditions give rise to a line integral i n the energy product. The element matrices can be obtained by substitution of the approximate f i n i t e element solution into (5.17), which would s t i l l be symmetric. However, some of the zero entries of the matrix (4.19) are now replaced by the boundary contribution from the non-zero line integral. The energy convergence would be ensured i f the coordinate functions W are complete in H . Further, since the boundary conditions prescribed on S and S^ contain stress terms they are therefore natural. The coordinate functions then need not satisfy these boundary conditions whereas the homogeneous kinematic boundary conditions on S are enforced in the same manner as discussed in section 5.1. 83 5.3 Nonhomogeneous Boundary Conditions The nonhomogeneous boundary conditions for a plane stress problem are usually prescribed as 0 u, = on S (5.18a) I i u 0 x..n. = T. on S_ (5.18b) i j J i T 0 x..n. + au. = C. on S„„. (5.18c) i] J i i M These can be incorporated by changing variables in such a way that the problem reduces to one with homogeneous boundary conditions. Thus, T assume that there exists functions A'=<u', T'> , i.e. functions u! and x x.'. , both continuous, such that i j u! = u? on S x x u T i j n j = T i ° n ST ( 5 - 1 9 ) T ! . n . + au! = C? on S„. i] ] x x M Define new variables u 2 and x 2. in the following way ± ± 3 u'. X X . . * . . 13 i3 13 (5.20) or A" = A - A' . (5.20a) Substituting equations (5.20) into (5.2), (5.3) and (5.18) yields -x'.'. . = f. + x'. . . = f.1 (5.21) 13 »J i 13 ,3 i u1.' = 0 on S X u x^n. =0 on S_ (5.23) T'.'.n. + au*.f = 0 on S.„. 13 3 i M 84 Thus the problem in A" has homogeneous boundary conditions, and the only difference in (5.2), (5.3) and (5.21), (5.22), respectively, is the introduction of the term T \ . . on the right hand side of (5.21) and [-1/2 I J >J (u! ,+u! .)+C... , T ' , ] in (5.22). If these terms were known, i t would be i , j j , i l j k l k l j possible to obtain an approximate solution for A" and convergence would be w ensured i f the coordinate functions were complete in the space and the associated energy product as given in equation (5.17). — • — T Let the f i n i t e element approximations for A"= <u" v"> be M k-l cb.u" ; 1=1,2. k=l 1 l k (5.24) or where i j A" T -M E I | ) . . T V . , ; i , j = l , 2 . I = 1 i J i J l T *N N <P w -M -M .k ,1 (5.25) (5.26) and N -M <u'' T ' . * . > " ik 131 Then the required element matrix equation, from section 4.2, is . * M ] A = (P ,*M) here (5.27) T P = f.' 0 1 0 s« Thus a solution for A" can be obtained such that |A" - A" I -> 0, as M and N ». (5.28) It i s worth noting that here the kinematic nonhomogeneous boundary conditions are conveniently incorporated through the load vector. Such a procedure cannot be achieved in the displacement approach because i t would 85 require displacements to satisfy the nonhomogeneous boundary conditions and have continuous f i r s t derivatives. However, a solution may be obtained directly for u. and T,. 1 i j without actually knowing A1. This is accomplished by assuming approximations for u. and x.. as M -k u. = X i> .u., ; i=l,2. I , , I ik k=l (5.29) N -2 T i i = 2t± V i J i J 1 ' j ' = 1 ' 2 ' A = $ w -M -M T TN 7k -1 where $ w = <d>, \b . . > -M T i T i j -N T and w_, = <u.. x. > . -M ik i j I The sets of functions {<^ } and {ij;^ } are complete in H^ space considered with respect to the nonhomogeneous boundary conditions. Therefore as the function A"+A' satisfies such nonhomogeneous boundary conditions, i t follows that there exist u., and x.., such that ik i j 1 |A - (A"+A')| -> 0, as M and N + ~. (5.30) Now the function A-(A"+A') has homogeneous forced boundary conditions; N therefore the energy product of i t with any arbitrary function 0 ^ in the space w H^ , with homogeneous forced boundary conditions, would also vanish in lieu of (5.30). Therefore IA - (A"+A'), 0 ^ J A = 0. (5.31) Since 0^ is arbitrary, i t can be replaced by 0^ in (5.26) and equation (5.31) can be written as 86 [A - ( A , , + A ' ) , * J J ] A = 0. (5.32) From linearity of the operator A, (5.33) and substituting this into equation (5.27) yields (5.34) But or / T A r (2 • •M> " /ft fV 0 1 o % l . r ( f ^ k ) ( g M . *"?.), + i >dft. d n Therefore from (5.21) /ft /ft <V + T i_,_>*_ d n and integrating by parts the second term on the right hand side yields / n ( f i ' < ) d f i - /ft ( f_*_ - T i j ^ f j > d n +* T _ J n J * i d s -Using equations (5.19) and the fact that <Ju=0 on gives /ft ( f i » ^ d n = /ft < fi*_ - ^ i , ^ + / s _ v _ d s + / s ( c _ - « u _ > * i d s M (5.35) and from (5.22) /ft <*__*__ > d n = /ft { - 1 / 2 ( u i , j + u j , i > + c u k _ T k i } * q i . d n - ( 5 ' 3 6 ) Therefore <ET.»S> = </, + / S t V _ d * + J C^ kds) - [ A ' , $ A (5.37) T M because, from (5.17) 87 / . T ! .cb k .dft + f au : cp kd M . J {i / 2 ( u ! . + u! .) - c . . . . T' }t|>?.dn. (5.38) Now equation (5.37) enables (5.34) to be written as 0 L- . 0 X f±*lda + /sT V i d s + U V i d s T M (5.39) 0 — — — T which are the mixed Galerkin equations that govern a s o l u t i o n f o r A=<u x> w i n when the boundary conditions are nonhomogeneous. Writing out i n f u l l the equation (5.39) k V i , j d f i + j s M a V i d s = k f ± < d n + u T>J D S + /s c i * i d s ' M T M M k=l , 2 , . . . N. M. (5.40) (5.41) The equations (5.40) and (5.41) require that the approximate s o l u t i o n f o r u^ given by (5.29) should s a t i s f y the nonhomogeneous forced boundary conditions and the coordinate functions <j> the homogeneous forced boundary conditions on S . u In the f i n i t e element method i t i s not necessary to introduce N -N d i f f e r e n t approximations A corresponding to <i> and * coordinate functions. - . -M -M N The coordinate functions of a f i n i t e element approximation associated with the degrees of freedom that do not l i e on s a t i s f y homogeneous condi-tions on S^, i . e . vanish on S^. Therefore the equations (5.40) and (5.41) can be solved by spec i f y i n g the values of u ^ of (5 .29) , the nodal degrees of freedom on S^. The remaining coordinate functions s a t i s f y homogeneous forced boundary conditions and thus only one set of coordinate functions need be introduced. 88 To show that the approximate solution A of (5.29) converges in the energy norm, rewrite (5.31) from linearity of the operator A in the form [_>§M]A = t A " + A , , ^ J ] A . (5.42) Also from equation (5.20a) [_>_M]A = [ ^ ' ^ M ] A ' ( 5 ' 4 3 ) Subtracting (5.42) from (5.43) and from linearity of the operator A; Using property ( i i i ) of theorem 3.4.2 (Schwarz inequality) But from (5.28), |A"-A"|-K); as M and N^; therefore [A-A,e*J]A + 0, as M and N-*». - - -MA N Since 0 ^ is an arbitrary function with homogeneous forced boundary conditions, i t may be set equal to A-A. Then U-A.A-A] 0, as M and N-*», which implies that l_~_lA °> a s M a n d N~*"- (5.44) That i s , the approximate solution to the nonhomogeneous boundary condition problem converges i n the energy norm of H^. One of the advantages that the mixed method offers li e s i n different ways of incorporating the boundary conditions. So far the natural boundary conditions have been associated with stresses, a consequence of extracting the boundary integrals from the equilibrium equations. This led to constrain-ing of forced boundary conditions on u on S^ through the nodal variables. It w i l l be demonstrated here that in fact a l l nonhomogeneous 89 boundary conditions can be incorporated through boundary integrals from both equilibrium and constitutive-kinematic equations through the Hellinger-Reissner's mixed variational principle since i t yields the same equations as the mixed Galerkin method for T^T theory (Chapter 3). The Hellinger-Reissner mixed variational principle for plane stress with unit thickness and zero body forces f and f can be written as x y I = /_ [ T U , + T v, +T (u, +v, ) - - ^ x 2 +x2 - 2 V T x +2(1+V ) T 2 } ] dft Ju xx x yy y xy 'y 'x 2E xx yy xx yy xy - J- (up +vp" )ds - / [(u-u)p +(v-v)p Ids (5.45) o 2 x y o 2 x y where p . = T . . n . . S J is the portion of the boundary S on which the stresses T . . or p. are prescribed. The part So has the displacements u. (u and v) prescribed on i t . The stress boundary conditions of (5.18a) and (5.18b) on S „ and S „ can be considered to be on part S i while S would coincide with T M u S 2 . Therefore p. = x..n. = T° on S, i i j J i > p. = x n. = C? - au. on S_, i i ] J i i M (5.46) and the f i r s t variation of I with respect to u. and x.. gives & 1 - /fl [ \ j & U ± , 3 ] d n + /fl ^ ^ ^ I J ^ . ^ ^ i j k i ^ 1 ^ 1 1 " " jST T° (5.47) 6 V S " /sM < C!hV 6 ui " /Su ( V U S ) 6 ( T l j B J ) d 8 " J S u T U n j f i u ± " °-Now assume approximations for u_^ and i as M k u. = £ «J).u ; 1=1,2. (5.48) 1 k=l 1 l k N I and x = E * , , T ; i,j=l,2. (5.49) .^J .2=1 90 Therefore M k - N 1 6u. = Z <t>.6uM ; fix.. = E I|>..6T.,., (5.50) 1 k=l 1 1 3 1=1 1 J i j N and fiu. . = E d> .6u., . (5.51) i, J k = 1 i . J ik Substitution of (5.48) to (5.51) into (5.47) yields 5 1 " j x l / n V i , i d ^ s _ TS*ids"/sM ( C i " ^ i ^ i d s - / s u ^ i jV i d s ] 6 u i k + qL { 1 / 2 ( ^ i , J + ^ , i ) - C i J k i\i }*i j d f i -/s u " u / i j n j d s + / s u u _ * i j n j d 8 ] fix..q = 0. (5.52) Since 61 vanishes for arbitrary variations 6u., and fix.. , the following ik i j q ' equations are obtained: k V i , j d n + 'sM a V_ d s ~ j s ^ V i d s = J s T T?*i d s + 'sM c i * i d s ; M u J J T M k-1,2, . . . M. (5.53) q-1,2, . . . N. (5.54) The equations (5.53) and (5.54) contain 2M+3N equations for 2M+3N unknowns with a symmetric matrix of coefficients. It is interesting to note that the nonhomogeneous forced boundary conditions are applied through the displace-ment vector in (5.54) and hence need not be constrained as was done in the previous case. Except for boundary integrals over these equations are exactly the same as (5.40) and (5.41). 91 5.4 Boundary Residual Concept An alternate procedure for incorporating boundary conditions in the mixed Galerkin method is similar to the boundary residual concept presented by Finlayson and Scriven [8] and Finlayson [7] in which the domain residual together with boundary residuals are made orthogonal to the shape functions of the approximate solution. Consider the plane stress linear e l a s t i c i t y problem of equation (5.9) with nonhomogeneous boundary conditions of equations (5.18) and approximations for u_^ and i of equations (5.48) and (5.49). The substitution of u^ and x into the f i e l d equations (5.9) and boundary conditions (5.18) yields the following residuals: R e i = [-\i,ff±] i n f i (5'55) R c i = [ 1 / 2 ^ , i + ~ U i A ) - % ^ l \ l ] ± n B °-56) R u i " " [ V U i J o n Su ( 5 - 5 7 ) *r±= [ V r T i ] O N S T ( 5 - 5 8 ) R... = [ T . .n.+au.-C1?] on S„ (5.59) Mi l j ] i i M where R . and R . are the domain residuals for the equilibrium and kinematic-el c i constitutive equations, respectively; R^ , R^, and R^ are the boundary residuals on S , S m and S w, respectively, u T M If the residuals R^ from equation (5.55) along with R^ and R^ from equations (5.58) and (5.59) are made orthogonal to the shape functions for u_^ and residuals R^ from (5.56) along with residual R^ from (5.57) to the shape functions for x „ , the following equations result: k ^ i j f j - f i ^ i d n + l s S v T i ) d s + Un <Vj + ° v ci>*i d s • 0 ; k=l,2, . . . M. (5.60) 92 /ft [1/2('Z±,3+'U3,± ) ~ W i j ^ V " " ^ <VUI>*i_ q=l,2, N. (5.61) Assuming f.=0; and applying Gauss' theorem to (5.60) yields M k=l,2, M. (5.62) /ft [ 1 /« ui, j 4» J li)"Vu ]* ,/^s u ..*?.n.ds = L u q=l,2, N. (5.63) The equations (5.62) and (5.63) are the same as equations (5.53) and (5.54) in the previous section. Therefore in linear e l a s t i c i t y , the equations obtained by the mixed Galerkin method, the mixed variational principle and the boundary residual concept are the same. In the displacement approach, Hutton [15] showed that the equations obtained from the Galerkin Method and the boundary residual concept for approximations from wider class would be identical i f the forced boundary conditions were either homogeneous or were identically satisfied by the f i n i t e element approximations when nonhomogeneous. However, the f l e x i b i l i t y offered by the mixed methods in incorporating the boundary conditions, forced or natural, provides a wider equivalence. 93 CHAPTER 6 EIGENVALUE ANALYSIS OF THE ELEMENT MATRIX An eigenvalue-eigenvector analysis of the element matrix arising from the mixed method i s presented in this chapter. Various combinations of displacement and stress approximations over a triangular and a rectangular element are considered. Two problems are included, namely the linear elasticity plane stress and the linear part of the Navier-Stokes equations. It i s anticipated that the analysis of eigenvalues w i l l provide some insight as to choice of approximations for the dependent variables involved so that completeness is achieved. 6.1 Linear Elasticity Problem Consider the matrix equation (5.9) for the plane stress problem with zero body forces, i.e. f =f =0, J x y 0 0 3 ~3x 0 3 "3y u 0 0 0 0 3 "3y 3 ~3x V 0 JL 3x 0 1 "E v_ E 0 T XX = 0 0 3 3y E 1 ~E 0 T yy 0 3 3 0 0 " -2(l+v) 0 dy 3x E T (6.1) The variational principle for (6.1) with homogeneous boundary conditions (section 5.1) can be written as I = /_ [x u, +T v, +T (u, +v, )-^{T2 +T2 -2VT T +2(1+V)T2 }]dn. (6.2) J Q xx 'x yy y xy y x 2E xx yy xx yy xy Since I represents strain energy, inspection of the right hand side in (6.2) suggests the following three rigid body modes which yield zero strain energy: 94 T = T T = X X yy xy T = T = T = X X yy xy T = T = T X X yy xy (i) u = constant, v = 0 ; = = 0 xx ( i i ) u = 0 , v = constant; x ^ = X t t = x^_ = 0 (6.3) ( i i i ) u = -cy , v = cx These r i g i d body modes are expected to be removed by the specified kinematic boundary conditions. Furthermore, i t is required by the functional of (6.2) that the displacements satisfy the kinematic boundary conditions while the stresses emerge as natural boundary conditions. Therefore the fi n i t e element approximations should be in compliance with these requirements, i.e. the matrix of equations (4.19) should exhibit the ri g i d body modes of (6.3) . The independently chosen approximations for the stresses and the displacements have to comply with the completeness requirement (i) of section 4.4. The mean convergence of strains from the assumed stresses to the strains derived from the assumed displacements would be assured for a fi n i t e number of degrees of freedom i f the former contains a l l the strain modes and perhaps more than the strain modes in the latter. It is assumed that the displacements possess a l l the rigid body and constant strain modes and that the stresses possess a l l the constant stress modes. It is now asserted that the violation of the completeness requirement (i) results in a hypersingular element matrix, i.e. the number of zero eigenvalues greater than the rigi d body modes expected in a problem. The eigenvectors for the extra zero eigenvalues correspond to mechanisms which are defined as the kinematic freedoms possible when the material has no elastic stiffness. This is illustrated by the following example. e * Consider one element domain fi and the approximate solutions for u,v,x ,x and x as xx yy xy 95 u = ax + by + c x 2 + dxy + ey 2 v = bx + fy + gx 2 + hxy + i y 2 x x x = 3 (6.4) x = k yy T = I. xy The polynomials f o r u and v are so chosen that the r i g i d body modes have been eliminated by s a t i s f y i n g the kinematic boundary conditions. Therefore e = u, = a + 2cx + dy xx x £ = v, = f + hx + 2iy (6.5) yy x Y = u, + v, = 2b + (d+2g)x + (2e+h)y. xy y x From the equations (6.5) the s t r a i n s derived from u and v are complete l i n e a r polynomials. Therefore the mean convergence of constant s t r a i n s from x , x and T i n (6.4) to the s t r a i n s i n (6.5) would not occur and the com-yy xy • pleteness requirement ( i ) i s v i o l a t e d . The parameters i n (6.4) are to be determined from the v a r i a t i o n a l formulation. The s u b s t i t u t i o n of the ex-pressions i n equations (6.4) into (6.2) y i e l d s I = [j(Aa+2ac+Bd) + k(Af+ah+2gi) + 1(2Ab+ad+2ag+23e+3h)] A 2E [j2+k 2-2vjk+2(l+v)I 2] (6.6) vhere A = / dfi, a = / xdft and 3 = / ydft. ft6 ft6 QE rhe system of equations governing the one element domain i s obtained by making 1 stationary with respect to the unknowns a,b,c,d, . . . I. This i s 96 0 0 0 0 0 0 0 0 1 A 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 2a 0 0 0 0 0 0 1 e 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 symmetric 0 1 0 0 0 0 a 0 28 A 0 0 2a A vA "E E 0 -I « -2(l+v)A (6.7) J k I It is observed that the f i r s t and the third rows; second, f i f t h and the seventh rows; and sixth and the ninth rows are the same except for some multiples, while the fourth and the eighth rows are linear combinations of the f i r s t and the second rows and the second and the sixth rows, respectively, rherefore only six of the twelve equations in (6.7) are linearly independent. Hence the element matrix has a rank of 6 instead of 12 and therefore is singular. As a consequence c,d,e,g,h and i are indeterminate and these are the coefficients of the quadratic terms i n the polynomials for u and v, equations (6.4). For a matrix of the form (4.19), 0 (m+n)x(m+n) mxm T a nxm mxn - b nxn (6.8) 97 w h i c h i s r e a l a n d s y m m e t r i c , t h e e i g e n v a l u e s a r e r e a l a n d b e c a u s e o f t h e i n d e f i n i t e n e s s o f t h e m a t r i x t h e s e c a n b e e i t h e r p o s i t i v e , z e r o o r n e g a t i v e . F u r t h e r , i f t h e n x n b s u b m a t r i x i s p o s i t i v e d e f i n i t e a n d r t h e r a n k o f t h e m a t r i x S, i t i s s h o w n i n A p p e n d i x B t h a t t h e e i g e n v a l u e d i s t r i b u t i o n i s o f t h e f o l l o w i n g t y p e : ( i ) m p o s i t i v e a n d n n e g a t i v e e i g e n v a l u e s i f r=m+n; ( i i ) ( r - n ) p o s i t i v e , ( m + n - r ) z e r o a n d n n e g a t i v e e i g e n v a l u e s i f r<m+n. T h u s t h e m a t r i x o f e q u a t i o n ( 6 . 7 ) w h i c h h a s r = 6 , m=9, n = 3 , w o u l d y i e l d t h r e e p o s i t i v e , s i x z e r o a n d t h r e e n e g a t i v e e i g e n v a l u e s . A l s o f r o m t h e f u n c t i o n a l i n e q u a t i o n ( 6 . 6 ) w i t h i n d e t e r m i n a t e c , d , e , g , h a n d i , f o r e i g e n v a l u e s t o b e z e r o , t h e s t r e s s e s j , k a n d I m u s t b e z e r o . I t i s , t h e r e f o r e , o b v i o u s t h a t w i t h t h e d i s p l a c e m e n t s a n d s t r e s s e s o f e q u a t i o n s ( 6 . 4 ) , t h e e l e m e n t c a n s t r a i n w i t h z e r o s t r e s s e s ( i . e . f o r m i n g m e c h a n i s m s ) . T h e n o n - s t r e s s i n g s t r a i n m o d e s a r i s e f r o m t h e q u a d r a t i c t e r m s i n u a n d v g i v i n g r i s e t o t h e l i n e a r t e r m s i n s t r a i n s , e q u a t i o n s ( 6 . 5 ) , w h i c h a r e n o t c o n t a i n e d b y t h e a s s u m e d s t r e s s e s . T h i s v i o l a t e s t h e c o m p l e t e n e s s r e q u i r e m e n t ( i ) a n d i t i s t h i s v i o l a t i o n w h i c h l e a d s t o m e c h a n i s m s . F u r t h e r t h e n u m b e r o f m e c h a n i s m s c o r r e s p o n d s t o t h e n u m b e r o f t e r m s p r e s e n t i n t h e d e r i v e d s t r a i n s , e q u a t i o n s ( 6 . 5 ) , b u t w h i c h a r e n o t i n c l u d e d i n t h e s t r a i n s o b t a i n e d f r o m t h e a s s u m e d s t r e s s e s u s i n g t h e c o n s t i t u t i v e l a w s . I t w a s s h o w n i n t h e p r o o f o f t h e o r e m 4 . 5 . 1 t h a t t h e e l i m i n a t i o n o f s t r e s s e s f r o m t h e m a t r i x e q u a t i o n ( 4 . 1 8 ) , i . e . "o e i a r en u r en P eT a - b 6 e T — 0 ( 6 . 9 ) g i v e s -1 T a e b a 6 u 6 = kV ( 6 . 1 0 ) 98 where the matrix k e is identical to the element stiffness matrix one would obtain from the displacement f i n i t e element method using the same approxima-ting polynomials for the displacements, provided the assumed displacements and stresses in the mixed method were complete. In the example considered here, obviously the completeness is violated. However, the static condensation of the matrix equation (6.7) by eliminating stress degrees of freedom j,k and I is performed to obtain the following matrix equation: 1 - V A 0 2a 3 0 vA 0 va 2vg (l-v)A 0 (l-v)a 2(l-v)8 0 2(l-v)a (l-v)B 0 4a 2 .2a3 n n 2va 2 4vag 0 2va 0 B 2 , ( l - v ) a 2 (l-v)qg 0 ( l - v ) a 2 (l+v)a3 2v3 2 A 2A ~ v 6 A " A 2A A 2(l-v)3 2 0 2(l-v)o3 (l-v)3 2 symmetric A A A O a 2(l-v)a 2 (l-v)a3 0 23 0 = 0. (6.11) 99 The very same linear dependence of equations is observed in (6.11) as was for (6.7) and again c,d,e,g,h and i are indeterminate. The rank of the 9x9 matrix is 3 instead of 9. Therefore i t cannot be the same matrix as one would obtain from the displacement approach using the same approxima-tionsfor u and v as in equations (6.4), since i t is well known that the stiffness matrix i s positive definite after the rigid body modes have been eliminated. Now i f the quadratic terms in u and v were dropped, then the completeness requirement (i) is satisfied thus resulting in the follow-ing non-singular equation with rank 6. 0 0 0 A 0 0 0 0 0 0 0 2A A 0 0 0 0 0 A o o . | ^ o 0 0 A f i 0 -2(l+v)A a b f 3 k 1 = 0. (6.12) 0 2A 0 0 0 The matrix of equation (6.12) would yield three negative and three positive eigenvalues. The static condensation of the matrix by elimination of the stress degrees of freedom j,k and 1 gives l-v< A 0 vA 0 2(l-v)A 0 vA O A = 0. (6.12a) The matrix of equation (6.12a) is exactly the same as the stiffnex matrix one would obtain from the displacement approach for constant stress triangles. To further demonstrate what has been explained so far, various combinations of interpolations for the displacements and the stresses over a triangular and a rectangular element are considered. The typical node 100 numbering and nodal degrees of freedom for both triangular and rectangular elements are shown in Figures 4 and 5. The derivation of the element matrices is demonstrated in Appendix A. An eigenvalue routine, which uses Householder transformations, is then used to solve the following eigenvalue problem. - X = 0. (6.13) T T T T T T T where 6 =<u v > and T =<T T T > and [S]= - - - - -xx -yy -xy 0 a T z a b as i n (6.8), and [I] is the identity matrix. The qualitative description of the eigenvalues and the composition of the eigenvectors for a l l the combinations of interpolations used for the displacements and stresses over a triangular element appears in Table I and that for a rectangular element is l i s t e d i n Table II. For both triangu-lar and rectangular elements, the number of negative eigenvalues corresponded to the number of stress degrees of freedom, whereas the three zero eigen-values for th_ expected rig i d body modes are obtained only for the displace-ment stress combinations which satisfy the completeness requirement ( i ) . In the cases where more than three zero eigenvalues are obtained, the number of extra zeroes corresponded to the number of modes present i n the strains derived from the assumed displacements that were not contained in the strains from the assumed stresses. The eigenvectors are composed of the same dis-tribution as the assumed approximations for displacements and stresses i n a l l cases except for rigid body modes where u and v satisfied u,x=0, v, =0, u, +v, =0 and stresses were zero, y x It is essential for convergence in the energy sense that a mixed fin i t e element formulation conform to the completeness requirements of 1 0 1 section 4.4, and be able to represent r i g i d body modes and the constant s t r a i n s (with completeness requirement (i) s a t i s f i e d and assumed dis p l a c e -ments that possess constant s t r a i n s , i t i s implied that the corresponding assumed stresses would contain the constant s t r e s s e s ) . To include r i g i d body modes and constant s t r a i n s i n the assumed displacements i s a simple matter. However for c e r t a i n combinations of assumed displacements and stresses, i t i s not quite obvious that completeness i s achieved, e s p e c i a l l y for incomplete polynomials. For example using biquadratic displacements and b i l i n e a r stresses over a rectangular element y i e l d s four zero eigen-values. A scheme to check completeness requirement (i) and trace the terms i n the polynomials used f o r the displacements which correspond to s t r a i n s not included i n the assumed stresses i s presented here. The poly-nomials considered are the ones mentioned i n the example above. The assumed biquadratic displacements u and v are u = ai+a2X+a3y+at +x 2+a5xy+a5y 2+a7X 2y+a3xy 2 v = bi+b2X+b3y+bi +x +b5Xy+bgy +byx^y+b8xy z and the b i l i n e a r stresses are Tx x = C! + c 2 x + c 3 y + c 4xy (6.14) T = d 1 + d 2x + d 3y + d^xy (6.15) x x y = e! + e 2x + e 3y + ei+xy. The s t r a i n s derived from u and v are e = u, = a2+2aitx+a5y+2a7xy+agy2 (6.16a) XX X e y y = v, = b 3+b 5x+2b 6y+b 7x 2+2b 8xy (6.16b) Y„„ = u, +v = (a3+b2)+(a5+2b 1 +)x+(2a 6+b 5)y+a 7x 2+2(a 8+b 7)xy+b 8y 2. (6.16c) xy y x 102 Now using the generalized Hooke's Law for the plane stress problem, the strains corresponding to the assumed stresses in (6.15) are Sex = E[Txx~V^yy-' = E [ ( cl- v dl)+(c2-vd 2)x+(c3-vd 3)y+(c l t-vdt t)xy] (6.17a) Sy = E 1"~ V^xx + Xyy ] = f [ ( d l - v c l ) + ( d 2 - v c 2 ) x + ( d 3 - v c 3 ) y + ( d ' t - v c i + ) x y ] (6.17b) 2(l+v) r 2(l+v) r , . . . fe. i t \ Yxy = — E Txy = — E Ui+e^egy+e^xy]. (6.17c) If the displacements were assumed to be a complete polynomial of degree p and stresses to be a complete polynomial of degree (p-1) , then completeness would be achieved. However, the polynomials considered here are not complete and comparison of equations (6.16) with (6.17) shows that certain terms in the derived strains are not contained in the corresponding strains from the assumed stresses, i.e. a 8y 2 in (6.16a), b 7x 2 (6.16b), a 7x 2 and b 8y 2 in (6.16c). But a 7 and b 8 as coefficients of the bilinear terms in the derived strains e and e match with the coefficients of the bilinear terms in e and e xx yy xx yy of (6.17), respectively; while (a 8+a 7) appearing as coefficient of the bilinear term in the derived shear strain Y matches with the coefficient xy of the bilinear term in Y of (6.17). Therefore only one of the coefficients xy a 7,a 8,b 7 and b 8 is indeterminate, hence i t leads to only one mechanism besides the three rigid body modes. This is confirmed by the results obtained for biquadratic u and v and bilinear stresses over a rectangular element, Table II. For the same combination, the mode shape for the mechanism after elimination of the rigid body modes, is illustrated in Figure 6. The static condensation by elimination of stresses was performed for a l l of the combinations of displacements and stresses appearing in Tables I and II. The condensed element matrix is found to be exactly the same as the stiffness matrix one would obtain from the displacement method using identical assumed displacements over an element except for the combinations 103 where the completeness requirement (i) i s violated which possess the same number of mechanisms. 6.2 Linear Part of the Navier-Stokes Equations The basic equations governing the two dimensional, steady, incom-pressible flow are the well-known Navier-Stokes equations p(uu, + vu, ) + p, - uV2u =0 in tt (6.18) x y x p(uv, + w, ) + p, - yV 2v = 0 in tt (6.19) x y y u, + v, = 0 in tt (6.20) x y where u,v are the x,y components of velocity, respectively, p is the f l u i d density, p, the pressure, u , the dynamic viscosity, and tt the open domain of the problem. In terms of deviatoric stresses directly, the equations (6.18) and (6.19) can be written as p(uu, + vu, ) + (p, - T - T ) = 0 in tt (6.21) x 'y r'x xx,x xy,y p(uv, + w, ) + (p, - T — T ) = 0 i n tt (6.22) x 'y y xy,x yy»y where T and T are the normal deviatoric stresses in the x and y directions, xx yy respectively; and T is the shear stress. The equations relating deviatoric stresses to velocities for a Newtonian flu i d are Txx = - 3^(u'x + V + 2^ U'x ( 6 ' 2 3 ) Tyy = " 3 ^ ' x + v » y > + 2 y V ' y ( 6 " 2 4 ) x = p(u, + v, ). (6.25) xy y x These equations can be put into an alternate form by solving for the velocity gradients as 104 1. 1 , = — ( T - -T ) u xx 2 yy v y 2 xx yy in ft in ft (6.26) U , + V , = —T y x y xy in ft. Now the equations (6.21), (6.22), (6.20) and (6.26) can be put into the matrix operator form for mixed formulation as t 9 _, 9 \ p ^ T x ^ o f - -9_ _3_ 3x 9x 0 8_ "3y 3 3x 3_ 3x 9_ 3y /> 9 x 9 \ p ( V x ^ 3 3y 0 3_ 3y 3_ 3x 3_ 3y .L. J _ 3y 3x 1 M 1 2y ~ 0 2y ~y 0 0 0 1 y X X yy xy = 0 in ft. (6.27) For steady creeping flow, a special case of incompressible, steady Newtonian flow, the nonlinear part of the matrix operator above, which makes i t non-symmetric, can be dropped. Thus the following f i r s t order, linear differen-t i a l equations result and involve a symmetric matrix d i f f e r e n t i a l operator 0 0 0 0 9_ 3_ "3x ~9y 0 9_ _9_ 3x ~9x f 0 3y 0 0 0 - — 3y _3_ _3_ "3y 9x 0 1 0 1_ 2y 0 0 v X X = 0 i n ft; (6.28) 105 comprising six equations for six unknowns. Here, in lieu of T*T theory 1* = -T = 9-- 9- 0 - a -3x 3x 3y 1- o 3y 3y 3x T T and for u=<u v> and x=<p x x x > , the analogous form'of equations xx yy xy (6.28) to the plane stress elasticity equations (5.10) is " 0 u = 0 T -F X or AA = 0 where F = ' 0 0 0 0 " 0 1 u 2y 0 -1_ 1 2u y 0 0 0 0 1 y (6.29) (6.29a) (6.30) which i s a positive semidefinite matrix. In tensor notation, the equations (6.28) take the form p,. - x.. . = 0; i=j=l,2; in ft i I J ,2 u. . = 0; i=l,2; i n 2 ( u i . i + U i . i ) - C n i T k i = °i i=J=k=I=l,2; in ft 2 v u i , j T U j , i ; " " i j k l l k l and can be subjected to some boundary conditions analogous to equations (5.18); on S (6.31a) (6.31b) (6.31c) u. = uV l l (-P<$. . + x. .)n. _ T 0 on S„ (6.32) (-p<5. . + x. .)n. + au. = C° i j i j 2 i I on S. M 106 where S +Sm+S =S, the boundary of domain fi. u T M J The matrix equation (6.28) is similar to matrix equation (6.1) except for pressure p, the continuity equation as zero divergence of velocities and the deviatoric normal stresses. The mixed variational principle for these equations for homogeneous boundary conditions can be derived as I = /_. [-p(u, + V , ) + X U, + T V, + T (U, + V , ) ; fi x y xx x yy y xy y x - ~ { T 2 + T 2 - T T + T 2 }]dfi (6.33) 2u xx yy xx yy xy which is also similar to the mixed variational principle for the linear elasticity plane stress problem in (6.2) except for the term p(u,x+v,^) and requires velocities to satisfy the kinematic boundary conditions. The variational principle of (6.33) also gives the three rigid body modes as in (6.3); (i) u = constant, v = 0; x = x X X ( i i ) u = 0, v = constant; x = x = x =0 (6.34) xx ( i i i ) u = -cy; v = cx; = x while pressure can be arbitrary, since the incompressibility leads to the divergence of u and v to be zero rather than be proportional to pressure. If the f i n i t e element approximations for the variables involved in (6.33) are chosen over an element domain fi as u = E d> .u. 1=1 (6.35a) m v = E cb. v. X = 0 yy xy X 0 yy xy X = 0 yy xy 107 _ n x = E ib. x XX . T l XXI 1=1 T = E \b.T yy ' i = 1 i yyi (6.35c) T = E \b ,T x y 1=1 x y Then the substitution into the functional I of (6.33) and setting i t s f i r s t variation 61 to zero for stationarity yields the following matrix equation: 0 0 e a 0 0 f 0 0 0 0 0 u-0 f , 2y-0 0 0 Ol b " u b a V Ol 0 p U 2y-0 T -XX - i d y-0 T -yy Ol - i d y- T -xy = 0. (6.36) Here u=<uj u 2 . . . u^ > , v=<vi v 2 v m> , P=<Pi P2 T =<T ,T „ -xx xxi xx2 T > , etc., are the linear vectors of nodal degrees xxn ° of freedom. The submatrices a,b,d,e and f are obtained in the following manner: 3 i j " / ne*i,x*J b i j " /ne+i.y^ dn 1=1,2, e 3=1,2, dVi i 5 j = l , 2 , d i i = / n e * i * j e j=l,2, . f i j = -la**i,y!U dQ, m; n. . n. m; 1. (6.37) 108 Note that the matrix d is symmetric and positive definite. The matrix of equation (6.36) is symmetric, therefore the eigenvalues of the matrix are real. Further, i t is indefinite and i f the rank of the matrix is (2m+I+3n), then from Appendix B i t has 2m positive and (I+3n) negative eigenvalues. The choice of polynomials for u,v,x , XX x and s t i l l has to comply with the same completeness requirements set out in the previous section. However the requirements on the pressure f i e l d are dubious because i t is not related to the strains. The obvious question that arises here i s : what are the completeness requirements on p? Further i t is not possible to ask for mean convergence of pressure to the volumetric strain as was done for the stresses in the plane stress linear elasticity problem because of the incompressibility constraint. In the f i n i t e element application to the Navier-Stokes equations (6.18) to (6.20) using the primitive dependent variables u,v and p, a similar situation was faced by Taylor and Hood [34] and Olson and Tuann [26], One of the possible variational principles for the linear part of these equations used i n the reference [26] i s J(U'V'P) " / f t [ ^ ' x S ^ » y + V » / 1 " p ( U ' x + V ' y ) ] d - / (Xu + Yv)ds (6.38) T, where (X,Y) are the specified traction on the boundary S^ _. The term P( u> x + v,^)dft appears in both variational principles, I of (6.33) and J(u,v,p) of (6.38). Therefore the requirement that the pressure interpolation should be one degree less than those for the velocity components, as found by Olson and Tuann [26], is also expected here and indeed confirmed by the numerical results. A different explanation for such a requirement is pre-sented here. 109 Since the divergence of u and v i n equation (6.20) i s zero, i t acts as a constraint equation. If the approximations used for u and v involve complete polynomials of degree m, then u,^ and v, would involve complete polynomials of degree (m-1). Therefore the unknowns corresponding to the (m-1) complete polynomial for u,^ must be r e l a t e d to the unknowns corresponding to (m-1) complete polynomial f o r v, i n order to s a t i s f y continuity i n the discr e t e sense. Assuming that U,V,T ,T and x have been chosen properly, xx yy xy r r i . e . do not y i e l d any mechanisms except for the r i g i d body modes, then the d i s c r e t i z e d continuity equations < e T f T > j " } = 0 (6.39) should have a rank not less than m+"'"C2, i . e . combinations of (m+1) taken 2 at a time. However i n the f i n i t e element formulation, equation (6.36), the number of d i s c r e t i z e d continuity equations i s associated with the number of degrees of freedom for pressure thus l i m i t i n g i t to m+~'"C2. As a consequence, the pressure should have degrees of freedom not more than C2y which implies a complete polynomial of degree not higher than (m-1). Consider the following schematic representation of the complete polynomials for u,v and p. V e l o c i t y u 1 Constant a\ x y Linear a2 a 3 x 2 xy y 2 Quadratic a^ a$ ag x 3 x 2y xy 2 y 3 Cubic a 7 a 8 ag a 1 0 x 4 x 3y x 2 y 2 xy 3 y 4 Quartic a n a\z al3 alh a15 etc. etc. etc. (a) 110 Thus i f the a.'s are the coefficients of the polynomial for u, the complete quadratic for u can be written as u = ai + a 2x + a 3y + a^x2 + a$y.y + agy 2. (6.40) Replacing a.'s by b^'s in scheme (a) for v; the complete quadratic for v is v = hi + b 2x + b 3y + bi+x2 + b 5xy + b 6y 2. (6.41) Similarly for pressure p; p = ci + c 2x + c 3y + c^x 2 + C5xy + C g y 2 . (6.42) The partial derivatives u,^ and v, can also be written schematically as complete polynomials, \° u'x \ a i \ \ 1 \ 0 Constant a 2 \a 3 • \ \ x y \ 0 Linear 2a!,. a^ \ag 2\ \ ( M x z xy y ^0 Quadratic 3aj 2aQ ag \a-io x 3 x 2y xy 2 y3^\0 Cubic 4an3a^ 2 23I 33T_I\3 15 etc. etc. etc. 0/ v, b l 7 / 7 / 0/ 1 Constsnt b 2/ b 3 / / 0/ x y Linear b^/bc, 2bg / / ( O 0/ x 2 xy y 2 Quadratic by/bg 2bg 3bj_o 0/^x3 x 2y xy 2 y 3 Cubic b^/bi 2 2b^ 33bji+4bi5 etc. etc. etc. Consider only the non-zero terms to obtain u, +v, 3s J x y (u,x+v,y) 1 Constant (a 2+b 3) x y Linear (la^+b^) (a5+2bg) x 2 xy y 2 Quadratic (3a7+bg) (2ag+2bg) (sg+3bio) x 3 x 2y xy 2 y 3 Cubic (4a n+b 1 2) (3a 1 2+2b 1 3) (2a 1 3+3b 1 1 +) (a l l f+4b 1 5) etc. etc. etc. (d) I l l and p; 1 x y Pressure p Constant Linear c2 c3 x 2 xy y 2 Quadratic c^ C5 eg x 3 x 2y x y 2 y 3 Cubic C 7 CQ Cg C^Q 4 3 2 2 3 4 (e) x H x ay x z y z xy 3 y H Quartic C)i Cj2 c13 c14 c15 etc. etc. etc. The d i s c r e t i z e d continuity equations can now be obtained over the domain i n the following manner: J o e P ( u » v + v, )dfl = 0. 9c. ^ n e F V " » x ' V V 1 ^ (6.43) Let a i j = Joe xVdft e; 1=3=0,1,2,3, 'fte then from schemes (d) and (e) the following matrix form for the d i s c r e t i z e d continuity equations i s obtained; a00 2 u 1 0 a01 3 a 2 0 2 a u <*02. . . a o o a10 2 a 0 i a 2 0 2 a n 3ao2-<*10 2a 2o « 1 1 3 a30 2 a 2 1 a 1 2 . • -aiO a20 2 a n «30 2 a 2 1 3 a 1 2 . <*01 2 a n <*02 3 a 2 i 2a 12 °<03' . -a 0 1 a x l 2 a 0 2 a 2 1 2 a 1 2 3a 03-"20 2 a30 a 2 1 3aito 2 a 3 l a 2 2 . . . a 2 0 C30 2 a 2 1 °"t0 2 a 3 l 3a2 2. an 2 a 2 1 a 1 2 3 a 3 l 2 a 2 2 C13- . . a n a 2 i 2 a 1 2 «31 2 a 2 2 3 a 1 3 . «0 2 2a 12 a03 3 a 2 2 2«13 ag u. . . a o z a 1 2 2ao3 a 2 2 2^13 3a 0 t t. C 3 0 2a 4 0 a31 3 a 5 0 2aiti a32- • ^ 30 aifO 2 a 3 i «50 2ai+l 3 a 3 2 . a21 2a31 a 2 2 3 a 4 1 2 a 3 2 a23- . . a 2 i a31 2 a 2 2 a41 2 a 3 2 3 a 23-a 1 2 2 a 2 2 « 1 3 3 a 3 2 2 a 2 3 • . . a 1 2 »2 2 2ai3 <*32 2 a 23 3 a 1 4 . <*03 2«13 a04 3 a 2 3 2a ih C Q 5 - • -C03 «13 2aQit a 2 3 2 a m 3a 0 5-(Nx2Q) a4 a5 a7 as a 9 b 3 b 5 b 5 b 8 bg b10 = 0. (6.44) 112 Here N= n + 2C2, where n i s the degree of complete polynomial for p; Q=in+^C2> m being the degree of polynomials used f o r v e l o c i t i e s u and v; e.g. for t h u and v cubic, m=3 and Q=6. I t can be observed that every (Q+j) column th i s either equal to or a simple multiple of the j column. Therefore the rank of the matrix of c o e f f i c i e n t s i n (6.44) i s at most Q. Further the rank of the matrix i s s t i l l Q even i f N i s greater than Q. Thus for N greater than Q, a l l the d i s c r e t i z e d continuity equations of (6.44) are not indepen-dent. Now consider the case where u,v and p are complete quadratics, i . e . m=n=2 and Q=3, N=6. The r e s u l t i n g d i s c r e t i z e d i n c o m p r e s s i b i l i t y constraints are «00 2 a 1 0 a o i a00 <*10 2a o f a 2 <*10 2 a 2 0 a n a10 a20 2on ait a01 2a u <*02 a01 a n 2 a 0 2 as «20 2<*30 a 2 i «20 «30 2 a 2 i b 3 a n 2 a 2 i a 1 2 a n a 2 i 2 a 1 2 b 5 c<02 2 a 1 2 a03 «02 a 1 2 2«0 3 b 6 0. (6.45) C l e a r l y the rank of the matrix of c o e f f i c i e n t s i n (6.45) i s 3. Therefore three of the constraints i n (6.45) are not independent. The corresponding contribution from p(u»x+v,^)dfi to the equi-l i b r i u m equations i s a00 aio «01 a 2 0 a n ao2 c l 2a io 2a 2o 2 a n 2 a 3 0 2 a 2 i 2ai2 c2 a01 a n ao2 a 2 i a 1 2 ao3 c3 aoo a10 a01 a20 a n a02 aiO a20 a n a30 a21 a 1 2 c5 2 a 0 1 2a n 2 a 0 2 2 a 2 i 2 a 1 2 2a 03 c6 = ac. (6.46) 113 where the matrix of coefficients is simply the transpose of the matrix in (6.45). This is the equivalent of e f part of the equilibrium equations in (6.36). In general, ac cannot be equal to zero, i.e. pressure cannot be in equilibrium by i t s e l f . Therefore ac = 0 (6.47) should not have any non-zero solutions. This can be possible i f the rank of matrix a is the same as the degrees of freedom c^'s, i.e. 6 in the example considered. But the rank of a is obviously 3, thus leading to indeterminate c^'s which causes a self-equilibrating pressure f i e l d and non-unique solutions. However, this situation can be avoided i f the rank of the matrix a is equal to the number of constraints. This is possible i f the pressure distribution is taken as linear for quadratic distributions in u and v. The equations (6.45) and (6.46) then reduce to arjo 2a 1 0 a 0 1 a 0 0 a 1 0 -2a0i a 1 0 2 c*20 a l l a 10 a 20 2 a l l a 0 1 2 a u a 0 2 a 0 1 a n 2 a 0 2 a 2 ak A5 b 5 0, (6.48) and a 00 "10 a 01 2 a 1 0 2 a 2 0 2 a n a 01 a l l «02 a 00 a 1 0 a 0 1 a 10 <*20 a n 2 a o i 2aii 2ao2 ^2 C3 = ac, (6.49) 114 respectively. It can now be concluded that in general N should not be greater than Q in equation (6.44) to avoid self-equilibrating systems in pressure and is equivalent to saying that the degree of interpolating polynomial for pressure should not be larger than (m-1) where m is the degree of complete polynomials used for the velocity components u and v. The pressure on the boundary is incorporated as a natural boundary condition i n the functional of (6.33). However, inspection of equations (6.31) reveals that pressure needs to be fixed at some point in the domain fi as adatum. If the approximating polynomial in the f i n i t e element formula-tion f a i l s to comply with the requirement concluded above, then to avoid self-equilibrating systems, the pressure needs to be specified at more than one point on the boundary depending on the number of self-equilibrating modes present. At the same time the completeness requirement (i) for mean convergence of the stresses to velocities should not be overlooked. Therefore a consistent formulation of the element matrix would have only three rigid body modes, as in equations (6.34). Various combinations of interpolations for the velocities, pressure and stresses (u,v,p,x , T and T ) over a triangular and a rectangular ' xx' yy xy & & element (Figures 4 and 5, expect for addition of pressure degree of freedom at the nodes) are considered. The following eigenvalue problem is 0 a B " "V 0 0 " T a 0 0 p - A 0 I -p 0 P = 0 3 T 0 -b T 0 0 I - T T (6.50) (2m+J+3n)x(2m+I+3n) (2m+l+3n)x(2m+If3n) 115 where 6 = <u 1x2m 2mx3n T T V > , X lx3n T T T = <T X X > -xx -yy -xy a 0 b e , a = 0 b a 2mxl f and b 3nx3n i-h o d 0 The submatrices £jb,d,e,f, and the sub-column vectors are the same as in the matrix equation (6.36) whose derivation is similar to that of plane stress problem (Appendix A), where as I. , U and I are the identity matrices. 2mx2m 1%1 3nx3n It was mentioned earlier in this section that i f the rank of this matrix is (2m+I+3n), then i t has 2m positive and (I+3n) negative eigenvalues. If there are q ri g i d body modes and r mechanisms which correspond to zero eigenvalues, and since these correspond to indeterminacies of the. u and v degrees of free-dom (2m), then only (2m-q-r) eigenvalues are positive. Similarly, i f there are p indeterminacies of pressure degrees of freedom (1), then (I+3n-p) eigenvalues are negative. The distribution of negative, zero and positive eigenvalues and the composition of eigenvectors for different combinations of interpolations for u,v,p and the stresses x's are presented in Tables III and IV for a t r i -angular and a rectangular element, respectively. It can be observed that the combinations of interpolations which do not comply with the requirements on pressure and stresses resulted i n more than three eigenvalues required for rigid body modes. The self-equilibrating modes in pressure are obtained when the pressure interpolation polynomial i s of the same degree as the polynomials for u and v and with the exception of linear u,v,p,xxx,x^^ and x over a triangular element, have the same distribution as the interpolation xy polynomial while u,v and the x's are zero. The mechanisms are obtained when u,v are quadratic and x's constant over a triangular element; u,v biqua-dratic, x's bilinear over, a rectangular element, as.in the plane stress 116 problem of the previous section. In addition self-equilibrating pressure modes result when the pressure is quadratic for the triangular element. Again, when the rigid body modes are eliminated, the mechanisms seem to possess the same distribution for u and v as the assumed interpola-tions while pressure and the stresses are zero. In recognizing the rigi d body modes of equations (6.34), the pressure was said to be arbitrary. However, the eigenvectors for zero eigen-values associated with the rigid body modes and mechanisms displayed zero pressure in the element. This can be explained by spli t t i n g the eigenvalue problem of equation (6.50) in the following manner: - XI 6 + ag + §T = 0 (6.51a) aT6 - XI p = 0 (6.51b) -p-BT6 - [b + XI ] •= 0. (6.51c) Now the rigid body modes consist of u=a-cy and v=b+cy, therefore i t is T T clear that a 6 and B 6 are zero since these involve derivatives u, . v, and - - - - x y (u,^+v,x). The stresses are zero from (6.34), hence from equation (6.51a) ap_ = 0. (6.52) If the rank of a is equal to the number of constrians, i.e. the degrees of freedom in pressure, then the equation (6.52) is only true i f p=0. There-fore pressue i s zero everywhere i n the domain. Hence the equations (6.34) for rigid body modes are modified here (i) u = Constant, v = 0; p = x = x xx ( i i ) u = 0, v = constant; p = x =x =x =0 (6.53) R XX V " ( i i i ) u = cy, v = cx; p = Txx = T Thus a consistent formulation should not have more then three zero eigen-values for the rigid modes of equations (6.53). Finally static condensation = X = 0 yy xy = X = 0 yy xy = X = 0. yy xy 117 of the matrix equation (6.36) by eliminating the stresses, for the cases where no mechanisms are present, gives exactly the same matrix equation as one would obtain from the functional of (6.38) using the same u,v and p interpolation polynomials. In fact, the results of the eigenvalue analysis for such cases are very similar to those presented by Olson and Tuann [25]. In concluding this chapter i t should be pointed out that for certain combinations of approximate displacements and stresses (which indeed comply with the completeness requirement (i) and exhibit proper eigenvalues and eigenvectors over one element domain) the assembled element matrices according to certain continuity requirements yield self-equilibrating modes over the f u l l domain. The typical example i s that of plane linear elasticity with quadratic displacements, linear stresses over a triangular element and both continuous across the interelement boundaries. In this case, the boundary integrals on the common boundaries amongst adjacent elements cancel each other. Thus the element matrices for the elements, which do not have edges coinciding with the boundary of the problem domain, can be formulated without extracting the boundary integrals from the energy product. Then the contribution to the energy product from the integrals like /g U . T . . .dfi is zero for the degrees of freedom at the vertex nodes. fie 1 ^ »3 Thus zero rows and columns are obtained for displacement degrees of freedom at a l l internal vertices of the triangular elements when the element matrices are assembled. This then gives self-equilibrating modes over the f u l l domain. To check the existence of such modes, the element matrices for a certain for-mulation can be assembled so that there is at least one internal vertex node for the triangular elements and a corner node for the quadilateral elements, and then an eigenvalue-eigenvector analysis performed on the resulting matrix. 118 CHAPTER 7 APPLICATIONS OF THE MIXED FINITE ELEMENT METHOD The applications of the mixed f i n i t e element method to beam bending, plane linear elasticity and the problems with stress concentrations and singularities are presented in this chapter. The strain energy conver-gence rates for various formulations are predicted and compared with the numerical results obtained from the f i n i t e element analysis. The energy convergence of the mixed f i n i t e element method for plane e l a s t i c i t y with stress singularities, established in section 4.6, is also demonstrated with a numerical example. Finally the stress intensity factor K^ . for plates with symmetric edge cracks and a central crack are calculated and compared with the nearly exact values available. 7.1 Beam Problem Using the nomenclature of Figure 7, the following four f i r s t order f i e l d differential equations for simple beam theory result: - g - q - 0 (7.1) - f - V - 0 (7.2) de M dx EI 0 (7.3) £ - e = 0 (7.4) where (7.1) and (7.2) are the equilibrium equations (7.3) is the constitutive relationship (E=Young's Modulus, I=moment of inertia) and (7.4) is the con-straint equation arising from the assumption of plane sections remaining plane after deformation. If equations (7.2) and (7.4) are satisfied exactly, V and 0 in equations (7.1) and (7.3) can be eliminated and two second order aquations are obtained 119 d M n d 2v M _ 7 " FT = °-(7.5) (7.6) dx z EI The mixed method is applied to both systems, namely the four f i r s t order equations (7.1) to (7.4) and the two second order equations (7.5) and (7.6). 7.1.A Two Second Order Equations The equations (7.5) and (7.6) can be put into the matrix form as (7.7) or where D= "0 D2' V q M 0 AA = f (7.7a) dx* Here the matrix operator A is symmetric, i.e. (AA,A) = (A,AA), and the energy product is given by (AA, A) = / X 2 [M'-v+v'^-^rldx. (7.8) Integration of the f i r s t and the second term by parts yields (AA,A) = M'v|X2 + v'M|X2 - | X 2 [2M'v'+^]dx. The mixed variational principle can now be expressed as (7.9) XM = [2M'V4| r2qv]dx. (7.10) It can be observed from (7.10) that the continuity requirement has been reduced by one order compared to the Potential Energy approach with the variational form as I D = JI2 [EIv,,2-2qv]dx. (7.11) 120 This allows one to use lower order polynomials for approximating v and M. For a stationary point, the f i r s t variation of I in (7.10) i s zero 61 = -/A2[2M'6v,+2v,6M,+^M-2q6v]dx = 0. M . hi Again integration of the f i r s t and the second term on the right hand side yields 6L, = / X2 2(M"-q)6vdx+/X2 2 (v"-^r) 6Mdx-M' 6v I X 2 - v ' 6M jXz = 0. (7.12) M J x 1 ^ ; xj x EI 'xi 'x! This indicates that the forced boundary conditions are implied on the variables v and M which are different from v and v', one would find from the variational principle in equation (7.11) (ordinary Potential Energy Theorem). The two boundary terms in (7.9) could have also been obtained by twice integrating by parts the f i r s t term in (7.8) giving a variational principle of the form J M = [ 2 v " M - i f - 2 q v ] d x ( 7 , 1 3 ) which involves a second derivative v, thereby requiring the same continuity requirement on v as the potential energy approach. It can be shown that the forced boundary conditions for the mixed variational principle in (7.13) are implied on v and v'. Thus the mixed method offers f l e x i b i l i t y not just in incorporating the boundary conditions as mentioned in Chapter 5, but also in continuity requirements when dealing with higher order operators. The mixed variational principle in (7.10) should be distinguished from the one in (7.13) in that the latter follows from the energy product in symmetric form as defined in equation (3.73). As i t i s advantageous to have reduced continuity requirements in f i n i t e element analysis for problems involving higher order operators, e.g. (7.7), the boundary terms have to be extracted from both the equilibrium (7.5) and the constitutive (7.6) equations. The effect of such a formulation on the convergence i s considered. In the simple beam bending theory the stress, which is the bending 121 moment M, is proportional to the curvature, the second derivative of displacement, equation (7.6). In order to obtain improved convergence in the strain energy the approximations for M and v have to comply with the requirement for mean square convergence of M to v", so that the error in strain energy can be e s t i -mated from equation (3.71). This can be rewritten as [A-A0,A-A0] = (M-M0,M-M0) (7.14) for two second order beam equations. Here Mo=EIV" is the exact bending moment distribution, M, the f i n i t e element approximation, and V Q , the exact solution for the deflection v. Let the approximations for v and M be (7.15) M v = E $ .v. i - l 1 1 N M = E Y.M., J - l 2 3 where $. and Y. e H^=VxM, the cross product space ($. e V, Y. e M). i J A l j The substitution of (7.15) into (3.68) gives p . _ , - a . , + , § _ « . , and minimization with respect to M,. yields N M E (Y.,Y.)M. = Z (Y.,$'.')v., i=l,2, . . ., N. (7.16) j = l 1 J J j = 1 I J J However, the mixed variational principle in (7.10) gives the following equation at extremum instead of (7.16), i.e. N N E . (Y.,Y.)M. = - E (Y!,$!)v., i=l,2, . . ., N. (7.17) j = l 1 J J j=l 1 J J The right sides of (7.16) and (7.17) imply that (4\ > v") = _ ( r ,v'), (7.17a) which means that the derivative of v', which is v", is taken as a generalized derivative. This is because v' is only required to be piecewise continuous and therefore would not possess an ordinary derivative everywhere. 122 As was done in Chapter 6, inspection of the functional I in (7.10), for g=0, reveals that the element matrix should have the following r i g i d body mode v = constant, M = 0. After trying out different polynomials for v and M, i t i s found that the inter-polations should be of the same degree in order to achieve the only r i g i d body mode as mentioned above. Such approximations also meet the requirement (i) of completeness. If the boundary conditions are homogeneous equations (7.16) and (7.17) are equivalent and the error in the energy product can be estimated from (7.14). A beam of uniform cross section with constant stiffness EI, length 1 is subjected to constant load per unit length q. The following four cases of boundary conditions are considered, Figure 8: (1) simply supported beam (S.S); v(0)=M(0)=v(J)=M(l)=0. (2) cantilever; v(0)=M(l)=0. (3) both ends clamped (fixed-fixed); v(0)=v(J)=0. (4) one end clamped and the other in a vert i c a l guide (fixed-guided); v(0)=0. (Tangent to the elastic curve at vert i c a l guide remains horizontal). Three different combinations for approximating v and M within the element are chosen (i) v-linear, M-linear; ( i i ) v-quadratic, M-quadratic; ( i i i ) v-cubic, M-cubic. The nodes per element and the number of degrees of freedom per node for the combinations ( i ) , ( i i ) and ( i i i ) are shown in Figure 9. The derivation of the element matrices i n a l l three cases is analogous to that of plane stress problem presented in Appendix A. In each case the beam is divided into 123 elements of equal length, 1^. The following quantities, where necessary, are tabulated for presentation purposes: 6 = deflection (v); 0 = rotation (v*); M = moment; V = shear (M') ; U = strain energy. The subscript M stands for middle, Q for quater point, E for end, RE for right hand end, and LE for l e f t hand end. (i) v-linear, M-linear The results are shown in Table V and the plots of quantities of interest versus the number of elements to show convergence appear in Figures 10. Linear approximations for v and M provide the continuity of v and M across the nodes and do not violate the completeness requirement ( i ) . Since the second derivative of linear approximation for v vanishes, equation (7.16) is satisfied only in a generalized sense, i.e. equation (7.17). For linear approximations within the element, the error in both v and M can be shown to be 0(1 2) and Q(lk) in I I M I I 2 , where I is the element length ( I =4-). From Figure 10(a), the relative error in strain energy converges as N~2 for the simply supported and cantilever configurations (cases 1 and 2, respectively), from below in case 1 (Table V(a)) and from above in case 2 (Table V(b)); and N - 4 from below for the fixed-fixed and fixed-guided ones (cases 3 and A, respectively) (Tables V(c) and V(d)). The expected rate of strain energy con-vergence is obtained when the moments are not forced to be zero (not as the forced homogeneous boundary condition on stress). Perhaps, for the cases when the bending moment is forced to be zero on the boundary, some lower 124 order error terms prevail. Figure 10(b) i l l u s t r a t e s the monotonic convergence of mid-point d e f l e c t i o n 6^ f o r the simply supported beam from below and t i p d e f l e c t i o n 6„ for the c a n t i l e v e r from above, whereas the r e l a t i v e error con-E verges as N - 2 , Figure 10(c). The bending moments at the nodes are exact for both cases 1 and 2. However, the reverse i s true for the other two cases, as can be observed from Tables V(c) and V(d). The computed de f l e c t i o n s at the nodes are exact, whereas the r e l a t i v e error i n the fixed moment con-verges as N - 2 , as i l l u s t r a t e d i n Figure 10(d). The other q u a n t i t i e s , appearing i n Tables V, seem to converge with increasing number of elements N f o r a l l four cases. ( i i ) v-quadratic, M-quadratic With v and M both quadratic, the completeness requirements f o r the energy convergence are s a t i s f i e d since the va r i a b l e s v and M are continuous across the nodes and have piecewise continuous f i r s t and second d e r i v a t i v e s . The errors i n v and M can be shown to be 0(1 3) which leads to an error of e O(l^) i n ||M|| 2. However, a quadratic approximation for M i s capable of representing the exact s o l u t i o n MQ for the constant load q along the beam length. Therefore the s t r a i n energy from the f i n i t e element s o l u t i o n i s expected to be exact. This i s confirmed by the r e s u l t s i n Tables VI for a l l four cases. Also the moments and the derived shears obtained are exact. I t i s i n t e r e s t i n g to note that the mid-deflections 6^ , i n a l l four cases, f o r an even number of elements obtained are exact and the r e l a t i v e error f o r an odd number of elements along the beam length converges as N _ t |, Figure 11; while the end computed d e f l e c t i o n 6 f or cases 2 and 4 are exact for odd RE and even number of elements. When a fa s t e r convergence i s observed for the displacement, the r e l a t i v e error i n the derived r o t a t i o n , i n a l l four cases appears to converge as N - 2 , Figure 11. Tables VI also i n d i c a t e convergence 125 of other nodal variables. ( i i i ) v-cubic, M-cubic With cubic approximations, both M and v and their f i r s t deriva-tives are continuous. Here again, the completeness requirements are satis-fied and like the previous approximations ( i i ) the strain energy is expected to be exact as well as the moments and the shears. The numerical results presented in Tables VII confirm this. The end deflections in the cases 2 and 4 are also exact. However, the relative error in the mid-deflection 6^ converges as N - ! + and in rotations as N~3 for a l l four cases; Figure 1 2 . 7.1.B Four First Order Equations The four f i r s t order beam equations (7.1) to (7.4) can be put into the matrix form as 0 D -EI V q 0 0 M 0 V 0 or AA = f where the matrix operator A is symmetric, i.e. (AA,A)=(A,AA). The energy product is given by / A » ,\ f r dV dM„ „ „, de,, M2,dvT7, , < A A . A ) = j1 [ - ^ - d ^ - 2 v e ^ - - + d - V ] d x . (7.18) (7.18a) (7.19) Integrating by parts the f i r s t and the fourth terms on the right hand side gives ,dv „dM„ (AA ,A) = -Vv|X2 + 9M|X2 + f [ 2 V ^ - 2 ^ 6 - 2 v e - ^ ] d x 'xj 'xi J l dx dx EI = -Vv| X 2 + 6M|X2 + [A,A] . x l (7.20) 126 Here [A,A] is the modified energy product of equation (3.73); The equation (7.20a) then leads to the mixed variational principle for the equations (7.18) with forced boundary conditions on v and M. The mixed Galerkin Method (by adding boundary residuals to the residuals of equations (7.1) and (7.3), section 5.4) and the mixed variational principle (7.21) would yield the same element matrix since the problem is linear and self-adjoint. Appendix C shows that the linear elasticity equation (5.9) yields the same energy product (7.19) as equations (7.18) when the basic assumptions of the simple beam bending theory are incorporated and the shear strain energy term, involving V 2, is neglected. Therefore the error in the energy product can be predicted from equation (3.71) provided approximations for the displace-ments (v,6) and stresses (M,V) are complete. Linear approximations are chosen for v,9,M and V. The element nodes and the nodal degrees of freedom are the same as illustrated in Figure 9, combination ( i i i ) . Further, the approximations chosen are complete and the error i n the energy product is governed by the mean square error in M, i.e. error in ||M|| 2, which is 0(1^) for the linear M distribution. Again four cases of boundary conditions (Figure 8) are considered. The results obtained from the f i n i t e element analysis are tabulated in Table VIII and convergence plots are shown in Figures 13. It is observed from Figure 13(a) that the relative error in strain energy converges as N _ l + for a l l four cases. Moments and shears at the nodes in cases 1 and 2 for both even and odd number of elements N are exact. However for cases 3 and 4, shears 127 are exact for both even and odd N but only end moments are exact for N even and they converge as N-1* for N odd in case 3. From Figure 13(b), the rela-tive error in rotation 0 converges as N - 2 in a l l cases. Deflections 6 at the free end for N even in case 2 and the guided end for N odd and even in case 4 are exact, converging as N - i +for odd N in case 2; whereas the mid-deflections 6 appear to converge to the exact value in an oscillatory manner. It has been mentioned in Chapters 3 and 5 and section 7.1.A that the mixed methods allow f l e x i b i l i t y in incorporating the boundary conditions. If the f i r s t and the second terms on the right hand side of equation (7.19) are integrated by parts, the following mixed variational principle results: X M = h [2V^+2Md|-2ve-g-2qv]dx (7.22) with forced boundary conditions on v and 0, the same as for the Potential Energy Theorem. Again linear v,0,M and V are used. Since the completeness requirements are not altered by shifting the forced boundary condition from M to 0, the convergence of strain energy is s t i l l expected to be N-t+. The strain energy for cases 1, 2 and 3 is computed and tabulated in Table IX. The relative error in strain energy versus N is then plotted in Figure 14 and in a l l cases the strain energy i s found to converge as N-t+. Next the shear strain energy term V 2 is also included and the mixed variational principle of (7.22) now includes an additional term (Appendix C) V - / I ^ a ^ - w o - g - ^ v ^ ^ ( 7 . 2 3 ) where v is poisson's ratio and h the height of the beam. Again the forced boundary conditions are on the variables v and 0. Using the same linear 128 approximations for v,9,M and V within an element, which s t i l l complies with the completeness requirements, cases 1, 2 and 3 are analysed and the results appear in Tables X. In the analysis, v is taken as 0.25 and h as ^ - The relative error i n strain energy converges as N _ t + i n a l l three cases as shown Figure 15(a). The shear V, in cases 1 and 3, is exact and converges as N~2 for the cantilever, case 2; the relative errors in the mid-moment for case 1, the fixed moment for case 2 and the end moments for case 3 converge as N"-2; the mid-deflection <5W for cases 1 and 3 converges as N - 2 while the free end deflection 6 of the cantilever converges as N~^ as shown in Figures 15. However, the end rotation 8 for the cases 1 and 2 is exact for N even (Tables X) and from Figure 15(b) i t appears to converge as N-** for N odd in case 2. Despite the several different convergence rates observed for basic variables v,9,M and V for the three different formulations considered above, the relative error in the strain energy in a l l cases converges as N-1+. Further, where the boundary integrals were taken out from the equilibrium equations, i.e. forced boundary conditions on v and 0, the strain energy converges from below. Alternately i t is from above in the case where the forced boundary conditions are on the displacement v and the moment M. 7.2 Plane Linear Elas t i c i t y The stress and the displacement are chosen to be linear within a three node triangular element with five degrees of freedom per node (Figure 4) and are forced to be continuous across the interelement boundaries by equating the nodal variables at common nodes. These approximations satisfy the completeness requirements. 129 The derivation of the element matrix and the consistent load vector arising either from body forces or nonhomogeneous stress boundary conditions are given in Appendix A. The f i n i t e element thus formulated is then applied to solve the following problems. 7.2.A Plane Stress: Square Plate with Parabolically Varying End Loads A square plate with parabolically varying end loads is shown in Figure 16. Since the problem is symmetric about the x and y axes, only a quarter of the plate ABCD is considered in the f i n i t e element analysis with the forced boundary conditions u=0 on AD and v=0 on AB. Since the displacements u, v and the stresses x , x and x xx yy xy are assumed linear, the error in the stresses is 0(2 2), where 1 is the e e largest diameter within the element. Further, since such approximations satisfy the completnness requirements, and displacements and stresses are continuous across the interelement boundaries, equation (3.71) holds and error in the strain energy is expected to be the mean square error in the stresses, i.e. 0(1^), which is 0(N_1+) for a uniform grid, Figure 16. The numerical results for some of the stresses and displacements at points A,B,C, and D and the strain energy from the mixed f i n i t e element analysis for various grids are presented in Table XI. Also presented in Table XI are the results from the displacement element, obtained by Cowper, Lindberg and Olson [5], (using f u l l cubics for u and v displacements over a triangular element with six degrees of freedom (u,u ,u ,v,v ,v ) at x y x y vertices and two (u,v) at the centroid) for comparison. Since no attempt is made to satisfy the stress boundary conditions, the correct values of stresses on the boundary are obtained only in the limit of N-^. The convergence plots for the mixed element are shown in Figures 17 130 and those for the displacement element in Figure 18, [5J. The energy con-vergence rate for the former appears to be very close to N - 1 + (Figure 17(a)) as predicted. Although the error in strain energy from the displacement element is much smaller than from the mixed, owing to the fact that the former is a much more refined element, the energy convergence rate is slightly less than N-^, Figure 18, which is lower than the predicted assymp-totic rate of N - 6. The other interesting observation about the energy convergence for the mixed f i n i t e element is the convergence from below when the forced boundary conditions are on the displacement variables which has also been observed in section 7.I.B. The convergence rate for the stresses from the mixed element, Figure 17(c), appears to be close to N~2 for N larger than 8. However, peculiar kinks are observed and can be associated with the fact that certain stresses were fortuitously close to their exact values for a certain grid; e.g. N g f° r N=6, etc. Figure 17(b) shows the convergence of displacements indicating faster convergence for u^ and (close to N-l+) than u^ and v^ (close to N — 2). This is also observed for the displacement element and N greater than 6, Figure 18. In the mixed element, the kinks are observed in the convergence plots for the displacement u,, at N=6 and v at N=6 and N=10, a C which can be associated with slightly larger errors for such grids. 7.2.B Cantilever (Plane Stress) The dimensions, loading and the material properties are detailed in Figure 19(a). Two types of boundary conditions are considered at the fixed end as indicated in Figures 19(b) and 19(c); B.C.I and B.C.2. The latter is used for comparison purposes since the solutions using various f i n i t e elements for B.C.2 are available in the literature, while the former is considered to investigate the energy convergence. 131 (a) Cantilever with Boundary Conditions B.C.I In an effort to investigate the energy convergence, the exact value of strain energy is necessary. Besides such boundary conditions being easily incorporated in the f i n i t e element analysis, i t is also possible to obtain an elasticity plane stress solution under the assumption of bilinear normal stress in the x-direction and the shear stress independent of x and quadra-t i c in y. The solutions for the displacements u and v along with the strain energy for the boundary conditions B.C.I are presented in Appendix D. Since the part of the boundary between A and F, and E and F can move in either direction and i f the stresses on the l e f t face were not in-cluded in the f i n i t e element analysis through a consistent load vector, a stress-free boundary w i l l be simulated. This violates the assumptions of the elasticity solution which shall put in doubt the validity of the exact strain energy to be used in the error analysis. Therefore the stresses on the l e f t end are included in the consistent load vector. The typical mesh used in the analysis is shown in Figure 22 and the numerical results are presented in Table XII. The results indicate that the stress T at x=12 inches and y=-6 inches, and the tip deflection xx are converging to the exact values in an oscillatory manner. However, the strain energy is converging to the correct value from below. The convergence plots are shown in Figures 20. In Figure 20(a), the plot of relative error in strain energy versus N the number of elements in the beam depth, the inclusion of N=6 leads to a kink in the plot. Note the grid for N=6 does not contain the previous grids for N=2 and N=4, and excluding this the energy appears to converge as N - L f as predicted. The same behaviour is also observed in Figure 20(c) for the relative error in the stress T versus N for N>4 ° xx (close to N - 2 without N=6). This is not surprising since the error in strain energy is governed by the mean square error in the stresses. It is gratifying 132 that the energy convergence rate is about double that for stress as predicted. Figure 20(b) does not indicate any definite convergence rate for the tip deflection 6 and perhaps, more data is required to establish any trend of convergence. (b) Cantilever with Boundary Conditions B.C.2 Unfortunately the exact solution for the boundary conditions B.C.2 is not available and hence the exact strain energy is not known. The problem is solved to compare with the solutions obtained by using different displacement f i n i t e elements for various grids. These are readily available in literature, e.g. Gallagher [9]. Here the boundary AE (Figures 19) is pre-vented from moving in either direction and leads to stress singularities at corners A and E. Furthermore the shear stress and normal stress distributions at the fixed boundary are not the same as assumed in the beam theory. How-ever the results are compared with the beam theory [35] which provides an upper bound for the tip deflection from the displacement f i n i t e element. The numerical results from the mixed f i n i t e element analysis for various grids (Figure 22) are tabulated in Table XIII. Again the stress T at x=12 inches and y=-6 inches, and the tip deflection appear to converge X X in an oscillatory manner. The strain energy for N=8 is slightly higher than the strain energy 1/2P6 obtained from the beam theory whereas in the previous examples, when the boundary integrals were extracted from the equilibrium equations, the energy converged from below. Since the exact numerical value is not known, the convergence of strain energy is rather d i f f i c u l t to establish. Table XIV shows the comparison of numerical results from the mixed f i n i t e element with those from the displacement models, e.g. constant stress triangle (C.S.T.), linear stress triangle (L.S.T.) and quadratic stress 133 triangle (Q.S.T.) elements, for the various grids shown in Figures 21. The mixed f i n i t e element used here appears to perform slightly poorer than the L.S.T., which uses quadratic approximations for the displacements and much better than the C.S.T., using linear displacements. The Figures 23(a) and 23(c) also indicate fast convergence of the tip deflection with increas-ing degrees of freedom and N, the number of elements in the beam depth, respectively, relative to the other elements. The graph of strain energy versus N shown in Figure 23(b) also shows rapid convergence. Finally the relative error in tip deflection is plotted against N for the mixed f i n i t e element and other displacement f i n i t e elements in Figure 24. It appears as in case (a) that more data is required to establish the convergence rate. However the plot .does exhibit fast convergence. It should be noted that the tip deflection from beam theory is used as exact solution in plotting these curves, and i t is in error i t s e l f . 7.2.C Stress Concentration around a Circular Hole (Plane Strain) A square plate (plane strain) with a circular hole in the middle (Figures 25) loaded by a uniform uniaxial stress to is considered. The diameter of the hole is one eighth of the plate width and the plate is of unit thickness. The plane strain state is analysed for both isotropic and orthotropic cases. It is demonstrated in Appendix A, how the element matrix for a plane stress isotropic case is modified for plane strain isotropic and orthotropic cases. The procedure is much simpler than for a displacement fi n i t e element. Because of symmetry only a quarter of the problem is con-sidered. The grid and the boundary conditions used in the fi n i t e element analysis are shown in Figure 26. This is essentially the same as used by Zienkiewiz, Cheung and Stagg [42] for constant stress triangular elements. 134 A comparison between the analytical solutions (for the isotropic case, Timoshenko and Goodier [35], and the orthotropic case, Savin [31]) for an in f i n i t e plate with a circular hole in the middle ( T along edge BC and T along edge AE) and the solutions obtained from the mixed f i n i t e element yy analysis is shown in Figure 27. A similar comparison with the solutions from the constant stress triangles [42] are shown in Figure 28. The mixed f i n i t e element solution shows excellent agreement with the analytical solution, and further the stresses are obtained directly at the nodes. The constant stress triangles also show good agreement with the exact solution, but the stresses are computed by averaging at nodes from the neighbouring elements, assuming the constant stress within the element to be the stress level at the node. Further the concentrations occuring at the boundaries are obtained by extrapo-lation. 7.3 Stress Singularities The strain energy convergence for plane stress elasticity with stress singularities, established in section 4.6.A, is demonstrated by a numerical example. The stress intensity factor K^. i s then determined from the method described in section 4.6.B for rectangular plates with symmetric edge cracks and a central crack (mode type I, Figure 1). 7.3.A Strain Energy Convergence The problem of a square plate with symmetric edge cracks (mode type I) is considered. Figure 29(a) shows the problem description and Figure 29(b) illustrates the f i n i t e element idealization of the quarter of the plate con-sidered because of the symmetry about the x and y axes. The problem is solved using mixed f i n i t e elements for various grid sizes for two cases. The stress x is kept continuous across point D (the crack tip) in the f i r s t case and 135 in the second case, an extra node is introduced along the x-axis next to the original one at D and only U , V , T and x degrees of freedom are equated at xx xy the two nodes, thus allowing x to be discontinuous across point D, the crack yy tip. (See distributions in Figure 32). The numerical results for both cases are presented in Tables XV. The strain energy is converging from below in both cases, while the peak stress T v v D a t the crack tip is about 28 percent higher when the normal stress i s discontinuous across the point D, than for the case when i t is continuous. The plots of strain energy versus the mesh size appear in Figures 30. The shape of the curves in both cases are very similar and they exhibit faster convergence than just linear as might have been expected. Figure 30(a) shows a comparison with the solutions obtained using various other elements. The present mixed element definitely shows a faster strain energy convergence than the constant stress triangles, the linear stress triangles, and the hybrid stress rectangles with cubic stress distribution within the element and quadratic displacements along the boundaries. The convergence rate is indicated by the plot of the relative error in strain x 2L 2 energy (exact U=3.228—, Tong and Pian [38]) versus N, the number of elements C i t along the edge OA, Figure 31. It can be observed that the convergence rate approaches N - 2 as N gets larger, for both cases. It is clearly faster than N - 1 indicating the cancellation of the errors in the energy product of equation (4.84) due to stress singular terms. Further, a slightly larger error is observed in the case of discontinuous normal stress at the crack tip. Finally the normal stress x is plotted along the edge OA in Figures 32. In both cases, a small zone of compression is observed on the stress free edge of the crack with a peak value of about XQ (the applied stress on edge BC) close to the crack tip. 136 7.3.B Evaluation of the Stress Intensity Factor K^ . Two plane s t r a i n problems, rectangular plates one with symmetric edge cracks and the other with a ce n t r a l crack are considered. The d e t a i l s for these are shown i n Figures 33 and the f i n i t e element i d e a l i z a t i o n i n Figure 34. The layout of the mesh, used i n both problems, i s analogous to the one used by Parks [27] with the exception that the present elements are tria n g u l a r . Also indicated i n Figure 34 are the ra t i o s of the r a d i i y„ r0 to the crack length a (0, 0.1, 0.2 and 0.5). These are then used to calcu-l a t e the p o t e n t i a l energy release rate for a crack extension of Aa i n the f i n i t e element analysis as described i n section 4.6.B. Because of symmetry about the x and y axes, only a quarter of the problem i s considered i n each case and th i s i s shown i n Figures 33 as shaded areas along with the respective boundary conditions. Although, i t i s s u f f i c i e n t to solve each problem once for the i n i t i a l crack length a (section 4.6.B), at present the f i n i t e element analy-s i s i s performed every time when the countour TQ i s translated along with the i n t e r i o r nodes by the amount Aa=5xl0 6 a i n the x- d i r e c t i o n . The p o t e n t i a l energy release rate, i n the d i s c r e t i z e d form can be expressed as A 7 IM G I = -AT <7-55> and An., = TL, - TL, M Mr Mg 1 0 where TL, i s the p o t e n t i a l energy associated with the i n i t i a l crack and TL, MQ Mr 1 0 when the crack t i p has been m o v e d by the amount Aa. Then the crack i n t e n s i t y factor i s calculated by 137 The numerical results for the plate with symmetric edge cracks are lis t e d in Table XVI and those for the plate with a central crack in Table VII. The crack intensity factor for the former is compared with the nearly exact K^. obtained by Bowie [2] and for the latter, with K by Bowie and Neal [3]. In both cases the results obtained are in excellent agreement with the references. It is seen that the least percentage error is obtained for the contour Tn with radius Y^ =0.la, and the worst for y p =0, i.e. only R 0 R 0 the crack tip node is translated. The former is also associated with the highest potential energy release rate which varies for different sets of nodes defining TQ. When the calculated value of G^. or J-integral is in fact indepen-dent of the particular set of nodes defining TQ, the mesh may be called optimal. Thus, i t suffices, for optimal meshes, to move only the exterior node defining the crack tip, thereby altering the boundary, regardless of the particular set of interior nodes comprisng the contour TQ. Alternatively, non-optimal meshes w i l l exhibit some path dependence in the calculated values of G^.. In such cases, personal judgement and experience can help determine the best value of G^.. In Table XVIII, a comparison with stress intensity factors obtained from the energy release rate by other authors is presented. It can be seen that excellent accuracy is obtained with much fewer mixed f i n i t e elements and degrees of freedom than the corresponsing displacement models. Finally the plots of normal stress on the cracked face OA (Figures 33) are shown in Figures 35. Note that a higher peak stress is obtained at the crack tip than the peak stress Indicated in Figure 32(a), probably because of the refined mesh near the crack tip. 138 CHAPTER 8 CONCLUSIONS A detailed investigation of the theoretical foundation and practical aspects of applying mixed methods of approximate analysis for continuum and fi n i t e element analysis has been presented. Mixed methods always involve Indefinite operators and consequently a new energy product and i t s associated energy norm had to be introduced for these special operators. The fields of definition of such operators were so restricted that when obeyed, the energy product was found to be positive definite and represented twice the strain energy. These new concepts then formed the bases for establishing the energy convergence of complete approximations for displacements and stresses. It was found that the energy convergence implied the mean square convergence of such approximations to the exact values. The completeness requirements for continuum analysis were defined in two steps: (i) the mean square convergence of the strains from the approxi-mate stresses to the strains derived from the approximate displacements; (i i ) convergence of the energy norm. The alternate form of the requirement (i) is as follows: The strains from the stress approximations should possess at least all the strain modes that are present in the strains derived from the displacement approximations. It was also concluded that a violation of this requirement leads to mechanisms and this was confirmed by the eigenvalue-eigenvector analysis of an element matrix. The presence of mechanisms also indicates the breakdown of positive definiteness of the energy product, hence re-quirement (i) is the prerequisite of ( i i ) . The error in the energy product was shown to be proportional to the mean square error in the stress approximation when completeness is satisfied. This leads to much faster convergence in the 139 strain energy calculated from the mixed method than obtainable from the corresponding displacement method, i.e. the latter with identical displace-ment approximations as used in the mixed method. The foregoing concepts for the continuum were then extended to the f i n i t e element method and the corresponding completeness c r i t e r i a were established. These were found to be as follows: a) the displacement approximations should include a l l rig i d body and con-stant strain modes; b) the stress approximations should include a l l constant stress modes; c) the same as requirement (i) for the continuum, i.e. a l l the displacement strain modes should be included in the stress (strain) approximations. It was also concluded that for complete approximations the strain energy convergence cannot be any faster than for that calculated from the corres-ponding displacement model, unless the stresses are made continuous across the interelement boundaries. In the example of beam bending with four f i r s t order equations the use of linear interpolations for the four basic variables resulted in a predicted mean square error in stresses of 0(N _ l f). Therefore the predicted error in strain energy was also of 0(N_ 1 +) and this was confirmed by the numerical examples. The plane stress triangular element using linear interpolations for both displacements and stresses yielded a predicted error in stresses of 0(N - Z) , and a mean square error and strain energy error of 0(N - t +). In the numerical applications of this element, the energy convergence rate was indeed found to be 0(N - 1 +) for the plane stress square plate with parabolic end loads and nearly the same for the plane stress cantilever. In comparison the corresponding displacement element (C.S.T.) yields a convergence rate of only 0(N - 2). A faster energy convergence rate (nearly 0(N - 2)) was also 140 observed for the plane stress square plate with symmetric edge cracks for which the strain energy converges only linearly with N, even with higher order displacement or hybrid-type f i n i t e elements, Tong and Pian [38]. Excellent accuracy was also obtained for the stresses around a circular hole in the middle of a square plate subjected to uniaxial com-pression for plane strain isotropic and orthotropic cases. Finally the crack intensity factors (K^) computed for plane strain rectangular plates, one with symmetric edge cracks and the other with a central crack, yielded errors of only 1.97% and 0.89%, respectively. The matrix equations to be solved in the mixed f i n i t e element analysis are always indefinite and have zeroes on the diagonals for the displacement degrees of freedom. The method of Gaussian elimination with partial pivoting was successfully employed to solve such equations. Hence i t is concluded that the indefinite nature of the mixed method equations presents no special d i f f i c u l t i e s . In general, methods involving indefinite operators preclude obtaining upper or lower bounds on energy. In the applications of the mixed fi n i t e element method discussed herein, the energy was observed to converge in some cases from above and in some from below. However, in the examples where the variational principle was formed by extracting the boundary integrals from the equilibrium equations, the strain energy always converged from below. In the examples solved, far more accurate results were obtained by using the mixed f i n i t e element method than the corresponding displacement method with the same displacement approximations. However, the mixed method required more degrees of freedom for the same number of elements. On the other hand, the results for the stress concentration and stress singular problems were generally more accurate even using fewer elements (and total 141 number of degrees of freedom) than for the displacement models. Hence i t seems f a i r to conclude that the mixed f i n i t e element method can produce more efficient solutions for problems involving stress concentrations or singularities. 142 BIBLIOGRAPHY 1. Anderson, G.P., Ruggles, V.L., and Stibor, G., "Use of Finite Element Computer Programs in Fracture Mechanics", International Journal of Fracture Mechanics, Vol. 7, No. 1, 1971, pp. 63-76. 2. Bowie, O.L., "Rectangular Tensile Sheet with Symmetric Edge Cracks", Journal of Applied Mechanics, 31, 1964, pp. 208-212. 3. Bowie, O.L., and Neal, D.M., "A Note on the Central Crack in a Uniformly Stressed Strip", Engineering Fracture Mechanics, Vol. 2, 1970, pp. 181-182. 4. Cook, R.D., Concepts and Applications of Finite Element Analysis, John Wiley & Sons, Inc., New York, 1974. 5. Cowper, G.R., Lindberg, G.M., and Olson, M.D., "A Shallow Shell Finite Element of Triangular Shape", International Journal of Solids Structures, Vol. 4, 1970, pp. 1133-1156. 6. Dunham, R.S., and Pister, K.S., "A Finite Element Application of the Hellinger-Reissner Variational Theorem", Proceedings of the Second Conference on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, AFFDL-TR68-150, 1968, pp. 471-487. 7. Finlayson, B.A., The Method of Weighted Residuals and Variational Principles, Academic Press, New York, 1972. 8. Finlayson, B.A., and Scriven, L.E., "The Method of Weighted Residuals— A Review", Applied Mechanics Review, 19, 1966, pp. 735-748. 9. Gallagher, R.H., Finite Element Analysis, Prentice-Hall, Englewood C l i f f s , New Jersey, 1975. 10. Hellan, K., "Analysis of Elastic Plates in Flexure by a Simplified Finite Element Method", Acta Polytechnica Scandinavica, C i v i l Engineering Series No. 46, Trondheim, 1967. 11. Hellwig, G., Differential Operators of Mathematical Physics, Addison-Wesley Publishing Company, London, (transl.) 1967. 12. Herrmann, L.R., "A Bending Analysis of Plates", Proceedings of the First Conference on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, AFFDL-TR66-80, 1966, pp. 577-604. L3. Herrmann, L.R., "Elasticity Equations for Incompressible and Nearly Incompressible Materials by a Variational Theorem", American Institute of Aeronautics and Astronautics Journal, Vol. 3, No. 10, 1965, pp. 1896-1900. L4. Huebner, K.H., The Finite Element Method for Engineers, John Wiley & Sons, New York, 1975. L5. Hutton, S.G., "Finite Element Method—A Galerkin Approach", Ph.D. Thesis, University of British Columbia, Canada, September 1971. 143 16. Hutton, S.G., and Anderson, D.L., "Finite Element Method: A Galerkin Approach", American Society of Civil Engineers, Engineering Mechanics Division, October, 1971, pp. 1503-1520. 17. Lorch, E.R., Spectral Theory, University Texts in the Mathematical Sciences, Oxford University Press, 1962. 18. Mikhlin, S.G., Variational Methods in Mathematical Physics, The Macmillan Company, New York, 1964. 19. Oden, J.T., Finite Elements of Nonlinear Continua, McGraw-Hill, New York, 1972. 20. Oden, J.T., "Generalized Conjugate Functions for Mixed Finite Element Approximations of Boundary-Value Problems", The Mathematical Foundations of the Finite Element Method—with Applications to Partial Differential Eguations, A.K.. Aziz edition, Academic Press, New York, 1972, pp. 629-669. 21. Oden, J.T. , and Reddy, J.N., "On Dual Complementary Variational Principles in Mathematical Physics", International Journal of Engineering Science, Vol. 12, 1974, pp. 1-29. 22. Oden, J.T., "A General Theory of Finite Elements I Topological Con-siderations", International Journal for Numerical Methods in Engineering, Vol. 1, No. 3, 1969, pp. 205-221. 23. Oden, J.T., "A General Theory of Finite Elements II Applications", International Journal for Numerical Methods in Engineering, Vol. 1, No. 3, 1969, pp. 247-260. 24. Oliveira, E.R.A., "Theoretical Foundations of the Finite Element Method", International Journal of Solids Structures, Vol. 4, 1968, pp. 929-952. 25. Oliveira, E.R.A., "Completeness and Convergence in the Finite Element Method", Proceedings of the Second Conference on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, AFFDL-TR68-150, 1968, pp. 1061-1090. 26. Olson, M.D., and Tuann, S.Y., "Primitive Variables versus Stream Function Finite Element Solutions for the Navier-Stoke 1s Equations", International Centre for Computer Aided Design—Proceedings of the Second International Symposium on Finite Element Methods in Flow Problems, S. Margherita Ligure (Italy), June 14-18, 1976. Also see: Tuann, S.Y., and Olson, M.D., "A Study of Various Finite Element Methods for the Navier Stoke's Equations", Structural Research Series, Report No. 14, Department of C i v i l Engineering, University of British Columbia, Vancouver, Canada, May 1976. 144 27. Parks, D.M., "A Stiffness Derivative Finite Element Technique for Determination of Crack Tip Stress Intensity Factor", International Journal of Fracture, Vol. 10, No. 4, December 1974, pp. 487-501. 28. Pian, T.H.H., and Tong, P., "Basis of Finite Element Methods for Solid Continua", International Journal for Numerical Methods in Engineering, Vol. 1, 1969, pp. 3-28. 29. Rice, J.R., "A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks", Journal of Applied Mechanics, Vol. 35, 1968, pp. 379-386. 30. Reddy, J.N., and Oden, J.T. , "Convergence of Mixed Finite Element Approximations of a Class of Linear Boundary-Value Problems", Journal of Structural Mechanics, Vol. 2, 1973, pp. 83-108. 31. Savin, G.N., Stress Concentrations around Holes, Pergamon Press,'1961. 32. Strang, G., Linear Algebra and its Applications, Academic Press, New York, 1976. 33. Strang, G., and Fix, G.F., An Analysis of the Finite Element Method, Prentice Hall, Englewood C l i f f s , New Jersey, 1973. 34. Taylor, C., and Hood, P., "A Numerical Solution of the Navier Stoke's Equations using the Finite Element Technique", Computer and Fluids, Vol. 1, 1973, pp. 73. 35. Timoshenko, S.P., and Goodier, J.N., Theory of Elasticity, 3rd edition, McGraw H i l l , New York, 1970. 36. Tong, P., "New Displacement Hybrid Finite Element Model for Continua", International Journal for Numerical Methods in Engineering, Vol. 2, 1970, pp. 73-83. 37. Tong, P., and Pian, T.H.H., "The Convergence of Finite Element Method in Solving Linear Elastic Problems", International Journal of Solids Structures, Vol. 3, 1967, pp. 865-879. 38. Tong, P., and Pian, T.H.H., "On the Convergence of the Finite Element Method for Problems with Singularity", International Journal of Solids Structures, Vol. 9, 1973, pp. 313-321. 39. Watwood, V.B., Jr., "The Finite Element for Prediction of Crack Behaviour", Nuclear Engineering Design II, 1969, pp. 323-332. 40. Wiinderlich, W., "Discretisation of Structural Problems by a Generalized Variational Approach", Proceedings of the International Association for Shell Structures, Pacific Symposium on Hydrodynamically Loaded Shells— Part I, Honolulu, Hawaii, October 10-15, 1971. 41. Zienkiewicz, O.C, The Finite Element Method in Engineering Science, McGraw H i l l , London, 1971. 145 Zienkiewicz, O.C. , Cheung, Y.K., and Stagg, K.G., "Stresses in Aniso-tropic Media with Particular Reference to Problems of Rock Mechanics", Journal of Strain Analysis, Vol. 1, No. 2, 1966, pp. 172-182. Balakrishnan, A.V., Applied Functional Analysis, Applications of Mathemati 3, Springer-Verlag, New York, 1976. Degrees of Sign of No. of Interpolation Freedom Eigen- Eigen- Composition of Eigenvectors j u,v x's* u,v x's Total values values (-) 3 X 's constant; u,v linear. Linear Constant 6 3 9 (0) (+) 3t 3 X X 's=u, =v, =u, +v, =0.t x y y x 's constant; u,v linear. (-) 9 X 's,u,v linear. Linear Linear 6 9 15 (0) (+) 3 3 X X 's=u, =v, =u, +v, =0. | x y y x j 1s,u,v linear. (-) 3 X 's constant; u,v quadratic. Quadratic Constant 12 3 15 (0) (+) 9** 3 X X 's=0;{u, =v, =u, +v, =0;u,v quadratic}, j x y y x 's constant; u,v quadratic. (-) 9 X 's linear; u,v quadratic Quadratic Linear 12 9 21 (0) (+) 3 9 X X s=u, =v, =u, +v, =0. x y y x 's linear; u,v quadratic. (-) 18 X 's,u,v quadratic. Quadratic Quadratic 12 18 30 (0) (+) 3 9 X X 's=u, =v, =u, +v, =0. x y y x s,u,v quadratic. *x's: A l l stresses x ,x and x have the same type of interpolation. xx yy xy **Extra zero eigenvalues are associated with mechanisms which have the same u,v distributions as the approximating polynomials. tRigid body modes (u, =0;v, =0;u, +v, =0). x y y x TABLE I: Eigenvalues and eigenvectors of element matrix for triangular elements using different combinations of interpolations for u,v and x's; linear elasticity plane stress. Interpolation Degrees of Freedom Sign of Eigen-values No. of Eigen-values Composition of Eigenvectors u ,v X ,x ,x xx yy xy u,v x's Total Bilinear Constant 8 3 11 (-) (0) (+) 3 3 x's constant; u,v bilinear. x's=0;{u, =v, =u, +v, =0;u,v bilinear.} x y y x x's constant; u,v bilinear. Bilinear Bilinear 8 12 20 (-) (0) (+) 12 3t 5 x's,u,v bilinear. x's=u, =v, =u, +v, =0.t x y y x x's,u,v bilinear. * Biquadratic; Bilinear 16 12 28 (-) (0) (+) 12 4** 12 x's bilinear; u,v biquadratic. x's=0;{u, =v, =u, +v, =0;u,v biquadratic}, x y y x x's bilinear; u,v biquadratic. I Biquadratic Biquadratic 16 24 40 (-) (0) (+) 24 3 13 x's,u,v biquadratic. x's=u, =v, =u, +v, =0. x y y x x's,u,v biquadratic. *Full quadratic in x and y plus x zy and xy . **Extra zero eigenvalues are associated with mechanisms which have the same u,v distributions as the approximating polynomials. tRigid body modes (u, =0;v, =0;u, +v, =0). ° J x y y x TABLE II; Eigenvalues and eigenvectors of element matrix for rectangular elements using different combinations of interpolations for u,v and x's; linear e l a s t i c i t y plane stress. Degrees of Sign of No. of 1 Interpolations Freedom Eigen- Eigen- Composition of Eigenvectors u,v P T'S u,v P x's Total values values (-) 10 p const.;U,V,T's linear. 2++3+0 3 p linear;u=v=x's=0. Linear Linear Linear 6 3 9 18 (0) (+) p=T s=u, =v. =u, +v, =0. x y y x p const.;U,V,T's linear. (-) 12 p,x's linear;U,V quadratic. Quadratic Linear Linear 12 3 9 24 (0) (+) 0+3*+0 9 P=T'S=U, =V, =U, +V, =0. x y y x p,x's linear;U,V quadratic. (-) 6 T'S const.; p,u,v quadratic. p quadratic;U=V=T's=0.^ f ] Quadratic Quadratic Constant 12 6 3 21 (0) (+) 3+3+6# 3 (P=T'S=U, =v, =U, +V, =0/. 1 x 'y 'y 'x )' p=x's=0;u,v quadratic. T'S const.;p,u,v quadratic. (-) 21 p linear; U,V,T'S quadratic. Quadratic Linear Quadratic 12 3 18 33 (0) (+) 3 9 P=T'S=U, =V, =U, +V, =0. I v » x 'y 'y 'x p linear;U,V,T's quadratic. (-) 21 P,T'S,U,V quadratic, p quadratic;U=V=T's=0. Quadratic Quadratic Quadratic 12 6 18 36 (0) (+) 3+3+0 9 p=T s=u, =v, =u, +v, =0. x y y x p,TTs,u,v quadratic. tSelf-equilibrating modes in pressure. *Rigid body modes (u, =0;v, =0;u, +v, =0). //Mechanisms. x y y x Note: The number of zero eigenvalues appear in the order; self-equilibrating, r i g i d body and M mechanisms respectively. Proper rigid body modes also have P=T'S=0. °° TABLE III: Eigenvalues and eigenvectors of element matrix for triangular elements using different combinations of interpolations for u,v,p and T'S (T ,T ,T ); two-dimensional, r r xx yy xy incompressible creeping flow (linear part of the Navier-Stokes equations). r Interpolations Degrees of Freedom Sign of Eigen-values (-) (0) (+) No. of Eigen-values 15 1^+3+0 5 Composition of Eigenvectors u,v P x's u,v P x's Total B i l i n e a r B i l i n e a r B i l i n e a r 8 4 12 24 p,x's,u,v b i l i n e a r , p bilinear;x's=u=v=0. p=x s=u, =v, =u, +v, =0. x y y x p,x's,u,v b i l i n e a r . Biquadratic B i l i n e a r B i l i n e a r 16 4 12 32 (-) (0) (+) 16 0+3+1// 12 p,x's b i l i n e a r ; u , v biquadratic. p=x's=u, =v, =u, +v, =0. { x y y x } p=x's=0;u,v biquadratic. p,x's b i l i n e a r ; u , v biquadratic. Biquadratict B i l i n e a r Biquadratic 16 4 24 44 (-) (0) (+) 28 0+3*+0 13 p b i l i n e a r , x ' s , u , v biquadratic. p=x's=u, =v, =u, +v, =0. x y y x p b i l i n e a r ; x ' s , u , v biquadratic. Biquadratic Biquadratic Biquadratic 16 8 24 48 (-) (0) (+) 26 2+3+0 13 p,x's,u,v biquadratic. p biquadratic;x's=u=v=0. p = x s=u > x=v, y=u, y+v, x-0. p,x's,u,v biquadratic. t S e l f e q u i l i b r a t i n g modes i n pressure. *Rigid body modes (u, x=v, y=u, y+v, x=0). //Mechanisms. Note: The number of zero eigenvalues appear i n the order: s e l f - e q u i l i b r a t i n g , r i g i d body and mechanisms, respectively. Proper r i g i d body modes also have p=x's=0. TABLE IV: Eigenvalues and eigenvectors of element matrix for rectangular elements using d i f f e r e n t combinations of inte r p o l a t i o n s for u,v,p and x's (x ,x ,x ): two-dimensional incompressible r xx yy xy creeping flow ( l i n e a r part of the Navier-Stokes equations). 150 j N ° -I o f 1026\„EI M 1029_EI hi i o v E 103UEI I I Elem. q l 3 q i q 2-f 5 I j 2 1.041667 2.083333 2.500 2.604167 1 4 1.236979 3.515625 3.750 3.743489 1 6 1.273148 3.858025 4.167 3.975909 8 1.285807 3.987630 4.375 4.058838 1 10 1.291667 4.050000 4.500 4.097500 I 12 1.294850 4.084680 4.583 4.118575 1 14 1.296769 4.105928 4.643 4.131308 I 16 1.298014 4.119873 4.688 4.139582 I 18 1.298868 4.129515 4.722 4.145260 [ 20 1.299479 4.136458 4.750 4.149333 J EXACT 1.302083 4.166667 5.000 4.166667 TABLE V(a): Simply supported beam; moments at the nodes are exact. No. of Elem. 1C 2 9 L E E I 1026 El M 1 0 e M E I ' 1 0 6 R E E I l o e ^ E i 102UEI q l 3 q l 4 q l 3 q l 3 2 — 4. 687500 — 1. 354167 1.770833 2 .864583 4 5. 33850 4. 492190 1. 601563 1. 276042 1.705730 2 .587891 6 3. 74228 4. 456020 1. 554784 1. 261574 1.685957 2 .538795 8 2. 88086 4. 443360 1. 531576 1. 256510 1.678060 2 .521770 10 2. 34167 4. 437500 1. 517500 1. 254167 1.674167 2 .513917 12 1. 97242 4. 434320 1. 508005 1. 252894 1.671971 2 .509659 14 1. 70372 4. 432400 1. 501154 1. 252126 1.670616 2 .507093 16 1. 49940 4. 431150 1. 495972 1. 251628 1.669723 2 .505430 18 1. 33887 4. 430300 1. 491912 1. 251286 1.669095 2 .504290 20 1. 20938 4. 429690 1. 488646 1. 251042 1.668646 2 .503474 EXACT 0. 00000 4. 427083 1. 45833 1. 250000 1.666667 2 .500000 | TABLE V(b): Cantilever; moments at the nodes are exact. TABLES V: Numerical results for two second order beam equations; v-linear, M-linear. Rotations and shears are derived from v and M, respectively. 151 No. of j Elem. 103617EI E 1036 EI M -10 ^ 104UEI j q l 2 q l 2 q 2j5 1 2 — 2.60417 6.25000 6.250 6.510417 | 4 5.859 2.60417 4.68750 7.813 6.917320 6 4.822 2.60417 4.39815 8.102 6.939086 8 3.988 2.60417 4.29688 8.203 6.942750 10 3.375 2.60417 4.25000 8.250 6.943750 12 2.918 2.60417 4.22458 8.275 6.944111 14 2.566 2.60417 4.20918 8.291 6.944264 16 2.289 2.60417 4.19922 8.301 6.944338 18 2.065 2.60417 4.19239 8.308 6.944380 j 20 1.880 2.60417 4.18750 8.313 6.944401 J EXACT 0.000 2.60417 4.16667 8.333 6.944444 TABLE V(c): Beam with both ends fixed, deflections at the nodes are exact. No. of Elem. 102 S L E E I 1C 2 6 M E I M 1 0 2 6 R E E I 1 0 2 e R E E I -1 0 M L E 1 0 MRE 102UEI q l 3 q l 4 q l 4 q l 3 q l 2 q l 2 q 2 i 5 2 — 2. 34375 4. 16667 — 3 .12500 1.87500 1.1067708 4 3. 1901 2. 34375 4. 16667 2.0182 3 .28125 1.71875 1.1108398 6 2. 3341 2. 34375 4. 16667 1.3696 3 .31019 1.68982 1.1110575 8 1. 8305 2. 34375 4. 16667 1.3467 3 .32031 1.67969 1.1110942 10 1. 5042 2. 34375 4. 16667 0.8292 3 .32500 1.67500 1.1111042 12 1. 2756 2. 34375 4. 16667 0.6920 3 .32755 1.67245 1.1111078 14 1. 1070 2. 34375 4. 16667 0.5937 3 .32909 1.67092 1.1111093 16 0. 9776 2. 34375 4. 16667 0.5198 3 .33008 1.66992 1.1111101 18 0. 8752 2. 34375 4. 16667 0.4623 3 .33076 1.66924 1.1111104 20 0. 7922 2. 34375 4. 16667 0.4162 3 .33125 1.66875 1.1111107 j EXACT 0. 0000 2. 34375 4. 16667 0.0000 3 .33333 1.66667 1.1111111 TABLE V(d): Beam with one end fixed and the other guided; deflections at the nodes are exact. TABLES V: Numerical results for two second order beam equations; v-linear, M-linear. Rotations and shears are derived from v and M, respectively. 152 fi No. of Elem. 10*6^1 10 2 e EI h 1036 QEI 10 3UEI q l * q l 3 1 1.250000 5.000000 — 4.16667 2 1.302083 4.791667 9.244792 4.16667 3 1.301440 4.506173 — 4.16667 4 1.302083 4.375000 9.277344 4.16667 5 1.302000 4.306667 — .4.16667 6 1.302083 4.266975 9.276942 4.16667 7 1.302062 4.241983 — 4.16667 8 1.302083 4.225260 9.277344 4.16667 15 1.302082 4.184197 — 4.16667 16 1.302083 4.183129 9.277344 4.16667 1 EXACT 1.302083 4.166667 9.277344 4.16667 j TABLE VI(a): Simply supported beam. No. of 1 Elem. 102 9. L EEI 10 26 EI M 106 EI M 1 0 6 R E E I 1 0 9 R E E I 102UEI q l 3 q l 4 q l 3 q l 3 q 2 2 5 1 5.0000 4.375000 1.250000 1.25 2.000000 2.50 2 1.6667 4.427083 1.562500 1.25 1.708333 2.50 3 0.8025 4.426440 1.435185 1.25 1.679012 2.50 4 0.4688 4.427083 1.484375 1.25 1.671876 2.50 5 0.3067 4.427000 1.450000 1.25 1.669333 2.50 6 0.2161 4.427083 1.469907 1.25 1.668210 2.50 7 0.1604 4.427062 1.454082 1.25 1.667637 2.50 8 0.1237 4.427083 1.464844 1.25 1.667318 2.50 15 0.0360 4.427082 1.457407 1.25 1.666763 2.50 16 0.0317 4.427083 1.459961 1.25 1.666747 2.50 | EXACT 0.0000 4.427083 1.458333 1.25 1.666667 2.50 TABLE VI(b); Cantilever TABLES VI: Numerical r e s u l t s f o r two second order beam equations; v-quadratic, M-quadratic. Rotations and shears are derived. Moments and shears are exact i n a l l cases. 153 No. of Elem. 1O30 EI E 1036 QEI 10 39 QEI 1036 EI M lO^UEI ql 3 ql1* qi 3 ql 4 q 2 ! 5 1 8.333 — — 2.083333 6.94444 2 6.250 1.432292 5.20833 2.604167 6.94444 3 3.395 — — 2.597737 6.94444 4 2.083 1.464844 9.11458 2.604167 6.94444 5 1.400 — — 2.603333 6.94444 6 1.003 1.464442 7.52315 2.604167 6.94444 7 0.753 — — 2.603950 6.94444 8 0.586 1.464844 8.13802 2.604167 6.94444 15 0.175 — — 2.604156 6.94444 16 0.155 1.464844 7.89388 2.604167 6.94444 EXACT 0.000 1.464844 7.81250 2.604167 6.94444 TABLE VI(c): Beam with both ends fixed. No. of Elem. i 0 2 e L E E i i o V 1 lO^JSI M 10 26 R EEI I O 3 G R E E I 102UEI ql 3 ql1* ql 3 ql1* ql 3 q 2 ! 5 1 5.0000 2.291670 4.16667 4.16667 33.3333 1.111111 2 1.6667 2.343750 7.29167 4.16667 04.1667 1.111111 3 0.8025 2.343107 6.01852 4.16667 01.2350 1.111111 4 0.4688 2.343750 6.51042 4.16667 00.5210 1.111111 5 0.3067 2.343667 6.16667 4.16667 00.2667 1.111111 6 0.2161 2.343750 6.36574 4.16667 00.1543 1.111111 7 0.1603 2.343728 6.20750 4.16667 00.0972 1.111111 8 0.1348 2.343750 6.27010 4.16667 00.0651 1.111111 15 0.0361 2.343749 6.24074 4.16667 00.0099 1.111111 16 0.0317 2.343750 6.26628 4.16667 00.0082 1.111111 EXACT 0.0000 2.343750 6.25000 4.16667 00.0000 1.111111 B^LE VI(d): Beam with one end fixed and the other guided. IBLES VI: Numerical results for two second order beam equations; v-quadratic, M-quadratic. Rotations and shears are derived. Moments and shears are exact in a l l cases. 154 No. of Elem. l O ^ E I 10 26 EI M 103UEI ql3 q^ q 2 ! 5 1 5.000000 — 4.166667 2 4.305556 1.307870 4.166667 3 4.209877 — 4.166667 4 4.185049 1.302594 4.166667 5 4.176092 — 4.166667 6 4.172123 1.302189 4.166667 EXACT 4.166667 1.302083 4.166667 TABLE VII(a) : Simply supported beam. Mo. of Elem. 10 39 T JEI Lh, 1026\,EI M 108WEI M 1 0 6 R E E I loe^Ei 102UEI J ql3 q^ ql3 ql4 ql3 q 2! 5 1 8.3333 — — 1.250 1.583333 2.50 2 1.3889 4.432870 1.458333 1.250 1.652778 2.50 3 0.4321 — — 1.250 1.662346 2.50 4 0.1838 4.427594 1.458333 1.250 1.664828 2.50 5 0.0943 — — 1.250 1.665724 2.50 6 0.0546 4.427189 1.458333 1.250 1.666121 2.50 EXACT 0.0000 4.427083 1.458333 1.250 1.666667 2.50 TABLE VII(b): Cantilever TABLES VII: Numerical results for two second order beam equations, v-cubic, M-cubic. Moments and shears are exact in a l l cases. 155 No. of Elem. 1039T,EI E 10 36 EI M lO^UEI ql * q 2! 5 1 8.3333 — 6.94444 2 1.3889 2.662037 6.94444 3 0.4321 — 6.94444 4 0.1838 2.609273 6.94444 5 0.0943 — 6.94444 6 0.0546 2.605228 6.94444 EXACT 0.0000 2.604167 6.94444 TABLE VII(c): Beam with both ends fixed. No. of Elem. io 3 e L EEi i o ' V 1 1029 EI M 10 26 R EEI I O 3 S R E E I 102UEI q l 3 q l * q l 3 q l * q l 3 q 2 ! 5 1 8.3333 — — 4.166667 -8.3333 1.111111 2 1.3889 2.349537 6.250 4.166667 -1.3889 1.111111 3 0.4321 — — 4.166667 -0.4321 1.111111 4 0.1838 2.344261 6.250 4.166667 -0.1838 1.111111 5 0.0943 — — 4.166667 -0.0943 1.111111 6 0.0546 2.343856 6.250 4.166667 -0.0546 1.111111 EXACT 0.0000 2.343750 6.250 4.166667 0.0000 1.111111 CABLE VII(d): Beam with one end fixed and the other guided. CABLES VII: Numerical results for two second order beam equations, v-cubic, M-cubic. Moments and shears are exact in a l l cases. 156 No. of Elem. 10 26 EI M i o2 e E I E 10 3UEI 1 ql 1* q l 3 1 2 1.852 5.5560 4.6300 3 — 4.5730 4.2337 4 1.360 4.5140 4.1956 5 — 4.3380 4.1766 6 1.338 4.3210 4.1723 j 1 EXACT 1.302 4.1667 4.1667 j TABLE VI I I ( a ) : Simply supported beam; nodal moments and shears are exact. No. of Elem. 1 0 9 L E E I ^ V 1 1 0 6 R E E I 1 0 6 R E E I 102UEI q l 3 q l 4 q l 3 q l 4 q i 3 q 2 ! 5 1 1.6667 — — 1.6667 1.6667 4.1667 2 0.5555 3.9350 1.2500 1.2500 1.9444 2.5463 3 0.2309 — — 1.2551 1.5387 2.5324 4 0.1389 4.4840 1.5625 1.2500 1.7361 2.5029 5 0.0844 — — 1.2507 1.6169 2.5043 6 0.0617 4.3470 1.4352 1.2500 1.6975 2.5006 EXACT 0.0000 4.4271 1.4583 1.2500 1.6667 2.5000 TABLE VIII(b): Cantilever; nodal moments and shears are exact. TABLES VIII: Numerical r e s u l t s f o r four f i r s t order equations, forced boundary conditions on v and M; v,0,M and V a l l l i n e a r . 157 No. of Elem. 1O20„EI h, 1036 KI M 10\ 10 ^ lO'+UEI q l 3 q l * q l 2 q l 2 q 2! 5 1 2 3 4 5 6 1.3889 0.9150 0.3472 0.3442 0.1543 4.62960 3.18290 2.57210 8.3333 4.1667 4.6296 -8.3333 -8.2305 -8.3333 -8.3333 -8.3333 -1.1574 8.4675 7.2338 7.1547 7.0016 EXACT 0.0000 2.60417 4.1667 -8.3333 6.9444 TABLE VIII(c): Beam with both ends fixed; shears are exact. No. of Elem. 102 6 L E E I 1C 2 6 M E I M 1C 2 6 R E E I 102 6RE E I " 1 0 M L E 1 0 MRE ' 102UEI q l 3 q l * q l * q l 3 q l 2 q l 2 q 2! 5 1 4. 1666 — 4. 16667 4. 1667 2.5000 2.50000 1.0466 2 5. 5555 2. 54629 4. 16667 2. 7778 3.3333 1.66667 1.1574 3 1. 3775 — 4. 16667 0. 4515 3.3230 1.67695 1.1263 4 1. 3888 2. 40158 4. 16667 0. 6944 3.3333 1.66667 1.1141 0. 5108 — 4. 16667 0. 1775 3.3320 1.66800 1.1132 6 0. 6172 2. 34050 4. 16667 0. 3086 3.3333 1.66667 1.1117 EXACT 0. 0000 • 2. 34375 4. 16667 0. 0000 3.3333 1.66667 1.1111 [ TABLE VIII(d):' Beam with one end fixed and the other guided; shears are exact. TABLES VIII: Numerical results for four f i r s t order equations, forced boundary conditions on v and M; v,0,M and V a l l linear. 158 No. of Elem. (1) Simply Supp. 103UEI (2) Cantilever 102UEI q 2 ^ (3) Fixed-fixed lO^UEI q 2 ^ 1 2 3.472222 1.388884 __ 3 4.034536 2.194787 . — 4 4.123264 2.365451 6.510467 5 4.133495 2.440016 6.521569 6 4.158093 2.467577 6.858711 8 4.163954 2.488878 6.917318 10 4.165556 2.495289 6.933333 12 4.166131 2.497662 6.939086 16 4.166497 2.499243 6.942750 j EXACT 4.166667 2.500000 6.944444 TABLE IX: Strain energy estimations for four f i r s t order beam equations with forced boundary conditions on v and 9; M,V,9 and v a l l linear. |No. of Elem. 1026 EI M 1C 2 e E E i 1 0 MM 103UEI q l 3 q i 2 q l 2 q 2 ! 5 1 2 2.43056 4. 16667 0.83333 0 .08333 6 .076389 3 — 4. 16816 — 0 .00953 6 .762226 4 2.17014 4. 16667 1.45833 0 .02083 6 .727431 5 — 4. 16793 — 0 .00397 6 .769543 6 2.13049 4. 16667 1.20370 0 .00926 6 .762251 8 2.10504 4. 16667 1.30208 0 .00521 6 .768121 10 2.10056 4. 16667 1.23333 0 .00333 6 .769722 12 2.09298 4. 16667 1.27315 0 .00231 6 .770298 16 2.08876 4. 16667 1.26302 0 00130 6 .770664 EXACT 2.08333 4. 16667 1.25000 0 00000 6 .770833 TABLE X(a): Simply supported beam, shears are exact. TABLES X: Numerical results for four f i r s t order beam equations; forced boundary conditions on v and 6; shear strain energy included; M,V,e,v a l l linear. 159 No. of Elem. 1 0 6 R E E I 1 0 9 R E E I 1 0 M L E VLE 102UEI q l * q i 3 q l 2 qi q 2 ! 5 1 0.84821 1.07143 1. 07143 1. 35714 3 .794643 2 1.43430 1.66667 3. 39744 1. 46154 3 .151709 3 1.54113 1.63422 4. 33756 1. 17524 3 .487009 4 1.55303 1.66667 4. 56439 1. 13636 3 .513652 5 1.55948 1.66195 4. 74883 1. 07071 3 .533663 6 1.56056 1.66667 4. 80288 1. 06272 3 .535970 8 1.56188 1.66667 4. 88839 1. 03571 3 .539845 10 1.56225 1.66667 4. 92835 1. 02299 3 .540917 12 1.56238 1.66667 4. 95016 1. 01601 3 .541304 16 1.56246 1.66667 4. 97192 1. 00904 3 .541552 EXACT 1.56250 1.66667 5. 00000 1. 00000 3 .541667 TABLE X(b): Cantilever. No. of Elem. 10 26 EI M 1 O 2 M M 10 2M^ 103UEI q l * q l 2 q i 2 q 2! 5 1 2 1.04167 0.00000 0.00000 2.604167 3 — — 6.51466 3.242249 4 1.12847 6.25000 6.25000 3.255208 5 — — 7.58171 3.289003 6 1.05024 3.70370 7.40741 3.290038 8 1.06337 4.68750 7.81250 3.295898 10 1.04500 4.50000 8.00000 3.297500 12 1.05131 4.39815 8.10185 3.298075 16 1.04709 4.29688 8.20313 3.298442 EXACT 1.04167 4.16667 8.33333 3.298611 1 TABLE X(c): Beam with both ends fixed; shears are exact. TABLES X: Numerical results for four f i r s t order beam equations; forced boundary conditions on v and 6; shear strain energy included; M,V£,v a l l linear. FINITE ELEMENT 10Etu_ 102Etu lOEtv 10Etv D -ION , | xxA ) N f )L (1-v Z)N nL ( l - V -)N0L (l-v ; -)N L Nn 1 GRID N Elem. A* Elem. Bt Elem. A Elem. B Elem. A Elem. B Elem. A Elem. B Elem. A Elem. B J l x l 2 1.054258 1.507941 -1.1435 2.1934 2.99336 1.31824 4.332647 5.085466 0.00000 1.44190 2 x 2 4 1.443990 1.519821 -0.7670 1.8684 1.45335 1.28574 4.849820 5.073595 2.01236 1.40137 3 x 3 6 1.463880 1.519773 0.6555 1.8046 1.34960 1.27936 4.974455 5.073633 1.41893 1.40559 4 x 4 8 1.498170 1.519862 1.3373 1.7896 1.28864 1.27787 5.016690 5.073544 1.58424 1.40789 5 x 5 10 1.508706 1.519900 1.5268 1.7852 1.28480 1.27742 5.036405 5.073507 1.28638 1.40880 6 x 6 12 1.512260 1.6488 1.27514 5.047228 1.49870 EXACT 1.519928 1.519928 1.7837 1.7837 1.27727 1.27727 5.073478 5.073478 1.40954 1.40954 j ION . yyA ION yyB 102N xxB 1 0 NxxC 10Et2U Nn No N 0 N 0 ( l - v 2 ) L 2 N 2 N Elem. A Elem. B Elem. A Elem. B Elem. A Elem. B Elem. A Elem. B Elem. A Elem. B 2 6.66667 8. 55810 4.04167 4.70735 6.2500 3.9928 0.0000 0.3181 2.553610 2.787981 4 9.16957 8. 59863 3.54268 4.17500 0.9848 0.4235 0.7266 -0.0005 2.746331 2.793366 6 7.75831 8. 59441 4.31117 4.11902 4.8662 0.0848 -0.4544 -0.0299 2.782543 2.793540 8 8.92269 8. 59211 4.02879 4.10971 3.0654 0.0329 -0.3742 -0.0285 2.789962 2.793562 10 8.39253 8. 59120 4.05786 4.10767 1.9439 0.0166 -0.2945 -0.0233 2.792055 2.793567 12 8.74161 4.13723 1.5472 -0.2530 2.792746 EXACT 8.59046 8. 59046 4.10670 4.10670 0.0000 0.0000 0.0000 0.0000 2.793570 2.793570 *Mixed Finite Element; displacements and stresses linear. tDisplacement Finite Element; u and v f u l l cubics [5]. TABLE XI: Numerical results for parabolically loaded plane stress problem. 1 No. of Elem. 8 in Beam Depth N Total No. of Elem. Degrees of Freedom Tip Deflec-t i o n 6 c v(L,0).(in.) T @ X X x=12",y=-6" k s i Longitudinal D e f l e c t i o n @ B=UB. (in.) St r a i n Energy U(k-in.) 2 32 131 0.336062 53.084 0.059717 6.731654 4 128 421 0.355121 60.501 0.064024 7.136889 6 288 871 0.355093 59.612 0.063980 7.140213 8 512 1481 0.355459 60.148 0.064038 7.145934 EXACT (ELASTICITY) 0.355333 60.000 0.064000 7.146667 TABLE XII: Numerical re s u l t s for the ca n t i l e v e r (plane stress) with boundary conditions B.C.I (Figure 19). Mixed f i n i t e element; displacements and stresses l i n e a r . No. of Elem. Total Degrees Tip Deflec- T @ X X Longitudinal S t r a i n i n Beam No. of of t i o n 5 = c x-12",y=-6" Def l e c t i o n @ Energy 1 Depth N Elem. Freedom v(L,0).(in.) k s i B=UB. (in.) U(k-in.) I 1 8 46 0.245120 39.107 0.042410 4.902405 I 2 32 129 0.3359427 52.694 0.059602 6.718577 4 128 415 0.355464 60.469 0.064064 7.109509 6 288 861 0.355698 59.734 0.064087 7.114006 8 512 1467 0.355952 60.125 0.064146 7.119015 J [ EXACT (BEAM THEORY) 0.355833 60.000 7.116667 J TABLE XIII: Numerical results f o r the ca n t i l e v e r (plane stress) with boundary conditions B.C.2 (Figure 19). Mixed f i n i t e element; displacements and stresses l i n e a r . 162 Element Type Mesh N Number of Degrees of Freedom Tip Deflection 6=vc (in.) x - Normal X X Stress @ x=12" y=-6". (ksi) A-1 2 48 0.19819 33.407 j C.S.T A-2 4 160 0.30556 51.225 A-3 8 576 0.34188 57.342 B-1 1 48 0.34872 L.S.T B-2 2 160 0.35506 59.145 B-3 4 576 0.35569 60.024 j C-1 1 68 0.35373* 58.973* I Q.S.T C-2 2 214 0.35506 59.843 I C-3 4 268 0.35580 59.993 MIXED M-1 1 46 0.24512* 39.108 DISPL. M-2 2 129 0.33594 52.694 AND M-4 4 415 0.35547 60.469 STRESSES M-6 6 861 0.35570 59.734 I LINEAR M-8 8 1467 0.35595 60.125 1 j BEAM THEORY 0.35583 60.000 1 *Average of values at y=6" and y=-6". TABLE XIV: Comparison amongst C.S.T., L.S.T., Q.S.T., and mixed f i n i t e element. Cantilever with boundary conditions B.C.2 (Figure 19). E U A E V A E v c E UD TxxD T _ TxxA EU* T 0 L T0 T0 2 0.4633 1.4778 1.1222 0.2440 0.7052 1.3835 0.6897 2.738912 4 0.4571 1.7771 1.1737 0.4271 1.4203 2.3751 -0.7849 3.122624 6 0.4507 1.7810 1.1771 0.4787 1.8202 3.0010 0.5045 3.166420 8 0.4739 1.8312 1.1849 0.5216 2.1230 3.4100 -0.2607 3.196404 10 0.4510 1.8186 1.1856 0.5364 2.4049 3.8175 0.0071 3.203674 12 0.4610 1.8272 1.1869 0.5546 2.6225 4.1357 0.0469 3.212178 TABLE XV(a): x continuous at the crack tip D. N E U A E V A E v c E uD T xxD TyyD TxxA EU* x 0 L T Q L x 0L x 0 L T0 T0 T0 2 0.0544 2.9301 2.8883 -0.7285 -0.1440 2.2921 -0.7936 1.704874 4 0.4568 1.7749 1.1719 0.4281 1.4257 3.0437 -0.7605 3.113852 6 0.4461 1.7709 1.1762 0.4817 1.8396 3.8743 0.4800 3.159226 8 0.4732 1.8273 1.1837 0.5252 2.1493 4.3664 -0.2376 3.191367 10 0.4996 1.8147 1.1848 0.5404 2.4361 4.8837 -0.0143 3.199594 12 0.4600 1.8241 1.1862 0.5585 2.6562 5.2798 0.0688 3.208874 TABLE XV(b): x discontinuous at the crack tip D. yy TABLES XV: Numerical results for the plane stress problem of square plate with symmetric edge cracks, Figure 29. UE O N Co *Exact value: U = 3.228 2 t 2-. T° ng and Pian [38] T n L t 164 E7TM 106EATL, M aE Air . M K I Error % a T 2hbt T^hbt T 2hbt Aa T 0/a 0.0 1.05942309 1.385584 0.274918 1.90402 10.98 0.1 1.05942339 1.680112 0.333358 2.09664 1.97 0.2 1.05942331 1.607250 0.318898 2.05067 4.12 J 0.5 1.05942333 1.626610 0.322740 2.06300 3.55 I n i t i a l Crack 1.05942171 EXACT K ; ref. [2] 2.13884 TABLE XVI: Stress intensity factors from the f i n i t e element analysis of the rectangular plate with symmetric edge cracks, Figure 33(a). (Mixed f i n i t e element; displacements and stresses linear). r r 0 E*M 106EATT M aE ATT„ M K I Error % a T 2hbt T 2hbt T 2hbt a T 0/a 0.0 1. 04790510 1.408614 0.279486 1.91978 8.98 0.1 1. 04790543 1.730814 0.343416 2.12804 0.89 0.2 1. 04790535 1.645374 0.326464 2.07485 1.63 0.5 1. 04790537 1.664358 0.330230 2.08679 1.06 I n i t i a l Crack 1 1. 04790370 EXACT K ; ref. [3] 2.10922 TABLE XVII: Stress intensity factors from the fi n i t e element analysis of the rectangular plate with a central crack, Figure 33(b). (Mixed f i n i t e element; displacements and stresses linear). Author(s) Number of Elements Degrees of Freedom Accuracy of K Error % Type of Element Watwood [38] 470 956 2.00 Triangular and rectangular Anderson et al. [1] 1470 3000 0.14 Quadrilateral Present result 174 505 1.97 Mixed triangles.* Symm. edge cracks Present result 174 505 0.89 Mixed triangles. Central crack TABLE XVIII: Comparison of stress intensity factors obtained from energy release rate using different elements and procedures. *Plane strain mixed f i n i t e element; displacements and stresses linear. 166'-FIG. 2 ! Typical contour for evaluation of J-integral 167 FIG.3' Accommodation of crack extension Aa by advancing nodes on the path TQ . 0 x FIG.4= Node numbers and degrees of freedom for a triangular element. H u 4 , v 4 ^4 4 4 T xx,Tyy,Txy 4 ? 8 * u„v, 1 xx >Lyy>t xy 7 -o--o-5 U 3 . V 3 3 3 3 Txx,Tyy,Txy • ° 3 46 u 2 ,v 2 2 2 2 1 xx >l y y>lxy 0 FIG.5-- Node numbers and degrees of freedom for a rectangular element. 169 FIG.7 : Fo rces act ing on an inf in i tes imal beam element 170 v(o) = M(o)=0 vU) = MU)=0 Cased) Simply supported (S.S.) v(o)=0 MU) = 0 1===--A Case(2) Cantilever v(o) = 0 A vU)=0 § Case(3) Both ends clamped (fixed- fixed) v(o)=0 i l « ^ Case(4) One end clamped and the other in a vertical guide (fixed - guided) Deflected elastic curve FIG.8 ! Forced boundary conditions on v and M for the beam problem. 171 V l v 2 M i M o 6 X, ( i ) v - l i n e a r , M - l i n e a r v, v 2 v 3 M| M 2 M 3 I <5 o cW 2 (ii) v - quad ra t i c , M - q u a d r a t i c 0, 0 2 M, M 2 Vl v 2 i <5 J p (iii) v - cub i c , M - c u b i c FIG.9= Degrees of f r eedom for the beam e l e m e n t - t w o s e c o n d order e q u a t i o n s . 172 I 2 4 6 8 10 20 40 60 80 100 Number of element along beam length"N" FIG.10(a) Two second order beam equations ; displacement linear and moment linear. FIG.IO(c):Two second order beam equations; displacement linear and moment linear. 1.75 FIG.IO(d):Two second order beam equations; displacement linear and moment linear. _ Relative error O Mid and end O ,L rotations ' 9LT 1 7 7 6 1 2 4 6 Number of elements N FIG.I2! Two second order beam equations, displacement cubic and moment cubic. Relative error in strain energy Rotation at guided end even no.of elements Fixed end rotation .odd no. of elements 2 4 6 8 1 Number of elements FIG. I3<b)= Four first order beam equations .forced boundary conditions on v and M; v, M, 0 and V all linear. v-2 h-1 00 O Numberof elements "N" Number of elements FIG. I4: Four first order beam equations,forced boundary conditions on v and 0;v,0, M and V all linear. FIG.15(a) Four first order beam equations,forced boundary conditions on v and 0;shear strain energy. FIG.15(b) Cantilever FIG. 15 (c) Beam clamped at both ends. FIGS. 15 : Four first order beam equations, forced boundary conditions on v and 6,shear strain energy included. o ZD Q O CL 1_ o I ro ro o O C L o CD o C L o o ro CD r*» co .O J > — ' —> Q CO -1 CD T3 D r T3 O CO CD Q 3 5 CO CD - • 3 — Q O ro o O CD nv CO CD tre -» J > Crt CD o CO 3 ce C D O X o CO CD CL -\ T J O CD JO CD o 3 bo CD — 3 a —t- Q Relative error in strain energy ~~ O i O 1 Re 3 , < la ti ve error in disp \ c >lacement 5 »78T 185 FIG.I7(c): Stress convergence. FIG.17 ' Parabolicolly loaded plane stress problem using mixed element; displacements and stresses linear. 186 FIG. 18 : Pa ra bo Ii ca lly loaded plane stress problem using displacement element (cubic displacements ), Ref. C5]. 187 Tota l load P = 4 0 D h = l2 T o t a l l oad P= 4 0 k E = 3 0 0 0 0 ksi v = 0 . 2 5 t = l " FIG. 19(a) Cantilever beam and load system L P a r a b o l i c a l l y v a r y i n g end shear X U _h_ 2 h_ 2 P a r a b o l i c s h e a r = £ ( # - / ) u ( 0 , 0 ) = v ( 0 , 0 ) = u(0,-§-) = u(0,-4) = 0 •xy g l v 4 FIG. 19(b) Fixed end boundary conditions B.C.I. x u 2 2 u( 0 , y ) = v ( 0 , y ) = 0 FIG.19(c) Fixed end boundary conditions B.C. 2 . FIGS.I9: Linear elasticity cantilever problem and forced boundary conditions used . 4 6 8 10 Number of elements in beam depth FlG.20(a) FIG.20(b) FIG.20(c) FlG.20= Convergence plots for the cantilever with boundary conditions B.C.I, using mixed finite ele me nts ; stresses and displacements all linear. 189 FIG.2l(a) ! Grids used for the constant stress triangular elements FIG.2l(b) : Grids used for the linear stress tr iangles 191 ~6-MESH C -3 - 42 Q.S.T's N = 4 FIG.21 (c) : Grids used for the quadratic stress triangles 192 ^ 4N - Ele me nts MESH M-4 128 MIXED FINITE ELEMENTS Fl G.22: Ty p i c a I mesh used for mixed finite ele ment ( N = 4 ). 0 200 400 600 800 1000 1200 1400 1600 Degrees of freedom FIG.23 (a): Plot of tip deflection vs. total degrees of freedom. 193 0 I 2 3 4 5 6 7 8 9 10 N = Number of elements in beam depth FIG.23(b) ! Plot of strain energy versus number of elements in beam depth. 285 270 266.875 I 2 3 4 5 6 7 8 9 10 N = Number of elements in beam depth FIG.23(c)-Plot of tip deflection versus number of elements in beam depth. FIGS.23 : Plots for the cantilever with boundary conditions B.C. 2,using mixed finite element; stresses and displacements all linear. 194 FIG.24=Convergence of tip deflection for C.S.T., L.S.T. , Q.S.T.,and mixed finite element with stresses and displacements linear for boundary conditions B.C.2. 195 2L r = 8 unit thickness E. i.o. v . o . i . s-zjjTT) FIG.25(a) : Isotropic case. 2L unit thickness E, = 1.0, E2=3.0,z/,=0.l,i/2=0,G|2=0.42 FIG.25(b)' Orthotropic case. FIGS.25 : Model for infinite plate by a square plate with a circular hole at centre, isotropic and orthotropic. 196 FIG.26: Finite element grid for the square plate with a circular hole in the middle. Isotropic and orthotropic cases. 197 Exact solution infinite plate Finite element solution for FIG.27(a)- Isotropic,E= 1.0,v = O.I and G = 2(1+1/) FIG.27(b):Orthotropic,E,= I.O,E2=3.0, •v,=O.I ,v2=0.0 and G | 2=0.42 FIGS. 27 : Comparison of theoretical and finite element (mixed ) results for infinite plate with a circular hole in the middle. E x a c t s o l u t i on for i n f i n i t e p l a te o F i n i t e e l e m e n t s o l u t i o n FIG.28(a)' Isotropic FIG.28(b) = Orthotropic FIGS.28= Comparison of theoretica and constant stress finite element, Zienkiewicz [42 .results with the same material properties as used in the mixed method,* T q = -I.O. 198 2 L — t t t t i tv) 2 L L x(u) " yy= To i i i i i i FlG.29(a) : Square plate with edge cracks under uniaxial tension , v - 0.3 . rIG.29(b) Finite element ~~ dealization of the quarter slate along with boundary :onditions ( N = 4 ). L IGS.29s Plane stress problem considered for investigation of energy convergence in the case of stress s ingu lar i t ies . FIG.30(a): Tyy continuous at the crack tip D and FIG.30(b): r y y discontinuous at the comparison with results from other crack tip D. elements. ( Tong and PianC381 ). FIGS.30 : Plots of strain energy versus the mesh size for the plane stress problem with symmetric edge eracks -Figure 29. 200 O ^"yy stress continuous at the crack tip D. A T y y stress discontinuous at the crack tip D. FIG. 31 • Strain energy convergence for the plane stress problem square plate with symmetric edge cracks-Figure 29. 201 -2 - 1 •2-1 FIG.32( a)-Tyy continuous at the crack tip D. FIG.32(b): r y y discontinuous at the crack tip D . FIGS.32: Normal stress distribution along the middle of a square plate with symmetric edge cracks-Figure 29(b),(N= 12). 202 FIG.33(a): Symmetric edge cracks Fl G. 33(b): Centra I crack FIGS.33 ; Rectangular plates with cracks used for determining the crack intensity factor Kj. The quarter plate considered for the finite element analysis along with boundary conditions shown as shaded area . FIG.34 : Finite element mesh used for determining the crack intensity factor K I f u s e d for both symmetric edge and central cracks. 204 - 2 - 1 - 2 - 1 FIG. 35(a):Sy mmetric edge cracks. FIG. 35(b ): Centra I Crack. FIGS.35* Normal stress distribution along the edge OA of the rectangular plates with cracks - Figure 33. Boundary conditions' u ( L ,0 ) = v ( L,0 )= u( L,c) = u ( L,-c )= 0 FIG.36 : Boundary tractions and conditions for cantilever linear elasticity solution (Appendix D). ho o Ln 2Q6. APPENDIX A CALCULATION OF THE FINITE ELEMENT MATRIX; LINEAR DISPLACEMENTS AND STRESSES OVER A TRIANGULAR ELEMENT The element matrix equation f o r plane stress l i n e a r e l a s t i c i t y i s derived here f o r l i n e a r displacements and l i n e a r stresses i n a tr i a n g u l a r element using i s o t r o p i c material. Also demonstrated are the modifications necessary to a l t e r the element matrix f o r plane s t r a i n i s o t r o p i c and orthotropic materials. The element geometry i s shown i n Figure 4. Using area coordinates Pl>P2>P3» [43-]; the l i n e a r approximations f o r u. and T . . can be written i l as u± = < D 1 Pz P3> T± j = <Pl P2 P3> u i l u i 2 ui3 . J. W. 1=1,2. (A.l) i=j=l,2. (A. 2) i j -3 i j Note t e n s o r i a l notation i s implied. Comparing equations (A.l) and (A.2) with (4.1) gives *k = *k = p k ; k=l,2,3. Consider the mixed v a r i a t i o n a l p r i n c i p l e for homogeneous boundary conditions i n equation (5.13) F(A) = [A,A] A - 2JQ fTudft. = jn [2x TTu-x TCx]dfi - 2 / f Tudfi where T = < T T T > xx yy xy u = <u v> = <ui U 2 > " T f = <f f > ; x y < T 1 1 T22 T 1 2 > ; T (A. 3) (A. 4) (A. 5) (A.6) 207 T = 9_ 9x 0 f -9_ 9y 9_ 9y 9_ 9x (A. 7) and 1 -v -v 0 1 0 Substituion of (A.l) and (A. 2) into Q y i e l d s "a 0 " F(A) = T 2T 0 b b a _ Ci - ~T T (A.8) T -c - vc 0 -vc c 0 0 0 2(l+v)c In (A.9) the submatrices a,b,c,d,e,x and tl are given by a. . = / p .p . dfi; i=j=3. e f n p.p. dsl; i=j=3. Jcl ,y J C i j = ~ I k p i p j d f i ' i = J = 3 ' d. = /. f p.dfi; I 1 n x i e. = /. f p.dfi; l ' U y l 1=1,2,3. 1=1,2,3. < T 1 T2 T3 T l T2 T3 T l T2 T3 >T xx xx xx yy yy yy xy xy xy T u = <Uj U2 U3 v^ V £ V3: Now for s t a t i o n a r i t y 3F(A) 9u (A. 9) (A. 10) (A.11) (A. 12) (A.13) (A.14) (A.15) (A.16) (A.17) 208 ^ - 0 (A. 18) which leads to the following matrix equation: , T T b a - v c -vc 2(l+v)c (A.19) or EA = g. (A. 20) It is only a simple matter to evaluate the submatrices [a], [b] and [c]. ttl If (x^y^) are the coordinates of the i node of a triangular element (Figure 4), then the matrices a and b in terms of nodal coordinates are [a] = 3 x 3 [b] = \ 3 x 3 6 Y2-Y3 Y3-Y1 Y1-Y2 Y2-Y3 Y3-yi Y1-Y2 Y2-Y3 Y3-Y1 Y1-Y2 x 3 " x 2 x l " x 3 x 2 ~ x l x 3-x 2 x r-x 3 x 2-x x |_x3-x2 X l - x 3 x2-xj_ (A.21) (A.22) and the symmetric matrix c in terms of the area A of the triangle and modulus of elasticity E is given by 2 0 9 [c] = -3 x 3 1 2 E 2 1 1 1 2 1 . (A.2 3 ) 1 1 2 Since [c] i s symmetric, the matrix of c o e f f i c i e n t s E i n equation (A.2 0 ) i s also symmetric as expected. Next the degrees of freedom A are so arranged by interchanging the corresponding rows and columns of E such that A = <u: Vl T 1 T 1 T 1 Uj V O T 2 T 2 T 2 Uq Vo T 3 T 3 T 3 > T - x x xx yy xy z z w \ r i r 3 3 . 1 L r 2 2 T 2 xx yy xy - T3 xx yy xy (A.2 4 ) and the matrix E becomes R) 0 E 1 5 x 1 5 a l l 0 b n 0 0 a 2 1 0 b 2 1 0 0 a31 0 b31 0 b u a n 0 0 0 b 2 1 a 2 1 0 0 0 b31 a31 c l l f 11 o a 1 2 0 c12 f l 2 0 a13 0 c13 f l 3 0 c n 0 0 b 1 2 f l 2 c12 0 0 b13 f l 3 c13 0 8 l l b 1 2 a 1 2 0 0 §12 b13 a13 0 0 813 0 0 a 2 2 0 b 2 2 0 0 a32 0 b32 0 0 b22 a 2 2 0 0 0 b32 a32 c 2 2 f 22 0 a2 3 0 c23 f23 0 Symmetric c 2 2 0 0 b23 f 23 c2 3 0 §22 b23 a2 3 0 0 823 0 0 a33 0 b33 0 0 b33 a33 c33 f 33 0 c33 0 833 (A.2 5 ) where [f] = -v[c] and [g] = 2 ( l + v ) [ c ] . The corresponding entries i n the load vector p are also interchanged and the modified load vector becomes 210 {p} = <di ei 0 0. 0 d 2 e 2 0 0 0 d 3 e 2 0 0 0>T. (A.26) The equation (A.20) alters to EA = p (A.27) The boundary conditions associated with the plane elasticity problems considered in chapter 7 are of the type u. = 0 on S l u T..n. = 0 on S (A.28) 1 ] J T T . . n . = T? on Sm i j J l T when the body forces f and f are zero. Here S and S are the portions x y u x of the boundary S where the displacements and stresses are zero, respectively; and the portion on which the tractions T_^ are prescribed. The equations (A. 28) are similar to equations (5.18) i f u(? = c?=a=0 and S same as S . Thus 1 1 M T equation (A.19) is similar to the equations (5.40) and (5.41) except that the former is expressed in the matrix form. Therefore the sub-load vectors d and e in the element matrix equation arise from the boundary integral / T?(j>.ds wh ere the element boundary coincides with S and at present <d.=p ^rjj I X J. I K . since the same shape functions are used for u and v. Hence the derivation T T of the load vector <d e > is identical to i t s generation in the displace-ment method. The procedure for deriving the element matrix E in (A.25), outlined above, is quite general for plane elasticity. The only change that needs to be introduced in switching from isotropic plane stress to isotropic or orthotropic plane strain l i e s in the compliance matrix C of (A.8). For isotropic plane strain, the compliance matrix C is given by 1 - 7 * - 0 l-\> l-v 2 0 0 T-=- , l-v J (A.29) 211 and for orthotropic plane strain C = (1-nv2) -v1(l+vz) - v i ( l + v 2 ) ± (l-vf) 0 1 '12 where n=^?'> elastic constants E 1 } v i and G 1 2 are associated with behaviour normal to the plane of strata; and the elastic constants E 2,v 2 and G 2 1 (G 2 1 is independent here) with the plane of strata, as shown in Figure 25(b). 212 APPENDIX B EIGENVALUE DISTRIBUTION FOR A REAL SYMMETRIC MATRIX A quadratic form in n-variables, x i , X 2 , . . . x^, is an expression of the type n n m T Q = I E a..x.x. = x ax (B.l) 1-1 J - l 1 J 1 J " " where a is a symmetric matrix of constants. If the coefficients a., and the variables x. are restricted to real values, then Q is real. Let v be I the rank of matrix a. Now there exists a nonsingular transformation x = ty, (B.2) Strang [32], such that the coefficients t can be chosen so that Q reduces to Q = yf+yf+ • • • + y i - y i + r • • • -y2.- ( B - 3 ) Equation (B.3) is called the canonical form of the quadratic form Q. The number of positive terms in (B.3), denoted by I, is called the index of the quadratic form. It is determined uniquely by the matrix a. The types of a quadratic form are determined by the rank r and the index I, as follows: (a) Q i s positive definite i f , and only i f , I=n. (b) Q is positive semidefinite i f , and only i f , I=r<n. (c) Q is negative definite i f , and only i f , 1=0 and r=n. (d) Q is negative semidefinite i f , and only i f , 1=0 and r<n. (e) Q i s indefinite i f , and only i f , 0<I<r. These conditions are obvious from the canonical form in (B.3). Further the types of the matrix a associated with the quadratic form Q in (B.l) corres-pond to the types of the quadratic form, i.e. the matrix a associated with the positive definite quadratic form is also positive definite. 213 Since the matrix of coefficients (4.19) of the matrix equation (4.18) is symmetric indefinite, therefore attention shall be focussed on the quadratic form of type (e) above. The quadratic form associated with the matrix (4.19) is T T T T Q = 2u av - v bv = <u v > 0 a u T (B.4) a -b _ V where - is a symmetric positive definite matrix, and - is a real n A n • mxn rectangular matrix. Assuming that the rank r of the real symmetric matrix 0 a l -j r i s m+n, then how many positive and negative real eigenvalues does a D this matrix possess? This necessitates determination of index I of the quadratic form (B.4). Consider the negative of the quadratic form in (B.4) T T Ql = -Q = <v u > p T - i r .. b - a V • — - a O 1 u (B.5) The (m+n)x(m+n) square matrix of equation (B.5) can be written as b l l b12 • • • b i n " a l l - a 2 i . . • _ aml b21 b 2 2 . . . b 2 n - a 1 2 - a 2 2 . . • _am2 • • • * • • • • • bnl bn2 • . . b nn ~a l n ~ a2n • • . -a mn -an " a12 • • # - a l n - a 2 i - a 2 2 . • -- a2n 0 . • • • • MxN _ aml - am2 • . .-a mn (B.6) The principal minors for the matrix (B.6) can be written as 214 Mi ; M 2 = b l l bi 2 ; M 3 = b l l b 1 2 B 1 3 B 2 1 b 2 2 B 2 1 B 2 2 B 2 3 B 3 1 B 3 2 B 3 3 etc. Since the submatrix b is positive definite, therefore a l l principal minors up to M r are positive. The M N + 1 principal minor as determinant of the matrix »11 '12 »21 D 2 2 3 n l Dn2 [_-an - a 1 2 b l n ~ a l l b 2 n _ a 1 2 nn 1 1 1 . -a ln (B.7) which has a zero on the diagonal, cannot be positive. This holds because B is not positive definite. Therefore M ,-, is either negative or zero. -n+1 r n+1 The latter cannot be true since i f m=l, i.e. only one degree of freedom in u; then the quadratic form is positive semidefinite which is not true. Hence M , 1 is negative. As for the principal minors M ,„ to M , , these n+1 ° r n+2 n+m can be either positive, zero or negative. However for rank r=m+n no two consecutive M^ 's can be zero; i f M^ and a r e z e r o then the rank of the matrix (B.6) is k or less. Further, any zero in an i n d i c i a l sequence lies between adjacent terms with opposite signs. Therefore the i n d i c i a l sequence can be arranged as (l,Mi,M2 . . . ,M ,M , .. , . . .,M , ). The index I of the quadratic form equals the n n+1 n+m number of permanences of sign in any i n d i c i a l sequence where any zero entry is given an arbitrary sign. Therefore the index I for Q, in (B.5) is n since the number of permanences in the i n d i c i a l sequence above is n; i.e. the signs of principal minors Mj to M^ are positive. Now the quadratic form Q, i n (B.5) through nonsingular transformation of the type (B.2) can 215 be expressed in the form of (B.3) as Ql = -Q = v2+v2+ . . . +v2-u2-u?- . . . -u 2; (B.8) 1 ^ n J- ^ m or Q = u2+u2+ . . . +u 2-v 2-v 2- . . . -v 2. (B.9) Hence i t can be deduced from (B.9) that the matrix (4.19) has m positive and n negative real eigenvalues. Next, consider the fi n i t e element matrix equation (6.36) for the linear part of the Navier-Stokes equations. This is expressed in a slightly different form as follows: 0 a £ u T a 0 0 P = 0 (B.10) / 0 "I. T where y i s a symmetric, positive definite (nxn) matrix. A rearrangement of (B.10) yields 0 3 a u SA = 3 T "Y 0 T = 0 T a 0 0 P-where T T T A = <u x p > (m+n+J)x(m+n+l) 0 3 a T 3 Y 0 T a 0 0 (B.ll) (B.12) (B.13) and a and 3 are (mxj) and (mxn) rectangular matrices, respectively. Let the rank r of the matrix S be (m+n+l). The matrix S is symmetric and indefinite, therefore a l l i t s eigenvalues are real and i t remains to be determined as to how many of these are positive or negative. The quadratic form of the mixed variational principle associated with (B.ll) is 216 Q = 2U TBT - T T Y T + 2u Tap. ( B . 1 4 ) The f i r s t two terms on the right hand side are i d e n t i c a l to the quadratic form i n ( B . 4 ) . Further, i t i s possible to diagonalize the submatrix 0 6] B T - Y through a transformation of the type (B.2) as was done for (B..6) Let such a nonsingular transformation be given by Qn Q12 .9.2 I Q22 T L — J (B.15) and Qn Q12 2 Q21 Q22 0 0 0 1 (B.16) to be complete, where I i s an ( i x l ) i d e n t i t y matrix. The s u b s t i t u t i o n of th i s transformation i n the quadratic form Q i n (B.14) (af t e r some algebraic manipulations) y i e l d s ; T <u _T ~T T p > 0 0 -X n T T 2 Ql1 « Ql2 0 u f p (B.17) where X and X are (mxm) and (nxn) diagonal matrices with p o s i t i v e e n t r i e s , r e s p e c t i v e l y . The f u n c t i o n a l form i n (B.14) for transformed variables then leads to the following matrix equation: SA . m 0 Q n a -A Q 1 2a T T a Qi1 a Q 1 2 0 = 0. (B.18) 217 The matrix S of the coefficients in (B.18) and S of the coefficients in (B. l l ) , both have the same number of positive and negative eigenvalues. This is because of the congruence transformation and the Law of Inertia. Therefore i t is only required to find the positive and negative eigenvalues for S. T T Let Qua=6 and Qjia=n, then the matrix S can be written as A? 0 0 A ™ 0 symmetric 0 0 0 0 (m*n) 1 L _ -X1 0 0 -An 0 -A-5 1 1 6 12 ( S21 ( 522 5ml 6m2 oTnnni2 0 jH21H22 0 n n l r l n 2 I I '11 >21 ml " I I 121 n l 0 ( i x l ) (B.19) The principal minors for the matrix (B.19) can be written as M i =UiI; M 2 = AT 0 , m A? 0 0 0 Af 0 etc. ,0 0 A m Since A™ and A^ are a l l positive real numbers, hence a l l the principle minors up to Mm are positive. The (m+l) t h principle minor is simply 218 -X1^!^. Therefore i t is negative and the sign alternates thereafter up to (m+n)1"*1. Beyond (m+n), the principal minors can be either positive, zero, or negative and for rank r=mfn+1, the zero entry would only appear between a positive and a negative entry in the i n d i c i a l sequence <l,Mj,M2, . . ., M ,M ._, . . .,M , ,M , , l S . . .,M , , >. Clearly the index I is m and m m+1 n+m n+m+1 m+n+1 hence only m eigenvalues of the matrix S are positive while (n+l) eigen-values are negative, which also holds for matrix S in (B.13). 219 APPENDIX C EQUIVALENCE OF ENERGY PRODUCTS FOR FOUR FIRST ORDER BEAM EQUATIONS AND PLANE LINEAR ELASTICITY EQUATIONS The energy products for the four f i r s t order beam equations (7.18) and the plane stress linear elasticity equations (5.19) after introducing the simple beam theory assumptions are shown to be the same when the shear energy contribution is dropped. The four f i r s t order beam equations 0 0 0 0 0 D _D -1 1_ ' E I 0 where D^J-; lead to the following energy product; ) V q e 0 ) M 0 V i_ J 0 ( C l ) (AlA,A) = ! [ _ v ^ - e | i - 2 v e + l ^ - ^ f ^ ] d x (C.2) - - - . J1 dx dx dx EI dx From equation (5.11), the energy product for plane stress linear elasticity (unit thickness) is (A2A,A) = / [u TT*T+T TTu-T TCx ] d n (C.3) ,d6 M2, „dv, for A = <u x > and where X* = -T = _9 A = 0 T* T -C - ~ 0 8_ "3y '3y 3_ '3x -v 0 -v 0 1 0 0 2(l+ v) 220 T T and u=<u v> ; 1~<XXX Tyy Txy > ' When t n e basic assumptions of plane sections remain plane after bending and small deflections, following the nomenclature of Figure 7, the following quantities are obtained: 0 (C.4a) yy x = _ Itt XX I V ,h2 2, T x y = 2 i ( r - y2> (C.4b) (C.4c) and u = -y9. (C.4d) He re I is the moment of inertia about the z—axis; h, the height of the beam, and M,V,0 and v are functions of x only. Since x =0: the matrix yy operators T* and T and the matrix C reduce to 3_ "9x 9_ "9y 3_ 3x (C4e) 5-1 1 0 0 2(l+v) (C.4f) Now the substitution of equations (C.4) into (C.3) yields the following. The f i r s t term on the right hand side of (C.3) becomes Ja u T*xdft = ji <-y6 v> 3_ _3_ 3x ~3y 0 3x My_ " i V ,h2 2 . dydx f f2 / Z 2 « d M y.2flw v d V ^ h 2 2SXA A j l hi I dx I 6 V-2I d^ (4- y ) } d y d x -2 h r2 Since for unit thickness, y2dy=I; therefore integration on y yields; "2 221 r uTT*T<m = ! (-efLev-^) tt - - - ' 1 dx dx (C.5) Now the second term on the right hand side of (C.3) takes the form 2 9 3x 0 < -y0 3 3 V _3y 3x dydx = h fh {i " s - i i ^ - y ) 42i (4 ~ - y )d^ } d y d x -2 Again integration on y gives T, L T XTufi = /, (M^-6V+V^)dx 1 tt - — J1 dx dx Finally 2 ' M v V (^-y 2)> I 2I V4 1 . 0 0 2(l+v) f My_ ~I dydx V ,h2 2 I M <rr-yz> f / l /I ^ + g 2 ( | - 2 - y 2 ) 2 > d y d x . r2 ry2., 9.h 2 ,h2 h. 2 This, after integration on y, yields a Adding the equations (C.5), (C.6) and (C.7) gives, ( C 6 ) (C.7) /. [u TT*x+T TTu-T TCT]dfi = f {-vf-0^-2V0+M^-^fV^+^4 2V 2]dx (C.S) tt - - - - — - — ' l dx dx dx EI dx 5EI where j1 5 E I h 2V 2dx is the contribution to the energy product from the shear stress due to flexure. This is only significant for short beams and when i t is dropped, the energy product in (C.3) becomes 222 <fei.A> - / 2 [ - v i - ^ f - ^ - E - I ^ ^ - ' «=•' This is exactly the same as in equation (C.2). 223 APPENDIX D ELASTICITY SOLUTIONS FOR A CANTILEVER An elas t i c i t y solution for the cantilever with the boundary tractions and the boundary conditions shown in Figure 36 is derived. It is also shown that the strain energy computed from the normal stress T and the shear stress x distributions is actually equal to half the xx xy J n total work done by the boundary tractions in moving through the displace-ments obtained from the elasticity solution on the boundaries. The stress distributions are taken as x ( D . l a ) XX I x = 0 (D.lb) yy Txy " - h ( c 2 ~ ? 2 ) ( D - 1 C ) where P is the total load due to the shear stress at the ends. The equations (D.l) identically satisfy the equilibrium Vl+(j)=0 for T x x = g ^ 2 > x = % and x =--—j— everywhere inside the cantilever, as well as yield yy 3x z xy . 3x3y J ' 3 the same stresses on the boundary as applied tractions. The corresponding strains, by applying Hooke's law, are 3u Txx Pxy , „ . £ = T:— = ~Tr~ = - -^r~ (D.2a) xx 3x E EI 3v xx vPxy . N eyy = ^ = ~Y- = -ET (D'2B) Y = Iii + = I S = _ i ( c2 2) ( D > 2 c ) Yxy 3y 3x G 2IG ^ y ; ^.u.zcj E where ^ =2(l+v)' ^ e s t r a ^ - n s derived in (D.2) also satisfy identically the 224 82e 32e 8 2 Y compatibility equation ( 9 y 2 + ^J? = dxd*J) • By integration of (D.2a) and (D.2b), u and v are obtained as u = - f f ^ + f i ( y ) ; v = ^ f f ^ + f 2 ( x ) (D.3) in which the functions f\ and f 2 are as yet unknown functions of y and x only, respectively. Substitution of u and v above into equation (D.2c) yields _ P x j dfl(y) + v P Z j + d f 2 ( x ) =__J_ 2 , 2EI dy 2EI dx 2IG ^ Y J K ' In equation (D.4) some terms are functions of x only, some are functions of y only, and one is independent of both x and y. These can be grouped as w x ) = _ P x j + df 2(x) * W 2EI + dx g ( y ) = y l l l + dfi(y) _ Pzi 2EI dy 2IG Pc 2 K = - (constant). Thus (D.3) becomes F(x) + G(y) = K. (D.5) It can be concluded from (D.5) that F(x) and G(y) must be constants. Denoting these by d and e, respectively, therefore d + e = k a n d df 2(x) = Pxj dfi(y) = _ vPyj. P y j . + a n d dx 2EI + d ' dy 2EI + 2IG 6" Integrating these yields the functions fi(y) and f 2 ( x ) ; Px^ f 2(x) = + dx + h. 225 Substitution in the expression for u and v in (D.3) gives Px 2y vPy3 _,_ Py 3 , . The four constants d,e,g and h can now be determined from the four boundary conditions (Figure 36) u(L,0) = v(L,0) = u(L,c) = u(L,-c) = 0 (D.8) and these are found to be d - - m [ 3 + ( 4 + 5 v ) i7 ] g = 0 and h = ||1 [1 + f(4 +5v) . Finally the substitution for d,e,g and h in the equations (D.6) and (D.7) and letting E,=j- and r\=2- gives J_i c u(5,n) = H^[^n(l-C2)-(^)^r1(l-n2)] (D.9) P T 3 T 1 ^ 1 r 2 v(c,n) = H I[i-fc:4ic- 34i^ 2-{3vc:n 2+(4+5v)(i-c:)}]. (D.IO) Therefore the tip deflection 6 at C (Figure 36) is given by PT 3 1 r - 2 6 = v(0,0) = [14j(4+5v)^2-] ( D . l l ) while the longitudinal deflection u at B, where 5=0 and n=-l, is obtained as UB " "V'-V = - l i T - <D-12> The strain energy in the cantilever can be computed by using the assumed stresses in (D.l) and the corresponding strains in (D.2) in the usual manner. 226 U = Tr /k f C [T e +T e +T Y ldydx 2 1 o -c xx xx yy yy xy xy After i n t e g r a t i o n , the s t r a i n energy i n the can t i l e v e r i s p2j 3 i o „2 U = -ggY" U+;p(l+v)^r]. (D.13) The work done on the boundary can be computed from the following l i n e i n t e g r a l W = f (x..n.u.)ds (D.14) where i s the portion of the boundary where the stresses are prescribed and n^ . , the unit outward normal. For the c a n t i l e v e r i n Figure 36, this i n t e g r a l takes the form W = C T x / ° > y > v ( ° > y ) d y + / ! C T x y ( L , y ) v ( L , y ) d y + jC_c T x x ( L , y ) u ( L , y ) d y = Wx + W2 + W3. (D.15) Here T x y ( ° ' y ) = " f l ( c 2~y 2> v(0,y) [ I 4 | ( 4 + 5 v ) ^ ] . Thus the f i r s t i n t e g r a l i n the r i g h t hand side of (D.15) i s given by 2 3 2 wi = C - 3EI5" ( c 2 _ y 2 ) ti+j<4+5 )^i2-]dy and a f t e r i n t e g r a t i o n W i = l i i ~ t i 4 i ( 4 + 5 v ) i > ] - (D-16) Next, i x ^ ( L , y ) i s the same as x (0,y) and v(L,y) i s given by , T N _ v PL 2 v(L,y) - "J g l y ' Then the second i n t e g r a l W2 i s 227 W2 = j _ c ~ 4 ^ f l y z ( c ^ - y 2 ) d y and leads to F i n a l l y to obtain W3 p 2 P 2 j W 2 = - i o — • <D-17> and u(L,y) = | * ( y 2 - c 2 ) Therefore which gives 3 J-c — 6 — ET2" y ( y c ) d y _ (2+v) P 2 c 2 L W3 - 1 5 E I . (D.18) Adding (D.16), (D.17) and (D.18) y i e l d s P T 3 i i r2 W = Wj + W2 + W3 = [14y^(l+v)^2-]. (D.19) Comparison of (D.19) with (D.13) gives W = 2U. (D.20) Therefore the s t r a i n energy computed from the stress d i s t r i b u t i o n s i n (D.l) i s exactly h a l f the work done by the boundary tractions i n going through boundary displacements as obtained from the e l a s t i c i t y solutions (D.9) and (D.10).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Convergence of mixed methods in continuum mechanics...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Convergence of mixed methods in continuum mechanics and finite element analysis Mirza, Farooque Aguil 1977
pdf
Page Metadata
Item Metadata
Title | Convergence of mixed methods in continuum mechanics and finite element analysis |
Creator |
Mirza, Farooque Aguil |
Date Issued | 1977 |
Description | The energy convergence of mixed methods of approximate analysis for problems involving linear self-adjoint operators is investigated. A new energy product and the associated energy norm are defined for such indefinite systems and then used in establishing the strain energy convergence and estimation of error for problems in continuum mechanics. In the process, the completeness requirements are laid out for approximate solutions. Also established is the mean convergence of the basic variable e.g. displacements and stresses. After accomplishing a new mathematical framework for the mixed methods in continuum, the theory is then extended to the finite element method. The completeness requirements, convergence criteria and the effect of continuity requirements on convergence are established. The flexibility offered by the mixed methods in incorporating the boundary con ditions is also demonstrated. For stress singular problems, the strain energy convergence is established and an energy release method for determining the crack intensity factor K. is presented. A detailed eigenvalue-eigenvector analysis of the mixed finite element matrix is carried out for various combinations of interpolations for the plane stress linear elasticity and the linear part of the Navier-Stokes equations. Also discussed is its relation to the completeness requirements. Finally, numerical results are obtained from applying the mixed finite element method to several examples. These include beam bending, a plane stress square plate with parabolically varying end loads, a plane stress cantilever and plane strain stress concentration around a circular hole. A plane stress example of a square plate with symmetric edge cracks is also solved to study the strain energy convergence. Lastly, two rectangular plates, one with symmetric edge cracks and the other with a central crack are considered to determine the crack intensity factor K. In most of the examples, the strain energy convergence rates are predicted and compared with the numerical results, and excellent agreement is observed. |
Subject |
Finite element method Continuum mechanics |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-21 |
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.0062681 |
URI | http://hdl.handle.net/2429/20637 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1977_A1 M57.pdf [ 11.4MB ]
- Metadata
- JSON: 831-1.0062681.json
- JSON-LD: 831-1.0062681-ld.json
- RDF/XML (Pretty): 831-1.0062681-rdf.xml
- RDF/JSON: 831-1.0062681-rdf.json
- Turtle: 831-1.0062681-turtle.txt
- N-Triples: 831-1.0062681-rdf-ntriples.txt
- Original Record: 831-1.0062681-source.json
- Full Text
- 831-1.0062681-fulltext.txt
- Citation
- 831-1.0062681.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062681/manifest