meee ut Third Edition AN INTRODUCTION TO THEAN INTRODUCTION TO THE FINITE ELEMENT METHODMcGraw-Hill Series in Mechanical Engineering Alciatore and Histand: Introduction to Mechatronics and Measurement Systems Anderson: Computational Fluid Dynamics: The Basies with Applications Anderson: Fundamentals of Aerodynamics Anderson: Introduction to Flight Anderson: Modern Compressible Flow Barber: Intermediate Mechanics of Materials Beer/Johnston: Vector Mechanics for Engineers Beer/Johnston/DeWolf: Mechanics of Materials Budynas: Advanced Strength and Applied Stress Analysis Cengel and Boles: Thermodynamics: An Engineering Approach Cengel and Turner: Fundamentals of Thermal-Fluid Sciences Cengel: Heat Transfer: A Practical Approach Cengel: Introduction to Thermodynamics & Heat Transfer Crespo da Silva: Intermediate Dynamics Dieter: Engineering Design: A Materials & Processing Approach Dieter: Mechanical Metallurgy 1: Measurement Systems: Application & Design Measurement & Data Analysis for Engineering & Science EDS, Inc.: /-DEAS Student Guide Hamrock/Schmid/Jacobson: Fundamentals of Machine Elements Heywood: Internal Combustion Engine Fundamentals Holman: Experimental Methods for Engineers Holman: Heat Transfer Hutton: Fundamentals of Finite Element Analysis Kays/Crawford/Weigand: Convective Heat and Mass Transfer Meirovitch: Fundamenials of Vibrations Norton: Design of Machinery Palm: System Dynamics Reddy: An Introduction to the Finite Element Method Schaffer et al.:. The Science and Design of Engineering Materials Schey: Iniroduction to Manufacturing Processes Shames: Mechanics of Fluids Shigley/Mischke/Budynas: Mechanical Engineering Design Smith: Foundations of Materials Science and Engineering Suryanarayana and Arici: Design and Simulation of Thermal Systems Turns: An Introduction to Combustion: Concepts and Applications Ugural: Mechanical Design: An Integrated Approach Ullman: The Mechanical Design Process White: Fluid Mechanics Viscous Fluid Flow Mastering CAD/CAMAn Introduction to the Finite Element Method Third Edition J.N. REDDY Department of Mechanical Engineering Texas A&M University College Station, Texas, USA 77843 wry . . va Higher Education Boston Burr Ridge, IL Dubuque, !A Madison, WI New York San Francisco St. Louis Bangkok Bogota Caracas KualaLumpur Lisbon London Madrid Mexico City Milan Montreal New Delhi Santiago Seoul Singapore Sydney Taipei Torontoey Higher Education AN INTRODUCTION TO THE FINITE ELEMENT METHOD, THIRD EDITION Published by McGraw-Hill, a business unit of The McGraw-Hill Companies, Inc., 1221 Avenue of the Americas, New York, NY 10020. Copyright © 2006, 1993, 1984 by The McGraw-Hill Companies, Inc. All rights served. No part ofthis publication may be reproduced or distributed in any form or by any means, of stored in a database of retrieval system, without the prior written consent of The McGraw-Hill Companies, Inc. including, but not limited to, in any network or other electronic storage or transmission, or broadcast for distance learning. ‘Some ancillaries, including electronic and print components, may not be available to customers, ‘outside the United States. ‘This book is printed on acid-free paper. 1234567890 QPFIQPF0987654 ISBN 0-07-246685-5 Senior Sponsoring Editor: Suzanne Jeans Developmental Editor: Amanda J. Green Senior Marketing Manager: Mary K. Kittell Project Manager: Peggy S. Lucas Senior Production Supervisor: Kara Kudronowicz Media Technology Producer: Eric A. Weber ‘Senior Coordinator of Freelance Design: Michelle D. Whitaker Cover Designer: Rebidas Design & Associates Lead Photo Research Coordinator: Carrie K. Burger ‘Supplement Producer: Brenda A. Ernzen Compositor: International Typesetting and Composition ‘Typeface: 10 Computer Modern (cmar 10) Printer: Quebecor World Fairfield, PA (USE) Cover Images: Pinched Cylinder image is generated by the data included in the paper “Least-squares finite element formulation for shear-deformable shes,” by J. P. Pontaza and J. N. Reddy, Computer Methods in Applied Mechanics and Engineering, 2004. Connecting Rod image courtesy of MSC.Software® Corporation. www.mscsoftware.com Library of Congress Cataloging-in-Publication Data Reddy, J. N. Gunuthula Narasimha), 1945— ‘An introduction to the finite element method / J. N. Reddy —3rd ed, p- cm.—(McGraw-Hill series in mechanical engineering) Includes bibliographical references and index. ISBN 0-07-246685-5 1. Finite element method. I. Title. Il. Series. TA47.FSR4 2006 620,001'51825—de22 2004058177 ‘www.mhhe.comABOUT THE AUTHOR J.N. Reddy is a Distinguished Professor and the inaugural holder of the Oscar S. Wyatt Endowed Chair in the Department of Mechanical Engineering at Texas A&M University, College Station, Texas. Prior to his current position, he worked as a postdoctoral fellow at the University of Texas at Austin, as a research scientist for Lockheed Missiles and Space Company, and taught at the University of Oklahoma and Virginia Polytechnic Institute and State University, where he was the inaugural holder of the Clifton C. Garvin Endowed Professorship. Professor Reddy is the author of over 320 journal papers and 14 text books on theoretical formulations and finite element analysis of problems in solid and structural mechanics (plates and shells), composite materials, computational fluid dynamics and heat transfer, and applied mathematics. Professor Reddy is the recipient of numerous awards including the 1984 Walter L. Huber Civil Engineering Research Prize of the American Society of Civil Engineers, the 1985 Alumni Research Award at Virginia Polytechnic Institute, the 1992 Worcester Reed Warner Medal and the 1995 Charles Russ Richards Memorial Award of the American Society of Mechanical Engineers, the 1997 Melvin R. Lohmann Médal from Oklahoma State Univer- sity, the 1997 Archie Higdon Distinguished Educator Award from the Mechanics Division of the American Society of Engineering Education, the 1998 Nathan M. Newmark Medal from the American Society of Civil Engineers, the 2000 Excellence in the Field of Com- posites Award and the 2004 ASC Outstanding Research Award from the American Society of Composite Materials, the 2000 Faculty Distinguished Achievement Award for Research, and the 2003 Bush Excellence Award for Faculty in International Research Award from Texas A&M University, and the 2003 Computational Structural Mechanics Award from the U.S. Association for Computational Mechanics. He is a fellow of the American Academy of, Mechanics, the American Society of Civil Engineers, the American Society of Mechanical Engineers, the American Society of Composites, the International Association of Computa- tional Mechanics, the U.S. Association of Computational Mechanics, and the Aeronautical Society of India. Dr. Reddy is the Editor-in-Chief of the journals Mechanics of Advanced Materials and Structures, International Journal of Computational Methods in Engineering Science and Mechanics and International Journal of Structural Stability and Dynamics; he also serves on the editorial boards of over two dozen other journals.To My teachers and studentsCONTENTS Preface Introduction 1.1 General Comments 1.2. Mathematical Models 1.3. Numerical Simulations 1.4 The Finite Element Method 14.1 The Basic Idea 1.4.2 The Basic Features 1.4.3 Some Remarks 1.44 A Brief Review of History and Recent Developments 1.5 The Present Study 1.6 Summary Problems References for Additional Reading Mathematical Preliminaries, Integral Formulations, and Variational Methods 2.1 General Introduction 2.1.1 Variational Principles and Methods 2.1.2 Variational Formulations 2.1.3 Need for Weighted-Integral Statements 2.2 Some Mathematical Concepts and Formulae 2.2.1 Coordinate Systems and the Del Operator 2.2.2 Boundary Value, Initial Value, and Eigenvalue Problems 2.2.3 Integral Identities 2.2.4 Linear and Bilinear Functionals 2.3. Elements of Calculus of Variations 2.3.1 Introduction Variational Operator and First Variation Fundamental Lemma of Variational Calculus ‘The Buler Equations Natural and Essential Boundary Conditions Hamilton’s Principle 27 27 27 28 28 31 31 33 36 39 41 4L 41 47 54viii ‘CONTENTS. 24 Integral Formulations 2.4.1 Introduction 2.4.2. Weighted-Integral and Weak Formulations 2.4.3 Linear and Bilinear Forms and Quadratic Functionals 24.4 Examples 2.5 Variational Methods 2.5.1 Introduction 2.5.2 The Ritz Method 2.5.3 Approximation Functions 2.5.4 Examples 2.5.5. The Method of Weighted Residuals 2.6 Summary Problems References for Additional Reading Second-Order Differential Equations in One Dimension: Finite Element Models 3.1 Background 3.2 Basic Steps of Finite Element Analysis 3.2.1 Model Boundary Value Problem 3.2.2 Discretization of the Domain 3.2.3 Derivation of Element Equations 3.2.4 Connectivity of Elements 3.2.5 Imposition of Boundary Conditions 3.2.6 Solution of Equations 3.2.7 Postcomputation of the Solution 3.3. Some Remarks 3.4 Axisymmetric Problems 3.4.1, Model Equation 3.4.2 Weak Form 3.4.3 Finite Element Model 3.5 Summary Problems References for Additional Reading Second-Order Differential Equations in One Dimension: Applications 4.1 Preliminary Comments 4.2. Discrete Systems 4.2.1 Linear Elastic Spring 4.2.2 Torsion of Circular Shafts 4.2.3 Electrical Resistor Circuits 4.2.4 Fluid Flow through Pipes 4.3 Heat Transfer 4.3.1 Governing Equations 4.3.2 Finite Element Models 4.3.3. Numerical Examples 358 58 58 64 66 74 74 74 16 n 1 97 98 102 103 103 105 105 106 108 125 132 132 134 141 146 146 147 148 150 151 154 155 155 156 156 158 159 161 162 162 166 1664.4 Fluid Mechanics 4.4.1 Governing Equations 4.4.2 Finite Element Model 4.5. Solid and Structural Mechanics 4.5.1 Preliminary Comments 45.2. Finite Element Model of Bars and Cables 45.3 Numerical Examples 4.6 Plane Trusses 4.6.1 Introduction 4.6.2 Basic Truss Element 4.6.3. General Truss Element 4.6.4 Constraint Equations: Penalty Approach 4.6.5 Constraint Equations: A Direct Approach 4.7 Summary Problems References for Additional Reading Beams and Frames 5.1 Introduction 5.2. Euler-Bernoulli Beam Element 5.2.1 Governing Equation 5.2.2 Discretization of the Domain 5.2.3 Derivation of Element Equations 5.2.4 Assembly of Element Equations 5.2.5 Imposition of Boundary Conditions 5.2.6 Postprocessing of the Solution 5.2.7 Numerical Examples 5.3. Timoshenko Beam Elements 1 Governing Equations 5.3.2. Weak Form 5.3.3, General Finite Element Model 5.3.4. Consistent Interpolation Elements 5.3.5. Reduced Integration Element 5.3.6 Numerical Examples 5.4 Plane Frame Elements 5.4.1 Introductory Comments 5.4.2 Frame Element 5.5 Summary Problems References for Additional Reading Eigenvalue and Time-Dependent Problems 6.1 Eigenvalue Problems 6.1.1 Introduction 6.1.2 Formulation of Eigenvalue Problems 6.1.3 Finite Element Formulation 6.2 Time-Dependent Problems 6.2.1 Introduction 6.2.2 Semidiscrete Finite Element Models contents ix 181 181 181 183 183 184 185 194 194 194 195 202 2u1 214 215 231 233 233 233 233 234 234 243 245 247 248 261 261 262 264 266 270 m 274 274 274 281 282 290 291 291 291 292 295 314 314 316Parabolic Equations Hyperbolic Equations 5 Mass Lumping 6.2.6 Applications 6.3 Summary Problems References for Additional Reading Computer Implementation 7. Numerical Integration 7.1.1 Background 7.1.2. Natural Coordinates 7.1.3 Approximation of Geometry 7.1.4 Isoparametric Formulations 7.1.5 Numerical Integration 7.2. Computer Implementation 7.2.1 Introductory Comments 7.2.2. General Outline 1.23. Preprocessor 7.24 Calculation of Element Matrices (Processor) 7.2.5 Assembly of Element Equations (Processor) 7.2.6 Imposition of Boundary Conditions (Processor) 7.2.7 Solving Equations and Postprocessing 7.3. Applications of Program FEMID 7.3.1 General Comments 7.3.2 Mlustrative Examples 74 Summary Problems References for Additional Reading Single-Variable Problems in Two Dimensions 8.1 Introduction 8.2 Boundary Value Problems 8.2.1 The Model Equation 8.2.2 Finite Element Discretization 8.2.3. Weak Form 8.2.4 Finite Element Model 8.2.5 Derivation of Interpolation Functions 8.2.6 Evaluation of Element Matrices and Vectors 8.2.7 Assembly of Blement Equations 8.2.8 Postcomputations 8.2.9 Axisymmetric Problems 8.3 A Numerical Example 8.4 Some Comments on Mesh Generation and Imposition of Boundary Conditions 8.4.1 Discretization of a Domain 8.4.2 Generation of Finite Element Data 8.4.3 Imposition of Boundary Conditions 318 324 326 328 337 337 342 343 343 343 345 346 347 348, 356 356 397 359 360 363 365 367 370 370 370 401 401 406 409 410 410 4u 412 415 417 425 436 44i 442 453 453 455 456| covets xi 8.5. Applications 458 8.5.1 Conduction and Convection Heat Transfer 458 8.5.2. Fluid Mechanics 472 8.5.3 Solid Mechanics 485 8.6 Eigenvalue and Time-Dependent Problems 490 8.6.1 Introduction 490 8.6.2. Parabolic Equations 491 8.6.3 Hyperbolic Equations 499 8.7 Summary 504 Problems 504 References for Additional Reading 522 9 Interpolation Functions, Numerical Integration, and Modeling Considerations 525 9.1 Introduction 525 | 9.2. Blement Library 525 9.2.1 Triangular Elements 525 9.2.2 Rectangular Elements $32 9.2.3. The Serendipity Elements 537 9.2.4 Hermite Cubic Interpolation Functions 539 9.3 Numerical Integration 540 9.3.1 Preliminary Comments 540 9.3.2 Coordinate Transformations 543 9.3.3 Integration over a Master Rectangular Element 549 9.3.4 Integration over a Master Triangular Element 557 9.4 Modeling Considerations 361 9.4.1 Preliminary Comments 561 9.4.2 Element Geometries 562 9.4.3 Mesh Generation 563 9.4.4 Load Representation 567 9.5 Summary 569 ' Problems 570 | References for Additional Reading 575 10 Flows of Viscous Incompressible Fluids 577 | 10.1 Preliminary Comments 317 | 10.2 Governing Equations 577 10.3 Velocity-Pressure Formulation 319 10.3.1 Weak Formulation 579 10.3.2. Finite Element Model 381 10.4 Penalty Function Formulation 583 | 10.4.1 Preliminary Comments 583 i 10.4.2. Formulation of the Flow Problem as a Constrained Problem 583 10.4.3 Lagrange Multiplier Model 584 | 10.4.4 Penalty Modet 585 | 10.4.5. Time Approximation 588 | eexii it 12 ‘CONTENTS. 10.5 106 107 ‘Computational Aspects 10.5.1 Properties of the Matrix Equations 10.5.2. Choice of Elements 10.5.3 Evaluation of Element Matrices in the Penalty Model 10.5.4 Postcomputation of Stresses ‘Numerical Examples ‘Summary Problems References for Additional Reading Plane Elasticity Wal 12 113 4 is 11.6 7 18 Introduction Governing Equations 11.2.1 Plane Strain 11.2.2 Plane Stress 11.2.3 Summary of Equations Weak Formulations 113.1 Preliminary Comments 11.3.2 Principle of Virtual Displacements in Vector Form 11.3.3 Weak Form of the Governing Differential Equations nite Element Model 11.4.1 General Model 11.4.2 Eigenvalue and Transient Problems Evaluation of Integrals Assembly of Finite Element Equations Examples Summary Problems References for Additional Reading Bending of Elastic Plates 124 12.2 123 124 125 12.6 Introduction Classical Plate Theory 12.2.1 Displacement Field 12.2.2 Virtual Work Statement 12.2.3 Finite Element Model 12.2.4 Plate Bending Elements Shear Deformation Plate Theory 12.3.1 Displacement Field 12.3.2 Virtual Work Statement 12.3.3 Finite Element Model 12.3.4 Shear Locking and Reduced Integration Eigenvalue and Time-Dependent Problems Examples Summary Problems References for Additional Reading 588 588 589 590 591 591 602 603 605 607 07 607 607 608 610 612 612 612 613 614 614 617 617 620 622 629 629 633 635 635 637 637 638 642 643 646 646 648 650 652 653 655 663 663 66513 Computer Implementation of Two-Dimensional Problems 14 1B. 13.2 13.3 134 13.5 Introduction Preprocessor Element Computations (Processor) Applications of the Computer Program FEM2D 13.4.1 Introduction 13.4.2 Description of Mesh Generators 13.4.3 Applications (Illustrative Examples) Summary Problems References for Additional Reading Prelude to Advanced Topics 14.1 14.2 143 144 145 146 Introduction Alternative Finite Element Models 14.2.1 Introductory Comments 14.2.2 Weighted Residual Finite Element Models 14.2.3 Mixed Formulations ‘Three-Dimensional Problems 14.3.1 Heat Transfer 14.3.2 Flows of Viscous Incompressible Fluids 143.3 Elasticity 143.4 Three-Dimensional Finite Elements 14.3.5 A Numerical Example Nonlinear Problems 14.4.1 General Comments 14.4.2 Bending of Euler-Bernoulli Beams 14.4.3. The Navier-Stokes Equations in Two Dimensions 14.4.4 Solution Methods for Nonlinear Algebraic Equations 144.5 Numerical Examples Errors in Finite Element Analysis 14.5.1 Types of Errors 14.5.2 Measures of Errors 14.5.3 Convergence and Accuracy of Solutions ‘Summary Problems References for Additional Reading Index 667 667 669 669 675 675, 681 686 703 705 709 mL nm m1 1 12 m2 725 726 1 TB 731 735 736 736 736 738 39 740 743 a3 744 745 750 751 153 757PREFACE The third edition of the book, like the previous two editions, represents an effort to select and present certain aspects of the finite element miethod that are most useful in developing and analyzing linear problems of engineering and science. In revising the book, students taking courses that might use this book as their textbook have been kept in mind. This edition is prepared to bring more clarity to the concepts being discussed while maintaining the necessary mathematical rigor and providing physical interpretations and engineering applications at every step. ‘This edition of the book is a revision of the second edition, which was well received by engineers as well as researchers in the fields of engineering and science. Most of the revisions in the current edition took place in Chapters | through 6 and 14. Chapter 5 on error analysis from the second edition is eliminated and a section on the same subject is added to Chapter 14. Chapter 3 from the second edition is now divided into two chapters, Chapter 3 on theoretical formulation and Chapter 4 on applications. Chapter 4 on beams from the second edition now became Chapter 5 in the current edition, making the total number of chapters in both editions the same. Another change is the interchanging of Chapters 10 and 11 to facilitate the natural transition from plane elasticity to plate bending (Chapter 12). New material on three-dimensional finite element formulations and nonlinear formulations is added in Chapter 14. In all the chapters, material is added and reorganized to aid the reader in understanding the concepts. Another change in this edition is the removal of appendices containing the program listing of FEMLD and FEM2D, which would have taken up over 80 pages. The pro- gram listings, executables, and instructor resources are available on the text’s website at www.mhhe.com/reddy3e. The source files of these programs may be obtained from the author for a small fee. Most people who have used the earlier editions of the book liked the “differential equa- tions approach” adopted here. This is natural because everyone, engineer or scientist, is looking for a way to solve differential equations arising in the study of physical phenom- ena. It is hoped that the third edition of the book serves the student even better than the previous editions in understanding the finite element method as applied to linear problems of engineering and science. The author has benefited by teaching an introductory course from this book for many years. While it is not possible to name hundreds of students and colleagues who have contributed to the author’s ability to explain concepts in a clear manner, the author xv aXVi PREFACE expresses his sincere appreciation to the following people for their helpful comments on this edition: Ronald C. Averill, Michigan State University Mahesh Gupta, Michigan Technological University HUN. Hashemi, Northeastern University Stephen M. Heinrich, Marquette University Paul R. Heyliger, Colorado State University John Jackson, University of Memphis Mehdi Pourazady, University of Toledo Robert L. Rankin, Arizona State University N. Sukumar, University of California, Davis Kumar Tamma, University of Minnesota Chi-Tay Tsai, Florida Atlantic University JN. Reddy jnreddy @tamu.edu Tejashwina vadheetamasthu (May what we study be well studied)Chapter 1 INTRODUCTION 1.1 GENERAL COMMENTS, One of the most important things engineers and scientists do is to model physical phe- nomena. Virtually every phenomenon in nature, whether aerospace, biological, chemical, geological, or mechanical can be described, with the aid of the laws of physics or other fields in terms of algebraic, differential, and/or integral equations relating various quantities of interest, Determining the stress distribution in a pressure vessel with oddly shaped holes and numerous stiffeners and subjected to mechanical, thermal, and/or aerodynamic loads, finding the concentration of pollutants in lakes and seawater or in the atmosphere, and simulating weather in an attempt to understand and predict the formation of thunderstorms and tornadoes are a few examples of many important practical problems that engineers deal with. Analytical descriptions of physical phenomena and processes are called mathematical ‘models. Mathematical models of a process are developed using assumptions concerning how the process works and using appropriate axioms or laws governing the process, and they are often characterized by very complex differential and/or integral equations posed on geometrically complicated domains. Consequently, the processes to be studied, until the advent of electronic computation, were drastically simplified so that the governing equations can be solved analytically. Over the last three decades, however, computers have made it possible, with the help of suitable mathematical models and numerical methods, to solve many practical problems of engineering. The use of a numerical method and a computer to evaluate the mathematical model of a process and estimate its characteristics is called a numerical simulation. There now exists a new and growing body of knowledge connected with the development of mathematical models and use of numerical simulations of physical systems, and it is known as computational mechanics {1}. Any numerical simulation, such as the one by the finite element method, is not an end in itself but rather an aid to design and manufacturing. There are several reasons why an engineer or a scientist should study a numerical method, especially the finite element method: 1. Most practical problems involve complicated domains (both geometry and material constitution), loads, and nonlinearities that forbid the development of analytical so- lutions. Therefore, the only alternative is to find approximate solutions using numerical methods.2 AN INTRODUCTION TO THE FINITE ELEMENT METHOD 2. A numerical method, with the advent of a computer, can be used to investigate the effects, of various parameters (¢.g., geometry, material parameters, or loads) of the system on its response to gain a better understanding of the process/system being analyzed, It is cost- effective and saves time and material resources compared to the multitude of physical experiments needed to gain the same level of understanding. 3. Because of the power of numerical methods and electronic computation, it is possible to include all relevant features in a mathematical model of a physical process without worrying about its solution by exact means. 4, Those who are quick to use a computer program rather than think about the problem to bbe analyzed may find it difficult to interpret or explain the computer-generated results. Even to develop proper input data to the computer program, a good understanding of the underlying theory of the problem as well as the numerical method (on which the ‘computer program is based) is required. 5. The finite element method and its generalizations are the most powerful computer oriented methods ever devised to analyze practical engineering problems. Today, fi- nite clement analysis is an integral and major component in many fields of engineer- ing design and manufacturing. Major established industries such as the automobile, aerospace, chemical, pharmaceutical, petroleum, electronics, and communications, as well as emerging technologies such as nanotechnology and biotechnology rely on the finite element method to simulate complex phenomenon at different scales for design and manufacture of high-technology products. 1.2 MATHEMATICAL MODELS, ‘A mathematical model can be broadly defined as a set of equations that expresses the essential features of a physical system in terms of variables that describe the system. The ‘mathematical models of physical phenomena are often based on fundamental scientific laws of physics such as the principle of conservation of mass, conservation of linear momentum, and conservation of energy, Below we consider three simple examples drawn from dynam- ies, heat transfer, and solid mechanics to illustrate how mathematical models of physical problems are formulated.‘HAPTER | oTRODUCTION 3 ma = mt me sy Figure 1.2.1 Simple pendulum, which states, in the present case, that the vector sum of externally applied forces on a system is ‘equial to the time rate of change of the linear momentum (mass times velocity) of the system: F < iny= a 1.2.1) where F is the vector sum of all forces acting on the system, 7 is the mass of the system, vis the velocity vector, and a is the acceleration vector of the system, To write the equation governing the angular motion, we set up a coordinate system, as shown in Fig. 1.2.1. Applying Newton's second law to the x-direction (note that the dynamic equilibrium of forces in the y-direction gives the reaction force R(W) in terms of the weight mg of the bob), we obtain where F, =—mg sin, v, = 04, 6 is the angular displacement (radians), vis the component of velocity (nvs) along the x coordinate, und 1 denotes time (s). Thus, the equation for angular motion becomes #0 8 =mg sind +4 sine =0 12.3) 3 he PRE A Equation (1.2.3) is nonlinear on account ofthe term sin @. For small angular motions (consistent with the goal of the study), sind is approximated as @. Thus, the angular motion is described by the linear differential equation (24) Equations (1.2.3) and (1.2.4) represent mathematical models of nonlinear and linear mo- tions, respectively, of a rigid pendulum, Their solution requires knowledge of conditions at time ¢ = 0 on @ and its time derivative 4 (angular velocity). These conditions are known as the initial conditions. Thus, the linear problem involves solving the differential equation (1.2.4) subjected to the initial conditions do. 9(0) =4, 0) = v0 2.5) at ‘The problem described by Egs. (1.2.4) and (1.2.5) is called an initial-value problem. The linear problem described by Eqs. (1.2.4) and (1.2.5) can be solved analytically. The ‘general analytical solution of the linear equation (1.2.4) (@ + 476 =0) is given by ow. A sind + B cos Ar (1.2.6)4 [AN INTRODUCTION 10 THE FINITE ELEMENT METHO! where 4=,/% und A and B are constants to be determined using the initial conditions in (1.2.5), We obtain + Bab 22 and the solution to the linear problem is Oc = FE sin At +A 08 24 1.2.8) For zero initial velocity and nonzero initial position 4, we have (1) =A cos dt (129) which represents a simple harmonic motion. If we were to solve the nonlinear equation (1.2.3) subject to the conditions in (1.2.5), we may consider using a numerical method because it is not possible to solve it exactly for large values of @. We will revisit this issue in the sequel. Example 1.2.2 (A Heat Transfer Problem) Here we wish to derive the governing equations (i., develop the mathematical model) of steady-state heat transfer through a cylindrical bar of nonuniform cross section. The bar is subject to a known temperature Tp (°C) at the left end and exposed, both on the surface and at the right end, to a medium (such as cooling fluid oF air) at temperature Ta. We assume that temperatute is uniform at any section of the bar, 7 = T(x). Due to the difference between the temperatures of the bar and the surrounding medium, there is convective heat transfer across the surface of the body and at the right end. The principle of conservation of energy (or the second law of thermodynamics) can be used to derive the governing equations of the problem. The principle of conservation of energy requires that the rate of change (increase) ‘of internal energy is equal to the sum of heat gained by conduction, convection, and internal heat generation (radiation not included). For a steady process, the time rate of internal energy is zero, Consider a volume element of length Ax and having an area of cross section A(x) (m?) normal to the x axis (see Fig. 1.2.2). If q denotes the heat flux heat flow per unit area, ‘Wm2), then [Aq], is the net heat flow into the volume element at x, [Aq]..ax is the net heat flow out of the volume element atx + Ax, and f PAx(Ts, ~ 7) is the heat flow through the surface of the rod into the body. Here f denotes the film (that is formed between the body and the medium around) conductance (WAm?-°C)], Tx. is the temperature of the sur- rounding medium, and P is the perimeter (m). We also assume that there is a heat source within the rod generating energy at a rate of g (W/m"). In practice, such energy source can be due to nuclear fission or chemical reactions taking place within the rod, or due to the pas- volume heating) . Then the energy balance sage of electric current through the medium ( gives i, [Ag], ~ [Aqhsyax + BP Ax(Too —T) +94 Ax (1.2.10)Convection from lateral surface (0, internal heat generation Maintained at temperature, Ty > Exposed to ambient temperature, 7 nar IS ast Figure 1.2.2 Heut transfer in a bar on, dividing throughout by Ax, [Agleax ~ LA], Ox + AP(Tx —T) + Ag=0 and taking the limit Ax > 0, we obtain Ag) + BP(Tx-T)+ Ag=0 2.4) We can relate the heat flux q (W/m?) to the temperature gradient, Such a relation is provided by the Fourier law (1.2.12) where k denotes the thermal conductivity [WAm °C)] of the material. The minus sign an the right side of the equality in the above equation indicates that heat flows from high temperature to low temperature, Equation (1.2.12) is known as a constitutive relation because it contains & material property Now using the Fourier law (1.2.12) in (1.2.11), we arrive at the heat conduetion equation S (ta which can also be written as 0 (1.2.13) ar A (tall) varer—ty=he ; (nonhomogeneous) second-order differential equation, which are Tor heat Aq (not both) at a boundary Equation (1.2.14) is alin can be solved with known conditions on tempera 5point (¢.g., ends of the bar). The known end conditions for the present case (recall the problem description at the start of this example) can be expressed as TO)=To, [eaSZ spac -r)| 0 (1.2.15) These conditions are called boundary conditions because they represent conditions at the boundary points of the bar. The first condition in (1.2.15) is obvious. The second condi- tion represents the balance of heat due to conduction [A(@7/d.x)] and convection [BA(T T,.)]. The problem described by Eqs. (1.2.14) and (1.2.15) is called, for obvious reasons, a boundary-value problem. This completes the mathematical model development of the problem Eqeations (1.2.14) and (1.2.15) ean be simplified for various special cases. First let Q=T-Tx, a=ka(Wem'O), cm pP(Wim’O)], fedg (Wm) (1.2.16) so that Eqs. (1.2.14) and (1.2.15) can be writ heat transfer problems) as (note that @ and are always positi (2.17) (1.2.18) Now if a and ¢ are constant (e.g, homogeneous bar of constant cross section) and g = 0 (no internal heat generation), Eqs. (1.2.17) and (1.2.18) become — +£0=0, where O Ay Stress, (AG reas Figure 1.2,3 Axial deformation of a bar.8 tensile force at x + Ax. Then setting the sum of the forces to zero (i second law in the x-direction) yields applying Newton's Aah + (Agr lrrax + 8A Ar =O where g(x) is the body force (1.e., weight of the body). Dividing throughout by Ax a the limit Ax > 0, we obtain id taking d qylde) + Ages) = (1.2.23) which represents the equilibrium of forces in the x-direction. ‘The stress o, can be related to the axial displacement using Hooke's law Gyo Eig, 1.2.24) where £ is Young’s modulus (N/m), u denotes the axial displacement (m), and ¢, is the axial strain (nvm). Again, we see that Eq. (1.2.24) is a constitutive equation, Note that a system can have several constitutive relations, each depending on the phenomenon being studied. The study of heat transfer in a bar (with nonuniform temperature) required us to employ Fourier’s, Jaw to relate temperature gradient to heat flow, and the study of deformation of the same bar (Gubjected to axial forces) requires us to use Hooke’s law to connect stress to displacement gradient Now using Eq, (1.2.24) in iq, (1.2.23), we artive at the equilibrium equation in terms of the displacement dx" dx ‘This second-order equation can be solved using known boundary conditions atx =Qand x = L. ‘The boundary conditions of a bar involve specifying either the displacement u or the force Aa ata houndary point ‘The known boundary condi + (24%) +gA=0, O 0 subject to the (initial) condition (0) = uo. We approximate the derivative at r, by Cc mi) Note that we replaced the derivative at =1, by its definition except that we did not take the limit A¢= 4,41 — 4 > 0; that is why itis an approximation, For increasingly small values of At, we expect the approximation to have decreasingly small error. Substituting (1.3.2) into (3.1) at udhiyn) ~ 6 13.2) we oblain wm FAL SUG), my =U, Bt= hy =H 3.3) 3) can be solved, starting from the known value uo of u(t) att =0,for uy =u (ti) = nethe values of uattimes! = At, 2At,...,.nAt tutta method), also known as Equation (1 (Ar). This process ean be repeated to de This is known as Buler’s explicit scheme (or first-order Rung the forward difference scheme, Of course there are higher-order finite difference schemes that are more accurate than Euler's scheme, and we will not discuss them as they are outside the scope of this study. Note that we are able to convert the ordinary differential equation (1.3.1) to an algebraic equation (1.3.3) that needs to be evaluated at different times to construct the time history of (t) We now apply Buler’s explicit scheme to the second-order equation (1.2.4) subjected to the initial conditions (1.2.5). To apply the procedure described above to the equation at hand, we rewrite Eq. (1.2.4) asa pair of coupled (i.e, one cannot be solved without the other) first-order equations D ie oy, Baie a34) a allwuuiuuuloisteniliu «A 39 —uuluuuula 00 0S 10 15 20 25 30 35 40 Time, t Figure 1.3.1 Analytical and numerical solutions @(2) of the simple pendulum.cuarrer ts wstropuction LI where /£, Applying the scheme of Eq. (1.3.3) to the equations at hand, we obtain BSH ACY vie = — At HO, ‘The expressions for @,,, and v,,1 in Eq. (1.3.5) are repeatedly computed solution (@,, v,) from the previous time step. At time f =0, we use the known initial values (@), 1%). Thus, one needs a computer and a computer language like Fortran or C++ to implement the logic of repeatedly computing 6,41 and vj, with the help of Eq. (1.3.5). ‘The numerical solutions of equation (1.3.5) for three different time steps, Af = 0.05, 0.025, and 0.001, along with the exact linear solution (1.2.9) (with €= 2.0, ¢= 32.2, é =7/4, and ¥» =0) are presented in Fig. 1.3.1. The accuracy of the numerical solutions is dependent on the size of the time step: the solution is more accurate with smaller time steps. For large time steps, the solution even diverges from the true solution. Example 1.3.2 (Heat Transfer Problem) Here we consider the boundary-value problem of Example 1.2.2. The finite difference ap- proximation of a boundary-value problem differs from that of an initial value problem. First we divide the domain (0, L) into a finite set of NV intervals of equal length Ax, as shown in Fig. 1.3.2. The ends of each interval is called a mesh point. Thus, there are N + | mesh points in the domain. Then we approximate the second-derivative directly using the centered difference scheme where the error is of order (Ax)? a ~(2 128+ Bran ae ae), (Ay oe Using the above approximation in Bg. (1.2.19), we obtain G4 = 20+ Oi41) + (AsV =0 or =O +12 + (OMAN)"I6, — ear (13.7) Equation (1.3.7) is valid for any mesh point x =x, =1,2,...,.N, and it contains values of @ both atx =). and x =x),1. Note that Eq, (1.3.7) is not used at mesh point x = x9 =0 ‘because the temperature is known there [see Eq. (1.2,20)]. However, use of Eq, (1.3.7) at mesh point x = xy = L requires the knowledge of the fictitious value. (as we see in the sequel, [sx Ambient temperature, @,= 0 a Aa OE creccecteamentwomnoed x= i rok oe Mesh points a Figure 1.3.2 Heat transfer in u bar and a typical finite difference mesh.12. axtvrropuction ro THe FINITE ELEMENT METHOD: ‘we never have to deall with such fictitious values in the finite element method). The forward finite difference approximation of the second boundary condition in (1.2.20) at mesh point x=L gives Ons = Oy 1 Foy = 0 Ax k from which we have (#28) es, 138 (Ce 036 Application of the formula (1.3.7) to mesh points at xy, 12, +. ay yields 6) + D6, — = 0 6; + Dé, ~ 0 = 0 fs + DPs — Ay = 0 ~By.1 + DBx — O41 = 0 (3.9) where D=(2+(mAx)*]; Eq. (1.3.9) consists of N’ equations in N unknowns, 6 8,211 On. As a specific example, consider a steel rod of diameter d h L = 0.05 m, and thermal conductivity k= 50 WAm=°C). Suppose that the temperature at the left end is Ty = 320°C, ambient temperature is T, = 20°C and film conductance (or heat transfer coeffi- ).02 m, len; cient) B = 100 Witm?-"C). For this data, we have B >_ BP _ Bird) £2 om = EEO = 400, =0(0) =300°C Kk m~ kA k(rd?/4) = For a subdivision of four intervals (W =4), we have Ax=0.0125 m and D=2+ (20 x 0.0125)? = 2.0625, For this case, there are four equations in four unknowns: 2.06256 6 = 300 @, +2.0625¢, % = 0 | +2.06250 ee) (hae fy +1.08750, 0 ‘The above tridiagonal system of algebra the Gauss elimination method (this is where we need a computer!). The solution is given by equations can be solved usin, {8} = (245.81, 206.98, 181.10, 166.52)" ‘The analytical solution at the same points is {0} = (248.75, 213.13, 190.90, 180.66)" The maximum error is about 7.8 percent. When the number of mesh points is doubled, the ‘maximum error is increased to 32. ,oes down to 4.2 percent, and it is 1 percent when the number of mesh points,Harner ts iwreopeeion 13; 1.4 THE FINITE ELEMENT METHOD 1.4.1 The Basic Idea ‘The finite element method is a numerical method like the finite difference method but is more general and powerful in its application to real-world problems that involve complicated physics, geometry, and/or boundary conditions. In the finite element method, a given domai is viewed as a collection of subdomains, and over each subdomain the governing equation is approximated by any of the traditional variational methods. The main reason behind seeking approximate solution on a collection of subdomains is the fact that it is easier to represent a complicated function as a collection of simple polynomials, as can be seen from Fig. 1.4.1. Of course, each individual segment of the solution should fit with its neighbors in the sense that the function and possibly derivatives up to a chosen order are continuous at the connecting points, These ideas will be more clear in the sequel. 1.4.2 The Basic Features The method is endowed with three distinct features that account for its superiority over other competing methods. First, a geometrically complex domain @ of the problem, such as the one in Fig. 1.4.2(a), is represented as a collection of geometrically simple subdomains, called finite elements [see Fig. 1.4.2(b)]. Each finite element Q, [Fig. 1.4.2(c)] is viewed s an independent domain by itself. Here the word “domain refers to the geometric region over which the equations are solved. Second, over each finite element, algebraic equations among the quantities of interest are developed using the governing equations of the problem. Third, the relationships from all elements are assembled (i.e., elements are put back into their original positions of the total domain) as indicated in Fig. 1.4.2(@), using certain interelement relationships. These statements may not make complete sense to the reader until some specific examples to illustrate the features are discussed in the coming pa Approximations enter engineering analyses at several stages. The division of the whole domain into finite elements may not be exact (i.e., the assemblage of elements, %, does not match the original domain ©), introducing error in the domain being modeled. The second stage is when element equations are derived. Typically, the dependent unknowns (w) of the problem are approximated using the basic idea that any continuous function can be ua) Laie Figure 1.4.1 Piecewise approximation of a function.14. ax ivrRoDUCTION 10 THE FINITE ELEMENT METHOD — Boundary, I Inter-element flux € 2 Representation of a two-dimensional domain by a collection of triangle laterals. Figure 1. represented by a linear combination of known functions , and undetermined coefficients ¢; (uw uy = Dcidy). Algebraic relations among the undetermined coefficients c; are obtained by satisfying the governing equations, in a weighted-integral sense, over each element. The proximation functions ¢, are often taken to be polynomials, and they are derived using concepts from interpolation theory. Therefore, they are termed interpolation functions. Thus, approximation errors in the second stage are introduced both in representing the solution « as well as in evaluating the integrals. Finally, errors are introduced in solving the assembled system of equations. Obviously, some of the errors discussed above can be zero. When all the errors are zero, we obtain the exact solution of the problet most problems). Next, the basic ideas and some terminol n (which is not the case in of the finite element method are introduced and the basic features rather via several simple examples. The reader should try to underst than to question the use of an approximation method to solve such a simple problem (which can possibly be solved exactly without using the finite element method) Example 1.4.1 (Approximation of the Perimeter of a Circle) This discussion is an expansion of an article written by the author in 1978 for a student magazine at the University of Oklahoma (Reddy, 1978). Consider the problem of determining the perimeter (quantity of interest) of a circle of radius R [see Fig. 1.4.3(a)]. We prerendeaapreR ts ireopuction 1516 2, in the mesh. For a typical element 2, he is given by [see Fig. 1.4.3(d)] 1 Resin 50, : sin 50. (4) where R is the radius of the circle and @, < 7 is the angle subtended by the line segment. The above equations are called element equations. Ancient mathematicians most likely made ‘measurements, rather than using Eq. (1.4.1), to find h, (they did not know 77). . Assembly of element equations and solution. The approximate value of the perimeter of the circle is obtained by putting together the clement properties in a meaningful way; this process is called the assembly of the element equations. It is based, in the present case, on the simple idea that the total perimeter of the polygon &, (assembly of elements) is equal to the sum of the lengths of individual elements: R= (4.2) Then P, represents an approximation to the actual perimeter, P. If the mesh is uniform, or ‘A, is the same for each of the elements in the mesh, then é, = 271/n, and we have (2K sin) (43) |. Convergence and error estimate. For this simple problem, we know the exact solution: P=2R. We can estimate the error in the approximation und show that the approximate solution P, converges to the exact value P in the limitasn + cc. Consider atypical element Q.. The error in the approximation is equal io the difference between the length of the sector and that of the line segment [see Fig. 1.4.3(d)] = 15, hel 4a) where 5, = Ré, is the length of the sector. Thus, the error estimate for an element in the mesh is given by (1.45) ‘The total error (called global error) is given by multiplying £, by n E=nk, -P, (1.4.6) 2R (x —nsin=) =27R- P= ‘We now show that E goes to zero as n > 00, Letting x = I/n, we have P,=2Rnsin = and (1.4.7)cuarrer 1: iwrropucron 1718 ax itRODUCTION To THE FINITE ELEMENT METHOD ‘The two examples considered above illustrates how the idea of piecewise approxima- tion is used to approximate irregular geometry and calculate required quantities. Thus, subdividing a geometrically complex domain into parts that allow the evaluation of desired quantities is a very natural and practical approach. The idea can be extended to approxi- mate functions representing physical quantities. For example, the temperature variation in a two-dimensional domain can be viewed as a curved surface, and it can be approximated over any part of the domain, ie., over a subdomain, by a polynomial of desired degree. A curved surface over a triangular subregion may be approximated by a planar surface, Approximate solution ‘True solution surface ‘on the subdomain ‘True solution surface ‘over the element domain “Approximate sofution —a plane (linear)— wy cover the clement domain Subdomain Original domain Figure 1.4.5 Approximation of a curved surface by a planar surface.CHAPTER 1: INrRODUCTION 19° as shown in Fig. 1.4.5, Such ideas form the basis of finite element approximations. The next example illustrates this idea for a one-dimensional continuous system described by a differential equation.20. aNivrRoDUCTION TO THE FINITE ELEMENT METHOD Equation (1.4.11) is subject to appropriate boundary conditions, for example, insulated atr = Ry and subjected to a temperature Ty at r= Ry’ Ri Ti=% a r=R (1.4.12) ity, which varies from layer to layer, R; and Ro are the ner and outer radii of the cylinder (Fig. 1.4,6(b)], and g is the rate of energy generation in the medium. When it is difficult to obtain an exact solution of Bs. (1.4.11) and (1.4.12), her because of complex geometry and material properties or because g(r) is a complicated function that does not allow exact evaluation of its integral, we seek an approximate one. In the finite element method, the domain (Rj, Ro) is divided into W subintervals (Fig. 1.4.6(0)], ‘and the approximate solution is sought in the form LT WPO, Ri srs Ri +h (first interval) Tey BO = OTP WP CO, Ri thi sr< Ri + hi + hy (second interval) wr) = OTM YW), Rit hy to + haya SoS Ro (Nth interval) (1.4.13) where A, denotes the length of the eth interval, 7; (the element label “e” is omitted in the interest of brevity) is the value of the temperature 7.(r) at the jth geometric point (called the node) of the eth interval, y; are the interpolation polynomials on the eth interval (finite element), and 7; are termed the nodes values (which are unknown). Thus, the continuous function T(r) is approximated in each finite element by a desired degree of interpolation polynomials, and the polynomial is expressed in terms of the values of the function ata selected number of points. in the interval. The number of points is equal to the number of parameters in the polynomial. For example, a linear polynomial approximation of the temperature over the interval requires two values, and hence two points, are identified in the interval. The endpoints of the interval are selected for this purpose because the two points, called nodes, also define the length of the interval (see Fig. 1.4.7). For higher-order polynomial approximations, additional nodes are Tr) To @ Figure 1.4.7 (a) Linear approximation T,(r) of a function 7'(r) over an element. (b) Quadratic, approximation 7.(r) of 7).‘cuarrer 1: iNrropucTION 21 1.4.3 Some Remarks Insummary, inthe finite element method a given domain is divided into subdomains, called finite elements, and an approximate solution to the problem is developed over each element. ‘The subdivision of a whole domain into parts has two advantages: 1, Allows accurate representation of complex geometry and inclusion of dissimilar material properties. 2. Enables easy representation of the total solution by functions defined within each element that capture local effects (e.g., large gradients of the solution), ‘The three fundamental steps of the finite element method that are illustrated via the examples are: 1, Divide the whole domain into parts (both to represent the geometry and solution of the problem). 2. Over each part, seek an approximation to the solution as a linear combination of nodal values and approximation functions, and derive the algebraic relations among the nodal values of the solution over each part. 3. Assemble the parts and obtain the solution to the whole. Although the above examples illustrate the basic idea of the finite element method, there are several other features that are either not present or not apparent from the examples. Some remarks are in order.28 . In general, the assemblage of finite elements AN INTRODUCTION TO THE FINITE ELEMENT METHOD ‘We can discretize the geometry of the domain, depending on its shape, into a mesh of more than one type of element (by shape or order). For example, in the approximation of an irregular domain, we can use a combination of rectangles and triangles. However, the element interfaces must be compatible in the sense that the solution is continuous. . If more than one type of element is used in the representation of the domain, one of each kind should be isolated and its equations developed. The governing equations of problems considered in this book are differential equations, In most cases, the equations cannot be solved over an element for two reasons. First, they do not permit the exact solution. It is here that the variational methods come into play. Second, the discrete equations obtained in the variational methods cannot be solved independent of the remaining elements because the assemblage of the elements is subjected to certain continuity, boundary, and/or initial conditions. . There are two main differences in the form of the approximate solution used in the finite element method and in the classical variational methods (i... variational methods applied to the whole domain). First, instead of representing the solution w as a linear combination (uj, = D>; ¢)¢;) in terms of arbitrary parameters c; as in the variational methods, in the finite element method the solution is often represented as a linear com- bination (w, =) ¢@s) im terms of the values w; of uy (and possibly its derivatives as well) at the nodal points. Second, the approximate functions in the finite element method are often polynomials that are derived using interpolation theory. However, the finite element method is not restricted to the use of approximations that are linear ‘combinations of nodal values 1; and interpolation functions y, that are algebraic poly- nomials. One may use nodeless variables and nonpolynomial functions to approximate a function (like in meshless or element-free methods). . The number and the location of the nodes in an element depend on (a) the geometry of the element, (b) the degree of the polynomial approximation, and (c) the weighted- integral form of the equations. By representing the required solution in terms of values at the nodes, we directly obtain the approximate solution at the nodes. The assembly of elements, in a general case, is based on the idea that the solution (and possibly its derivatives for higher-order equations) is continuous at the interelement boundaries. is subjected to boundary and/or initial conditions. The discrete equations associated with the finite element mesh are solved only after the boundary and/or initial conditions have been imposed, . There are three sources of error in a finite element solution: (a) those due to the ap- proximation of the domain (this was the only error present in the first two examples); (b) those due to the approximation of the solution; and (c) those due to numerical computation (¢.g., numerical integration and round-off errors in a computer). The es- timation of these errors, in general, is not simple. However, under certain conditions, they can be estimated for an element and problem, ). The accuracy and convergence of the finite element solution depends on the differential equation, its integral form, and the element used. “Accuracy” refers to the difference between the exact solution and the finite element solution, while “convergence” refers to the accuracy as the number of elements in the mesh is increased,‘CHAPTER |: INTRODUCTION — 23 Figure 1.4.8 Finite element discretization of a connecting rod and piston. 10, For time-dependent problems, a two-stage formulation is usually followed. In the first stage, the differential equations are approximated by the finite element method to obtain a set of ordinary differential equations in time. In the second, the differential equations time are solved exactly or further approximated by either variational methods or finite difference methods to obtain algebraic equations, which are then solved for the nodal values (see Chapter 6). 11. The desktop computers of today are more powerful than the supercomputers that existed when the finite element method was first implemented. Consequently, the analysis time is considerably reduced, provided the mesh used to model the problem is adequate (a finite clement mesh of a connecting rod and piston is shown in Fig. 1.4.8). Even the automatic mesh generation programs cannot guarantee meshes that are free of irregularly shaped elements or have sufficient number of elements in a region containing high gradients in the solution, both of which result in loss of accuracy or, in the case of nonlinear problems, nonconvergence of solutions. 12. Element-free methods that require no assembly (because there are no elements) are being developed. Such methods have applications in fracture mechanics and wave propagation problems where re-meshing is necessary. 1.4.4 A Brief Review of History and Recent Developments The idea of representing a given domain as a collection of discrete parts is not unique to the finite element method. It was recorded that ancient mathematicians estimated the value of x by noting that the perimeter of a polygon inscribed in a circle approximates the circumference of the latter. They predicted the value of zr to accuracy of almost 40 significant digits by representing the circle as a polygon of a finitely large number of sides (Reddy, 1978, 1993). In modern times, the idea found a home in aircraft structural analysis, where, for example, wings and fuselages are treated as assemblages of stringers, skins, and shear panels. In 1941, Hrenikoff [1941] introduced the so-called framework method in which a plane elastic medium was represented as a collection of bars and beams. The use of piecewise-continuous functions defined over a subdomain to approximate an unknown fune- tion can be found in the work of Courant (1943), who used an assemblage of triangular elements and the principle of minimum total potential energy to study the St. Venant tor- sion problem. Although certain key features of the finite element method can be found in eis