AN INTRODUCTION TO THE FINITE ELEMENT METHOD
McGrawHill Series in Mechanical Engineering Consulting Editors Jack P. H...
1152 downloads
3362 Views
22MB Size
Report
This content was uploaded by our users and we assume good faith they have the permission to share this book. If you own the copyright to this book and it is wrongfully on our website, we offer a simple DMCA procedure to remove your content from our site. Start by pressing the button below!
Report copyright / DMCA form
AN INTRODUCTION TO THE FINITE ELEMENT METHOD
McGrawHill Series in Mechanical Engineering Consulting Editors Jack P. Holman, Southern Methodist University John R. Lloyd, Michigan State University Anderson: Modern Compressible Flow: With Historical Perspective Arora: Introduction to Optimum Design Bray and Stanley: Nondestructive Evaluation: A Tool for Design, Manufacturing, and Service Culp: Principles of Energy Conversion Dally: Packaging of Electronic Systems: A Mechanical Engineering Approach Dieter: Engineering Design: A Materials and Processing Approach Eckert and Drake: Analysis of Heat and Mass Transfer Edwards and McKee: Fundamentals of Mechanical Component Design Gebhart: Heat Conduction and Mass Transfer Heywood: Internal Combustion Engine Fundamentals Hinze: Turbulence Howeil and Huckius: Fundamentals of Engineering Thermodynamics Hutton: Applied Mechanical Vibrations Juvinall: Engineering Considerations of Stress, Strain, and Strength Kane and Levinson: Dynamics: Theory and Applications Kays and Crawford: Convective Heat and Mass Transfer Kimbrell: Kinematics Analysis and Synthesis Marlin: Kinematics and Dynamics of Machines Modest: Radiative Heat Transfer Norton: Design of Machinery Phelan: Fundamentals of Mechanical Design Raven: Automatic Control Engineering Reddy: An Introduction to the Finite Element Method Rosenberg and Karnopp: Introduction to Physics Schlichting: BoundaryLayer Theory Shames: Mechanics of Fluids Shennan: Viscous Flow Shigley: Kinematic Analysis of Mechanisms Shigley and Mischke:. Mechanical Engineering Design Shigley and Uicker: Theory of Machines and Mechanisms Stifiler: Design with Microprocessors for Mechanical Engineers Stoecker and Jones: Refrigeration and Air Conditioning Ullman: The Mechanical Design Process Vanderplaats: Numerical Optimization: Techniques for Engineering Design, with Applications White: Viscous Fluid Flow Zeid: CAD/CAM Theory and Practice
Also Available from McGrawHill Schaum's Outline Series in Mechanical Engineering Most outlines include basic theory, definitions, and hundreds of solved problems and supplementary problems with answers. Titles on the Current List Include: ACOlLStics Basic Equations of Engineering Continuum Mechanics Engineering Economics Engineering Mechanics, 4th edition Fluid Dynamics, 2d edition Fluid Mechanics & Hydraulics, 2d edition Heat Transfer Introduction to Engineering Calculations Lagrangian Dynamics Machine Design Mathematical Handbfok of Formulas & Tables Mechanical Vibrations Operations Research Statistics & Mechanics of Materials Strength of Materials, 2d edition Theoretical Mechanics Thermodynamics, 2d edition
Schaum's Solved Problems Books Each title in this series is a complete and expert source of solved problems containing thousands of problems with worked out solutions. Related Titles on the Current List Include:
•
3000 Solved Problems in Calculus 2500 Solved Problems in Differential Equations 2500 Solved Problems in Fluid Mechanics and Hydraulics 1000 Solved Problems in Heat Transfer 3000 Solved Problems in Linear Algebra 2000 Solved Problems in Mechanical Engineering Thermodynamics 2000 Solved Problems in Numerical Analysis 700 Solved Problems in Vector Mechanics for Engineers: Dynamics 800 Solved Problems in Vector Mechanics for Engineers: Statistics
Available at your College Bookstore. A complete list of Schaum titles may be obtained by writing to: Schaum Division McGrawHill, Inc. Princeton Road, S1 Hightstown, NJ 08520
AN INTRODUCTION TO THE FINITE ELEMENT METHOD Second Edition
J. N. Reddy Oscar S. Wyatt Chair in Mechanical Engineering Texas A &M University College Station, Texas 77843
McGrawHill, Inc. New York St. Louis San Francisco Auckland Bogota Caracas Lisbon London Madrid Mexico Milan Montreal New Delhi Paris San Juan Singapore Sydney Tokyo Toronto
yr:; This book was set in Times Roman. The editors were John J. Corrigan and John M. Morriss; the production supervisor was Louise Karam. The cover was designed by Joseph Gillians. Project supervision was done by The Universities Press. R. R. Donnelley & Sons Company was printer and binder.
K
I
.
.
·'1"' .
LJ
TA347.F5 R4 1993 An lDtrodllcUoD to !he tIDItt ele.mtlltlDttho
1llllJlllfJ/lIlfmlflllflfIm/IIIII
'fl l ~mtlfI I IJl f
0020046676
269154
AN INIRODUCTION TO THE FINITE ELEMENT METHOD
Copyright © 1993, 1984by McGrawHill, Inc. All rights reserved. Printed in the United States of America. Except as permitted under the United States Copyright Act of 1976, no part of this publication may be reproduced or distributed in any form or by any means, or stored in a data base or retrieval system, without the prior written permission of the publisher.
234567890 DOCmOC 99876543 ISBN 0070513554 Library of Congress CataloginginPubllcatlen Data Reddy, J. N. (John N.) An introduction to the finite element method/J. N. Reddy.2nd ed. p. cm.{McGraw.Hill series in mechanical engineering) Includes bibliographical references and index. ISBN 007.QS13554 1. Finite element method. I. Series. TA347.F5R4 1993 9229532 620'.OOl'51535dc20
ABOUT THE AUTHOR
J. N. Reddy is the inaugural holder of the Oscar S. Wyatt, Jr, Endowed Chair
in Mechanical Engineering at Texas A&M: University. Prior to the current appointment, he was the Clifton C. Garvin Chaired Professor of Engineering Science and Mechanics at Virginia Polytechnic Institute and State University. After receiving his Ph.D. in Applied Mechanics in 1974, he joined Lockheed Missiles and Space Company in Huntsville, Alabama, as a research scientist. In 1975, he joined the School of Aerospace, Mechanical and Nuclear Engineering at the University of Oklahoma. In 1980, he was appointed as a professor of Engineering Science and Mechanics at Virginia Polytechnic Institute and State University. Dr. Reddy has authored and coauthored over 150 journal papers and seven books on the theory and applications of the finite element method in solid and structural mechanics, fluid mechanics and heat transfer, and on the development of refined theories and finite element models of laminated composite plates and shells. Dr. Reddy has taught many short courses on the finite element method to industry and government, and has received Certificates of Teaching Excellence from Virginia Polytechnic Institute and State University. Dr. Reddy serves on the editorial boards of a dozen journals, including Journal of Applied Mechanics, Meccanica, International Journal for Numerical Methods in Engineering, International Journal of Numerical Methods in Fluids, Journal of NonLinear Mechanics, and Journal of Composites Technology and Research. He has received the College of Engineering Research Award from the University of Oklahoma and the Alumni Research Award from the Virginia Polytechnic Institute and State University. Dr. Reddy was awarded the 1983 Walter L. Huber Research Prize of the American Society of Civil Engineers and 1992 Worcester Reed Warner Medal of the American Society of Mechanical Engineers. He is a Fellow of the American Academy of Meehanics, the American Society of Mechanical Engineers, the American Society of Civil Engineers, and the Aeronautical Society of India. vii
To MY FAMILY Aruna, Anita, and Anil
CONTENTS
Preface to the Second Edition Preface to the First Edition
xvii xix
Part 1 Preliminaries 1 Introduction 1.1 General Comments 1.2 Historical Background 1.3 The Basic Concept of the Finite Element Method 1.3.1 General Comments 1.3.2 Approximation of the Circumference of a Circle 1.3.3 Approximate Determination of the Center of Mass 1.3.4 Solution of Differential Equation 1.3.5 Some Remarks 1.4 The Present Study 1.5 Summary References for Additional Reading
2 Integral Formulations and Variational Methods 2.1 Need for WeightedIntegral Forms 2.2 Some Mathematical Concepts and Formulae 2.2.1 Boundary, Initial, and Eigenvalue Problems 2.2.2 Integral Relations 2.2.3 Functionals 2.2.4 The Variational Symbol 2.3 Weak Formulation of Boundary Value Problems 2.3.1 Introduction 2.3.2 WeightedIntegral and Weak Formulations 2.3.3 Linear and Bilinear Forms and Quadratic Functionals 2.3.4 Examples
3 3
5
5 5
6 8 10 13
15 15
16 18 18 20
20 22
26 27 28
28 28
33
35 ix
X
COlITENTS
2.4 Variational Methods of Approximation 2.4.1 Introduction 2.4.2 The RayleighRitz Method 2.4.3 The Method of Weighted Residuals 2.5 Summary Problems References for Additional Reading
40 40 40 51 57 59 63
Part 2 Finite Element Analysis of OneDimensional Problems 3 SecondOrder Boundary Value Problems 3.1 Introduction . 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 Postprocessing of the Solution 3.2.8 Radially Symmetric Problems 3.3 Applications 3.3.1 Heat Transfer 3.3.2 Fluid Mechanics 3.3.3 Solid Mechanics 3.4 Summary Problems References for Additional Reading
4 Bending of Beams 4.1 Introduction 4.2 The EulerBernoulli Beam Element 4.2.1 Governing Equation 4.2.2 Discretization of the Domain 4.2.3 Derivation of Element Equations 4.2.4 Assembly of Element Equations 4.2.5 Imposition of Boundary Conditions 4.2.6 Solution 4.2.7 Postprocessing of the Solution 4.2.8 Examples 4.3 Plane Truss and EulerBernoulli Frame Elements 4.4 The Timoshenko Beam and Frame Elements 4.4.1 Governing Equations 4.4.2 Weak Form 4.4.3 Finite Element Model 4.5 Inclusion of Constraint Equations 4.6 S u m m a r y 
67 67 70 70 72 72 89 95 95
96 103 105 105 117 123 127 128 141 143 143 143 143
144 144 151 154 156 158 160
167 177 177 177 178 187 191
CONTENTS
xi
Problems References for Additional Reading
192 198
5 Finite Element Error Analysis
199 199 200 201 202 207 207 208
5.1 Approximation Errors 5.2 Various Measures of Errors 5.3 Convergence of Solution 5.4 Accuracy of the Solution 5.5 Summary Problems References for Additional Reading
6 Eigenvalue and TimeDependent Problems 6.1 Eigenvalue Problems 6.1.1 Introduction 6.1.2 Formulation of Eigenvalue Problems 6.1.3 Finite Element Models 6.1.4 Applications 6.2 TimeDependent Problems 6.2.1 Introduction 6.2.2 Semidiscrete Finite Element Models 6.2.3 Time Approximations 6.2.4 Mass Lumping 6.2.5 Applications 6.3 Summary Problems References for Additional Reading
7 Numerical Integration and Computer Implementation 7.1 Isoparametric Formulations and 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 7.2.3 Preprocessor 7.2.4 Calculation of Element Matrices (Processor) 7.2.5 Assembly of Element Equations (Processor) 7.2.6 Imposition of Boundary Conditions (Processor) 7.2.7 Solution of Equations and Postprocessing 7.3 Applications of the Computer Program FEMIDV2 7.3.1 General Comments 7.3.2 Illustrative Examples 7.4 Summary Problems References for Additional Reading
209 209 209 210 213 216 224 224 225 227 232 233 241 241 245
246 246 246 248 249 251 251 258 258
259
260 262
265 267 269
270 270 271
286 287 291
xii
CONTENTS
Part 3 Finite Element Analysis of TwoDimensional Problems 8
SingleVariable Problems 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 Interpolation Functions 8.2.6 Evaluation of Element Matrices and Vectors 8.2.7 Assembly of Element Equations 8.2.8 Postprocessing 8.2.9 Axisymmetric Problems 8.2.10 An Example 8.3 Some Comments on Mesh Generation and Imposition of Boundary Conditions
8.3.1 Discretization of a Domain 8.3.2 Generation of Finite Element Data 8.3.3 Imposition of Boundary Conditions 8.4 Applications 8.4.1 Heat Transfer 8.4.2 Fluid Mechanics 8.4.3 Solid Mechanics 8.5 Eigenvalue and TimeDependent Problems 8.5.1 Introduction 8.5.2 Parabolic Equations 8.5.3 Hyperbolic Equations 8.6 Summary Problems References for Additional Reading
9 Interpolation Functions, Numerical Integration, and Modeling Considerations 9.1 Library of Elements and Interpolation Functions 9.1.1 Introduction 9.1.2 Triangular Elements 9.1.3 Rectangular Elements 9.1.4 The Serendipity Elements 9.2 Numerical Integration 9.2.1 Preliminary Comments 9.2.2 Coordinate Transformations 9.2.3 Integration over a Master Rectangular Element 9.2.4 Integration over a Master Triangular Element 9.3 Modeling Considerations 9.3.1 Preliminary Comments 9.3.2 Element Geometries 9.3.3 Mesh Generation 9.3.4 Load Representation9.4 Summary
295 295 2Cfl 297 298 299 301 303 311 318 322 323 324 334 334 336 339 340 340 353 365 370 370 372 379 384 386 402 404 404 404
405 411 417 421 421 423 429 433 439 439 439 441 446 448
COl'ITENTS
Problems References for Additional Reading
10 Plane Elasticity 10.1 Introduction 10.2 Governing Equations 10.2.1 Assumptions of Plane Elasticity 10.2.2 Basic Equations 10.3 Weak Formulations 10.3.1 Preliminary Comments 10.3.2 Principle of Virtual Displacements in Matrix Form 10.3.3 Weak Form of the Governing Differential Equations lOA Finite Element Model 1004.1 Matrix Form of the Model 1004.2 Weak Form Model 1004.3 Eigenvalue and Transient Problems 10.5 Evaluation of Integrals 10.6 Assembly and Boundary and Initial Conditions 10.7 Examples 10.8 Summary Problems References for Additional Reading
11 Flows of Viscous Incompressible Fluids 11.1 11.2 11.3 11.4
Preliminary Comments Governing Equations VelocityPressure Finite Element Model PenaltyFinite Element Model 11.4.1 Penalty Function Method 11.4.2 Formulation of the Flow Problem as a Constrained Problem 11.4.3 Lagrange Multiplier Formulation 11.4.4 Penalty Function Formulation 11.4.5 Computational Aspects 11.5 Examples 11.6 Summary Problems References for Additional Reading
12 Bending ofElastic Plates 12.1 Introduction 12.2 Classical Plate Model 12.2.1 Displacement Field 12.2.2 Virtual Work Statement 12.2.3 Finite Element Model 12.3 Shear Deformable Plate Model 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 12.4 Eigenvalue and TimeDependent Problems
xiii 448 453
455 455 456 456 457 459 459 459 461 461 461 464 465 465 468 469 476 476 480 482 482 483 484 488 488 489 491 491 492 494 502 503· 506 508 508 510 510 510 514 516 516 517 518 519 520
xiv
CONlENTS
12.5 Examples 12.6 Summary Problems References for Additional Reading
522 529 529 531
13 Computer Implementation
533 533 534 535
13.1 13.2 13.3 13.4
Introduction Preprocessor Element Computations: Processor Applications of the Computer Program FEM2DV2 13.4.1 Introduction 13.4.2 Description of Mesh Genrators 13.4.3 Applications (Illustrative Examples) 13.5 Summary Problems References for Additional Reading
540 540 548 551 563
570 575
Part 4 Advanced Topics 14 WeightedResidual Finite Element Models, and Finite Element Models of Nonlinear and ThreeDimensional Problems 14.1 Introduction 14.2 Alternative Formulations 14.2.1 Introductory Comments 14.2.2 WeightedResidual Finite Element Models 14.2.3 Mixed Formulations 14.3 Nonlinear Problems 14.3.1 General Comments 14.3.2 LargeDeflection Bending of (EulerBernoulli) Beams 14.3.3 Solution Methods for Nonlinear Algebraic Equations 14.3.4 The 20 NavierStokes Equations 14.4 ThreeDimensional Problems 14.5 Summary Problems References for Additional Reading
579
579 580 580 580 590 594
594 595 597
598 599
601 602
606
Appendixes 1 Fortran Listing of FEMIDV2
609
2 Fortran Listing of FEM2DV2
640
Index
679
PREFACE TO THE SECOND EDITION
This second edition has the same objectives as the first, namely, an introduction to the finite element method as applied to linear, one and twodimensional problems of engineering and applied sciences. The revisions are mainly in the form of additional details, expansion of the topics discussed, and the addition of a few topics to make the coverage more complete. The major organizational change from the first edition is the division of its five chapters into fourteen chapters here. These chapters are grouped into four parts. This reorganization should aid instructors in selecting suitable material for courses. Other organizational changes include putting problem sets at the ends of the chapters, providing a chapter summary for each, and reviewing pertinent equations and text in each chapter instead of referring to several chapters back. In addition, example problems in Chapters 3 and 8 are presented in separate sections on heat transfer, fluid flow, and solid mechanics. Additional details are provided on the construction of the weak forms, time approximations (e.g., accuracy and stability of schemes, and mass lumping), alternative finite element formulations, and nonlinear finite element models. The new topics include sections on trusses and frames, the Timoshenko beam element, eigenvalue problems, and classical plate bending elements. All these changes are also reflected in the revised computer programs FEMIDV2 and FEM2DV2 (revised versions of the FEMID, FEM2D and PLATE programs in the first edition). Therefore the sections on computer implementation and applications of FEMIDV2 and FEM2DV2 have also been modified extensively. These changes are accompanied by the addition of several figures, tables, and examples. These extensive changes have resulted in a second edition that is 60% larger. In the interest of keeping the cost of the book within reasonable limits xv
xvi
PREFACE TO THE SECOND EDITION
while retaining the basic approach and technical details, certain portions of the original manuscript have been omitted. More specifically, answers to selective problems have been included at the end of the problem statements themselves, rather than in a separate section. Interested readers and instructors can obtain a copy of the excutable programs on a diskette from the author. Fortran source programs can also be purchased from the author. ' There is no doubt that this edition is more complete and thorough than the first. It can be used as a textbook for an introductory and/or intermediate level course on the finite element method at senior undergraduate as well as graduate levels. Students of engineering and applied sciences should feel comfortable with the coverage in the book. The author gratefully acknowledges help in reading the manuscript and suggestions for constructive changes from several colleagues. These include: Hasan Akay, Purdue University at Indianapolis, Norman Knight, Jr, Clemson University; J. K. Lee, Ohio State University; William Rule, University of Alabama; Martin Sadd, University of Rhode Island; John Whitcomb, Texas A&M University, and the author's research students: Ronald Averill, Filis Kokkinos, Y. S. N. Reddy, and Donald Robbins. It is a great pleasure to acknowledge typing of the manuscript by Mrs Vanessa McCoy, without whose patience and cooperation this work would not have been completed. J. N. Reddy
PREFACE TO THE FIRST EDITION
The motivation which led to the writing of the present book has come from my many years of teaching finiteelement courses to students from various fields of engineering, meteorology, geology and geophysics, physics, and mathematics. The experience gained as a supervisor and consultant to students and colleagues in universities and industry, who have asked for explanations of the various mathematical concepts related to the finiteelement method, helped me introduce the method as a variationally based technique of solving differential equations that arise in various fields of science and engineering. The many discussions I have had with students who had no background in solids and structural mechanics gave rise to my writing a book that should fill the rather unfortunate gap in the literature. The book is designed for senior undergraduate and firstyear graduate students who have had a course in linear algebra as well as in differential equations. However, additional courses (or exposure to the topics covered) in mechanics of materials, fluid flow, and heat transfer should make the student feel more comfortable with the physical examples discussed in the book. In the present book, the finiteelement method is introduced as a variationally based technique of solving differential equations. A continuous problem described by a differential equation is put into an equivalent variational form, and the approximate solution is assumed to be a linear combination, ~ Cj4Jj, of approximation functions 4Jj. The parameters Cj are determined using the associated variational form. The finiteelement method provides a systematic technique for deriving the approximation functions for simple subregions by which a geometrically complex region can be represented. In the finiteelement method, the approximation functions are piecewise polynomials (i.e., polynomials that are defined only on a subregion, called an element). I
xvii
xviii
PREFACE TO DIE FIRST EDITION
The approach taken in the present book falls somewhere in the middle of the approaches taken in books that are completely mathematical and those approaches that are more structuralmechanicsoriented. From my own experience as an engineer and selftaught applied mathematician, I know how unfortunate outcomes may be arrived at if one follows a "formula" without deeper insight into the problem and its approximation. Even the best theories lead ultimately to some sort of guidelines (e.g., which variational formulation is suitable, what kind of element is desirable, what is the quality of the approximation, etc.). However, without a certain theoretical knowledge of variational methods one cannot fully understand various formulations, finiteelement models, and their limitations. In the present study of variational and finiteelement methods, advanced mathematics are intentionally avoided in the interest of simplicity. However, a minimum of mathematical machinery that seemed necessary is included in Chapters 1 and 2. In Chapter 2, considerable attention is devoted to the construction of variational forms since this exercise is repeatedly encountered in the finiteelement formulation of differential equations. The chapter is concerned with two aspects: first, the selection of the approximation functions that meet the specified boundary condtions; second, the technique of obtaining algebraic equations in terms of the undetermined parameters. Thus, Chapter 2 not only equips readers with certain concepts and tools that are needed in Chapters 3 and 4, but it also motivates them to consider systematic methods of constructing the approximation functions, which is the main feature of the finiteelement method. In introducing the finite element method in Chapters 3 and 4, the traditional solid mechanics approach is avoided in favor of the "differential equation" approach, which has broader interpretations than a single special case. However, when specific examples are considered, the physical background of the problem is stated. Since a large number of physical problems are described by second and fourthorder ordinary differential equations (Chapter 3), and by the Laplace operator in two dimensions (Chapter 4), considerable attention is devoted to the finiteelement formulation, the derivation of the interpolation functions, and the solution of problems described by these equations. Representative examples are drawn from various fields of engineering, especially from heat transfer, fluid mechanics, and solid mechanics. Since this book is intended to serve as a textbook for a first course on the finiteelement method, advanced topics such as nonlinear problems, shells, and threedimensional analyses are omitted. Since the practice of the finiteelement method ultimately depends on one's ability to implement the technique on a digital computer, examples and exercises are designed to let the reader actually compute the solutions of various problems using computers. Ample discussion of the computer implementation of the finiteelement method is given in Chapters 3 and 4. Three model programs (FEM1D, FEM2D, and PLATE) are described, and their application is illustrated via several examples. The computer programs are very easy to understand because they are designed along the same lines as the
PREFACE TO THE fIRST EDITION
xix
theory presentedin the book. The programs are available for mainframe and IBM PC compatibles from the author for a small charge. Numerous examples, most of which are applications of the concepts to specific problems in various fields of engineering and applied science, are provided throughout the book. The conclusion of the examples are indicated by the symbol II. At approprate intervals in the book an extensive number of exercise problems is included to test and extend the understanding of the concepts discussed. For those who wish to gain additional knowledge of the topics covered in the book, many reference books and research papers are listed at the end of each chapter. There are several sections that can be skipped in a first reading of the book (such sections are marked with an asterisk); these can be filled in wherever needed later. The material is intended for a quarter or a semester course, although it is better suited for a semester course. The following schedule of topics is suggested for a first course using the present textbook: Undergraduate
Graduate
Chapter 1
Selfstudy
Chapter 1
Selfstudy
Chapter 2
Section 2.1 (self) Section 2.2 Sections 2.3.12.3.3
Chapter 2
Section 2.1 (self) Section 2.2 Section 2.3
Chapter 3
Sections 3.13.4 Sections 3.63.7
Chapter 3
Sections 3.13.7
Chapter 4
Sections 4.14.4 Section 4.7 Sections 4.8.14.8.4
Chapter 4
Sections 4.14.8
Chapter 5 Term Paper
Due to the intimate relationship between Sections 3.5 and 4.6, 3.6 and 4.7, and 3.7 and 4.8, they can be covered simultaneously. Also, it is suggested that Sections 3.6 and 3.7 (hence, 4.7 and 4.8) be covered after Section 3.2. The author wishes to thank all those students and colleagues who have contributed by their advice and criticism to the improvement of this work. The author is also thankful to Vanessa McCoy for skillful typing of the manuscript, to Mr. N. S. Putcha and Mr. K. Chandrashekhara for proofreading the pages, and to the editors Michael Slaughter and Susan Hazlett for their help and cooperation in publishing the mauscript. J. N. Reddy Tejashwina uadheetamasthu (May what we study be well studied)
PART
1 PRELIMINARIES
1
CHAPTER
1 INTRODUCTION
1.1 GENERAL COMMENTS Virtually every phenomenon in nature, whether biological, geological, or mechanical, can be described with the aid of the laws of physics, in terms of algebraic, differential, 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 seawater or in the atmosphere, and simulating weather in an attempt to understand and predict the mechanics of formation of tornadoes and thunderstorms are a few examples of many important practical problems. Most engineers and scientists studying physical phenomena are involved with two major tasks:
1. Mathematical formulation of the physical process 2. Numerical analysis of the mathematical model The mathematical formulation of a physical process requires background in related subjects (e.g., laws of physics) and, most often, certain mathematical tools. The formulation results in mathematical statements, often differential
3
4
PRELIMINARIES
equations, relating quantities of interest in the understanding and/or design of the physical process. Development of the mathematical model of a process is achieved through assumptions concerning how the process works. In a numerical simulation, we use a numerical method and a computer to evaluate the mathematical model and estimate the characteristics of the process. While the derivation of the governing equations for most problems is not unduly difficult, their solution by exact methods of analysis is a formidable task. In such cases, approximate methods of analysis provide alternative means of finding solutions. Among these, the finite difference method and the variational methods such as the RayleighRitz and Galerkin methods are most frequently used in the literature. In the finite difference approximation of a differential equation, the derivatives in the latter are replaced by difference quotients (or the function is expanded in a Taylor series) that involve the values of the solution at discrete mesh points of the domain. The resulting algebraic equations are solved, after imposing the boundary conditions, for the values of the solution at the mesh points. In the solution of a differential equation by a variational method, the equation is put into an equivalent weightedintegral form and then the approximate solution over the domain is assumed to be a linear combination (L.jc/Pj) of appropriately chosen approximation functions if'j and undetermined coefficients, c/. The coefficients Cj are determined such that the integral statement equivalent to the original differential equation is satisfied. Various variational methods, e.g., the RayleighRitz, Galerkin, and leastsquares methods, differ from each other in the choice of the integral form, weight functions, and/or approximation functions. A more complete discussion of variational methods will be given in Chapter 2. They suffer from the disadvantage that the approximation functions for problems with arbitrary domains are difficult to construct. The finite element method overcomes the disadvantage of the traditional variational methods by providing a systematic procedure for the derivation of the approximation functions over subregions of the domain. The method is endowed with three basic features that account for its superiority over other competing methods. First, a geometrically complex domain of the problem is represented as a collection of geometrically simple subdomains, called finite elements. Second, over each finite element, the approximation functions are derived using the basic idea that any continuous function can be represented by a linear combination of algebraic polynomials. Third, algebraic relations among the undetermined coefficients (Le., nodal values) are obtained by satisfying the governing equations, often in a weightedintegral sense, over each element. Thus, the finite element method can be viewed, in particular, as an elementwise application of the RayleighRitz or weightedresidual methods. In it, the approximation functions are often taken to be algebraic polynomials, and the undetermined parameters represent the values of the solution at a finite number of preselected points, called nodes, on the boundary and in the interior of the element. The approximation functions are
INTRODUCTION
5
derived using concepts from interpolation theory, and are therefore called interpolation functions. One finds that the degree of the interpolation functions depends on the number of nodes in the element and the order of the differential equation being solved.
1.2 HISTORICAL BACKGROUND 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 lr by noting that the perimeter of a polygon inscribed in a circle approximates the circumference of the latter. They predicted the value of lr to accuracies of almost 40 significant digits by representing the circle as a polygon of a finitely large number of sides. 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 introduced the socalled framework method, in which a plane elastic medium was represented as a collection of bars and beams. The use of piecewisecontinuous functions defined over a subdornain to approximate an unknown function 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 torsion problem. Although certain key features of the finite element method can be found in the works of Hrenikoff (1941) and Courant (1943), its formal presentation is attributed to Argyris and Kelsey (1960) and Turner, Clough, Martin, and Topp (1956). The term "finite element" was first used by Clough in 1960. Since its inception, the literature on finite element applications has grown exponentially, and today there are numerous journals that are primarily devoted to the theory and application of the method. A review of the historical developments and the basic theory of the finite element method can be found in more than three dozen textbooks that are exclusively devoted to its introduction and application. The selective finite element books listed in References for Additional Reading at the end of this chapter are only for additional information on certain topics (e.g., threedimensional problems, shells, structural dynamics, plasticity, and mathematics of finite elements). For the beginner, it is not necessary to consult these; the present book provides complete details of the method as applied to linear field problems, with examples from fluid mechanics, heat transfer, and solid mechanics.
1.3 THE BASIC CONCEPT OF THE FINITE ELEMENT METHOD 1.3.1 General Comments The most distinctive feature of the finite element method that separates it from others is the division of a given domain into a set of simple subdomains, called finite elements. Any geometric shape that allows computation of the solution or its approximation, or provides necessary relations among the values of the
6
PRELIMINARIES
solution at selected points, called nodes, of the subdomain, qualifies as a finite element. Other features of the method include seeking continuous, often polynomial, approximations of the solution over each element in terms of nodal values, and assembly of element equations by imposing the interelemerit continuity of the solution and balance of interelement forces. Here the basic ideas underlying the finite element method are introduced via two simple examples: 1. Determination of the circumference of a circle using a finite number of line
segments 2. Determination of the center of mass (or gravity) of an irregular body The first example is an expansion of an article written by the author in 1978 for a student magazine at the University of Oklahoma. Ideas expressed in the second can be found in books on statics of rigid bodies.
1.3.2 Approximation of the Circumference of a Circle Consider the problem of determining the perimeter of a circle of radius R (see Fig. 1.1a). Ancient mathematicians estimated the value of the circumference
(a)
(b)
,,/ I'
I
Se
f... Oe~
~
Element length h,
:~.,
I
L.
L
\ \
I
R
\
" ... 
/ .,."
(e)
" J
,/
FIGURE 1.1 Approximation of the circumference of a circle by line elements: (4) Circle of radius R; (b) uniform and nonuniform meshes used to represent the circumference of the circle; (c) a typical element.
INTRODUCTION
7
by approximating it by line segments, whose lengths they were able to measure. The approximate value of the circumference is obtained by summing the lengths of the line segments used to represent it. Although this is a trivial example, it illustrates several (but not all) ideas and steps involved in the finite element analysis of a problem. We outline the steps involved in computing an approximate value of the circumference of the circle. In doing so, we introduce certain terms that are used in the finite element analysis of any problem.
1. Finite element discretization. First, the domain (i.e., the circumference of the circle) is represented as a collection of a finite number n of subdomains, namely, line segments. This is called discretization of the domain. Each subdomain (i.e., line segment) is called an element. The collection of elements is called the finite element mesh. The elements are connected to each other at points called nodes. In the present case, we discretize the circumference into a mesh of five (n = 5) line segments. The line segments can be of different lengths. When all elements (i.e., line segments) are of the same length, the mesh is said to be uniform; otherwise, it is called a nonuniform mesh (see Fig. 1.Ib). 2. Element equations. A typical element (i.e., line segment, Qe) is isolated and its required properties, i.e., length, are computed by some appropriate means. Let he be the length of element Qe in the mesh. For a typical element Qe, he is given by (see Fig. 1.Ie) h e=2Rsin!Be
(1.1)
where R is the radius of the circle and Be < st is the angle subtended by the line segment. The above equations are called element equations. Ancient mathematicians most likely made measurements, rather than using (1.1), to find he. 3. Assembly of element equations and solution. The approximate value of the circumference (or perimeter) of the circle is obtained by putting together the element 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 (assembled elements) is equal to the sum of the lengths of individual elements: (1.2) Then P" represents an approximation to the actual perimeter, p. If the mesh is uniform, or he is the same for each of the elements in the mesh, then Be = 2n/n, and we have.
P"
=
n(2R sin~)
(1.3)
4. Convergence and error estimate. For this simple problem, we know the exact solution: p = 2nR. We can estimate the error in the approximation and show that the approximate solution P" converges to the exact p in the
8
PRELIMINARIES
limit as n + 00, Consider the typical element g". The error in the approximation is equal to the difference between the length of the sector and that of the line segment (see Fig. l.le): (1.4) where S" = RO" is the length of the sector. Thus, the error estimate for an element in the mesh is given by
n)
2n . E,,= R ( 2sm
n
(1.5)
n
The total error (called global error) is given by multiplying E; by n:
E = 2R ( n  n sin ~) = 2nR  P"
(1.6)
We now show that E goes to zero as n + 00, Letting x = lin, we have
nx
. n sin P" = 2Rn sin  = 2 R  
n
x
and
nx) . (
nx)
. P = lim . ( 2R sin  = lim 2nR cos  = 2n:R hm n XloO X .x~O 1
(1.7)
n""""":O
Hence, En goes to zero as n+ 00. This completes the proof of convergence. In summary, it is shown that the circumference of a circle can be approximated as closely as we wish by a finite number of piecewiselinear functions. As the number of elements is increased, the approximation improves, i.e., the error in the approximation decreases.
1.3.3 Approximate Determination of the Center of Mass Another elementary example to illustrate the finite element concept is provided by the calculation of the center of mass of a continuous body. It should be recalled, from a first course on statics of rigid bodies, that the calculation of the center of an irregular mass or the centroid of an irregular volume makes use of the socalled method of composite bodies, in which a body is conveniently divided (mesh discretization) into several parts (elements) of simple shape for which the mass and the center of mass (element properties) can be computed readily, The center of mass of the whole body is then obtained using the moment principle of Varignon (a basis for the assembly of el~ment properties):
(m i + m2 + ... + mn)X = mix t + m2x2 + , , . + mnxn
(1.8)
where X is the x coordinate of the center of mass of the whole body, me is the mass of the eth part, and x" is the x coordinate of the center of mass of the eth
INTRODUCTION
9
part. Similar expressions hold for the y and z coordinates of the center of mass of the whole body. Analogous relations hold for composite lines, areas, and volumes, wherein the masses are replaced by lengths, areas, and volumes, , respectively. When a given body is not expressible in terms of simple geometric shapes (elements) for which the mass and the center of mass can be represented mathematically, it is necessary to use a method of approximation to represent the properties of an element. As an example, consider the problem of finding the centroid (X, Y) of the irregular area (region) shown in Fig. 1.2. The region can be divided into a finite number of rectangular strips (elements), a typical element having width he and height be. The area of the eth strip is given by A e = hebe. The area A e is an approximation of the true area of the element because be is an estimated average height of the element. The coordinates of the centroid of the region are obtained by applying the moment principle:
EA..xe 
e
x=EA ' e
EAeYe 
e
y=
EA e
e
e
where xe and Ye are the coordinates of the centroid of the eth element with respect to the coordinate system used for the whole body. When the center of mass is required, A e in the above equations is replaced by the mass me = PeAe, p, being the mass density of the eth element; for a homogeneous body, Pe is the same for all elements. It should be noted that the accuracy of the approximation will be improved by increasing the number of strips (decreasing their width) used. Rectangular elements are used in the present discussion for the sake of simplicity only; one may choose to use elements of any size and shape that approximate the given area to a satisfactory accuracy. For example, a trapezoidal element will require two heights to compute the area:
A e = !he(b e + be+ 1) I f/'
y
he
~
Gill ::;T
=ri"
:I!i
be
:.".:
1:

:
i; j,
f.. . i , 
•
r y
~
X x FIGURE 1.2 Approximate determination of the mass or geometric centroid of an irregular region by dividing it into a set of rectangular or trapezoidal subregions.
10
PRELIMINARIES
where be and b e+ 1 are the left and right heights, respectively, of the eth element. The two examples considered above illustrate how the idea of piecewise approximation is used to approximate irregular geometries and calculate required quantities. In the first example, the circumference of a circle is approximated by a collection of line segments, whose measure is available. In the second, the geometric centroid or mass centroid of an irregular domain is located by approximating its geometry as a collection of strips that allow computation of their areas. Rectangles and trapezoids provide examples of the element geometries. 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 approximate functions representing physical quantities. For example, the temperature variation in a twodimensional domain can be viewed as a curved surface, and it can be approximated over any part of the domain, i.e., over a subdornain or element, by a function of desired degree. Figure 1.3 shows a curved surface over a triangular subregion approximated by a planar surface, i.e., a linear polynomial. Such ideas form the basis of finite element approximations. The next example illustrates this idea for a onedimensional continuous system.
1.3.4 Solution of Differential Equation Consider the temperature variation in a composite cylinder consisting of two coaxial layers in perfect thermal contact (see Fig. 1.4). Heat dissipation from a wire "(With two insulations) carrying an electric current and heat flow across a thickwalled composite circular cylindrical tube are typical examples. The temperature T is a function of the radial coordinate r. The variation of T with r is, in general, nonuniform. We wish to determine an approximation 1;,(r) to T(r) over the thicknesses of the cylinder. The exact solution is determined by solving the differential equation
1rdr dj dr
d ( rk =q(r) 
(1. 9a)
Curved surface
FIGURE 1.3
Approximation of curved surface (or a nonuniform function) over a triangular region by a planar surface.
INTRODUCTION
R0IIFdJ~~
.. r
Element
(b)
11
FIGURE 1.4 (a) Coaxial (composite) cylinder made of two different materials. (b) Finite element representation of a radial line of the cylinder.
subject to appropriate boundary conditions, for example, insulated at r = R, and subjected to a temperature To at r = R o: dT krdr=O
atr=R;;
T(r)=1'"
atr=R o
(1.9b)
where k is the thermal conductivity, which varies from layer to layer, R, and R: are the inner and outer radii of the cylinder, and q is the rate of energy generation in the medium. Note that the temperature is independent of the circumferential coordinate (because of the axisymmetric geometry, boundary conditions, and loading), and it has the same variation along any radial line. When it is difficult to obtain an exact solution of the problem (1.9), either because of complex geometry and material properties or because q(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 (R i J R o ) is divided into N subintervals, and the approximate solution is sought in the form n
11(r) = 2: Tft/lJ(r) (R, ~ r ~ R, + hI; first interval) j=I
n
T2(r) = 2: TJt/lJ(r) . (R, + hI ~ r ~Ri + hI + h 2 ; second interval) (1.10) j~l
n
TN(r) =
2: Tft/lf(r) j=I
(R, + h, + ... + hN  1 ~ r ~Ro; Nth interval)
12
PRELIMINARIES
where he denotes the length of the eth interval, TJ is the value of the temperature T,,(r) at the jth geometric point of the eth interval, and 1J1j are polynomials on the eth interval. The continuous function T(r) is approximated in each interval by a desired degree of polynomial, and the polynomial is expressed in terms of the values of the function at a 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 also define the length of the interval (see Fig. 1.5a). For higherorder polynomial approximation, additional points are identified interior to the interval (see Fig. 1.5b). The intervals are called finite elements, the points used to express the polynomial approximation of the function are called nodes, T'j are called nodal values, and 1J1'j are called finite element approximation functions. The nodal values TJ are determined such that 1;,(r) satisfies the differential equation (1.9a) and boundary conditions (1.9b) in some sense. Usually, the differential equation is satisfied in a weightedintegral sense, and boundary conditions on the function itself are satisfied exactly. The piecewise (i.e., elementwise) approximation of the solution allows us to include any discontinuous data, such as the material properties, and to use meshes of many lowerorder elements or a mesh of few higherorder elements to represent large gradients of the solution. Polynomial approximations of the form (1.10) can be derived systematically for any assumed degree of variation. The satisfaction of the differential equation in a weightedintegral sense leads, for steadystate problems, to algebraic relations among nodal temperatures Tj and heats Qj of the element. The algebraic equations of all elements are assembled (i.e., related to each other) such that the temperature is continuous and the heats are balanced at nodes common to elements. The Linear, T,
T(r)
I~
T 'I [
i 19
h '
.I T 3) will be zero. Figure 2.4 shows a comparison of the RayleighRitz solution with the exact solution.
lmEGRAL fORMULATIONS AND VARIATIONAL METIlODS
49
2.0.. Threeparameter solution  Twoparameter solution 1.0  Oneparameter solution
0.0 IV
!..='
x llY
_ dw X 103
dx
dw
 X
1.0
dx
103
2.0
3.0
+~r~rr''"'"r.l
0.2
0.0
0.6
0.4
0.8
1.0
1.2
xlL FIGURE 2.4 Comparison of the RayleighRitz solution with the exact solution of a cantilever beam under a uniform transverse load (EulerBernoulli beam theory).
The next example deals with twodimensional heat conduction in a square region. Note that the dependent variable, namely the temperature, is denoted by T, consistent with the standard notation used in heat transfer books. Example 2.6. Consider the Poisson equation in a unit square region:
kV2T=qo in Q={(x,y):O«x,y)i = Xi,
=
0
the Ritz coefficients are
C1 =
62
PRELIMINARIES
2.13. Set up the equations for the Nparameter RayleighRitz approximation of the following equations associated with a simply supported beam and subjected to a uniform transverse load I = 10: 2
2
W) Ja
d ( EI ddx = dx 2 2 W
2.14. 2.15. 2.16.
2.17.
d2w
= EI dx 2
=
for 0 < x < L
0 at x
= 0, L
(a) Use algebraic polynomials. (b) Use trigonometric functions. Compare the twoparameter RayleighRitz solutions with the exact solution. Answer: (a) Ct = 7«L, C2= 24' (
fif = (X;?  (Xk? Y~
= (xi  xD,
(3.21)
82
FINITE ELEMENT ANALYSIS OF ONE· DIMENSIONAL PROBLEMS
1
2
¢
~.f
3
6
0
h.
~
(a)
/
12
L.
.J ethelement
1he
rJ;o., i .. ~ +1 " he+l
~(e
11+2
+
.. l)thelement
(c)
FIGURE 3.5 Onedimensional Lagrange quadratic element and its interpolation functions: (a) geometry of the element; (b) interpolation functions; (c) global interpolation functions corresponding to the quadratic interpolation functions. Here I denotes the global node number, e the element number, and i the element node number.
and (3.19) takes the form 3
U' (r) = tJ.'i(x )Ul + tJ.'Hx )u~ + tJ.'3(X )U3 =
2: tJ.'j(x )uj
(3.22)
}=1
where 1jJJ are the quadratic Lagrange interpolation functions, (3.23)
SECONDORDER BOUNDARY VALUE PROBLEMS
83
Here De denotes the determinant of the matrix in (3.20b), and ai, Pi, and ri are defined by (3.21). The subscripts used in (3.21) permute in a natural order: if i = 1 then j = 2 and k = 3
For example, a~,
P;,
if i = 2 then j
=3
and k = 1
if i = 3 then j
=1
and k = 2
rt are given by xt(xy2, P; = (Xr)2 
(3.24)
and
ai = x3(xiY 
(xi?,
r~ =
x;  xi
The quadratic interpolation functions can be expressed in terms of a local coordinate X, with origin fixed at node 1, the left end of the element. The global coordinate x is related to the local coordinate .i by the relation
x =xi +x
(3.25)
where xi = XA is the global coordinate of the first node of the element Q,e. For a quadratic element with the interior node, node 2, located at x = ah e , we have
e() 1 X( .i) 1/J2 X = a(l a) h 1 h e _ 1/Jix)=
(3.26)
a.i ( 11 X)
(1a)h
ah
where 0 < a' < 1 and xi = xi + oh.: For a' = !, i.e., when node 2 is placed at the midpoint of the element, (3.26) becomes
(3.27)
Plots of the quadratic interpolation functions are given in Fig. 3.5(b). The function 1/Ji is equal to 1 at node i and zero at the other two nodes, but varies quadratically between the nodes. All Lagrange family of interpolation functions satisfy the following
84
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
properties, known as the interpolation properties: (1)
o
s, == {1
1jJ'[(xJ) =
if i =1= • 'f' J.
I =: J
1
n
(2)
2: lPj(X)
=:
1,
hence
1'=1
n
d 1V
1'=1
dx
(3.28)
L ' =0
where n 1 is the degree of interpolation polynomials and xJ is the global coordinate of node j in the element Qe. It can be verified that the linear interpolation functions in (3.16) and quadratic interpolation functions in (3.26) and (3.27) satisfy the two properties in (3.28). The first is a direct result of the requirement Ve(x'j) = u'j, and the second comes from the inclusion of a constant term in the polynomial. For example, if the approximation tr is to represent a uniform state of solution, U" = V'O = constant, then all uj = V o, and we have n
in =
2: U'Q1JJJ(x)
1'=1
or 1=
L" lPj(x) j=l
The interpolation properties (3.28) can be used to construct the Lagrange interpolation functions of any degree. For example, the quadratic interpolation functions (3.27) can be derived using property (1) of (3.28). Since 1JJHi) must vanish at nodes 2 and 3, i.e., at i = ih e and x = he> it is of the form
The constant C is to be determined such that 1 = ceO
1he)(O
he)
1JJl is equal to 1 at i = 0:
or
C = 2/h;
This gives
which is the same as in (3.27). The other two interpolation functions can be derived in a similar manner. Although a detailed discussion is presented here on how to construct the Lagrange interpolation functions for onedimensional elements, they are readily available in books on numerical analysis, and their derivation is independent of the physics of the problem to be solved. Their derivation depends only on the geometry of the element and the number and location of the nodes. The number of nodes must be equal to the number of terms in the polynomial. Thus, the interpolation functions derived above are useful not only in the finite element approximation of the problem at hand, but also in all
SECONDORDER BOUNDARY VALUE PROBLEMS
85
problems that admit Lagrange interpolation of the variables, i.e., all problems for which the primary variables are the dependent unknowns of the governing equations. STEP 3: FINITE ELEMENT MODEL. The weak form (3.8) or (3.10) is
equivalent to the differential equation (3.1) over the element Qe and also contains the natural boundary conditions (3.7). Further, the finite element approximations (3.16a) or (3.22) satisfy the essential boundary conditions (3.14) of the element. The substitution of (3.16a) or (3.22) into (3.8) will give the necessary algebraic equations among the nodal values u1 and Q1 of the element Qe. In order to formulate the finite element model based on the weak form (3.8), it is not necessary to decide a priori the degree of approximation of U". The model can be developed for an arbitrary degree of interpolation: 11
U= U" =
2: uj1J1j(x)
(3.29)
j=1
where 1J17 are the Lagrange interpolation functions of degree n  1. When n > 2, the weak form in (3.8) must be modified to include nonzero secondary variables, if any, at interior nodes: 0=
f
XB (
XA
du adw  d +CWU ) dx dx x
L wqdx ~ w(x1)Q1 X
8
"
(3.30)
,=1
XA
where x1 is the global coordinate of the ith node of element ge. If nodes 1 and n denote the endpoints of the element then Qi and Q~ represent the unknown
point sources, and all other Q1 are always known (i.e., applied point sources). Following the RayleighRitz procedure developed in Section 2.4.2, we substitute (3.29) for u and wi, 1J12' ... , 1J1~ for w into the weak form (3.30) to obtain n algebraic equations: 0=
( d1J1e) (" ) t("'8 [ a d1J1e dx j~ U} d/ + c1J1i ~1 U7tPj(X) 1
11
f)
("
)
] 11 tPiq dx  j~ 1J1i(xj)Qj
]
1J . 0= L: oy, oxy) are the' total stress components (in N m") and (fx,fy) are the x and y components of the body force vector (measured per unit volume). The total stress components can be expressed in terms of the viscous stress components (.x, .y, .Xy) and the hydrostatic pressure P (in N m Z) : Ox
= l"x P,
(3.104)
120
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
The viscous components of stress are related to the velocity gradients by Newton's law of viscosity. For isotropic, Newtonian fluids, these are
au
'ry =
21t ay'
'r
xy
=It
au) (au ay+ax
(3.105)
where It is the viscosity (in kg m" S1) of the fluid. Combining (3.103)(3.105), we obtain
au P ) +a [It (au au)] a (21t+
(au au au) ax ax ay ay ax +fx=p +u+vat ax ay (3.106) a [ (au av) J a ( av) (av av av) +f,Y=p +u+vax It +ay ax +ay 21tP ay at ax ay" /
These are known as the NavierStokes equations. When the viscosity is zero, we have
ap (au au au) +fx=p +u+vax at ax ay sr =p (av au au) +f, ay Y +u+vat ax ay
(3.107)
which are known as the Euler equations. Conservation of energy for incompressible fluids is expressed by
aT et aj =k (efT efn pc(+u+vy2 +q+~ at ax ay ax 2+a
(3.108)
where c is the mean heat capacity at constant volume, q is the internal heat generation, k is the thermal conductivity of the (isotropic) fluid, and I, v = 1>~, v = cP~, and v = 1>'4, to obtain a set of four algebraic equations. The ith algebraic equation of the finite element model is (for v = 1>1) 4 o_ /=1 ~
(1
X /1 dx
_
e
Qj
(4. 12a)
or 4
L Kijuj 
Fi = 0
(4. 12b)
/=1
where (4. 12c)
Note that the coefficients (4.12b) can be written as
Kit
are symmetric:
KZ· = KJi'
In matrix notation,
(4.13)
This represents the finite element model of (4.1). Here [K e ] is the stiffness matrix and {P} is the force vector of a beam element. When a transverse point force F~ is applied at a point Xo inside the element, it is distributed to the element nodes by the relation [see Remark 5 in Chapter 3: (3.60)]:
Fi = Fo1>'l(xo)
(4.14)
which contains both transverse forces (F 1 and F3) and bending moments (F; and F'4).
BENDING OF BEAMS
151
For the case in which b (= EI) and f are constant over an element, the element stiffness matrix [KeJ and force vector {F"} have the specific forms (see Fig. 4.2 for the element displacement and force degrees of freedom)
6  3h h3 6 [ 3h
[ "l = 2b K
sn 6 2 2h 3h hZ
3h 6 3h
3h] hz 3h 2h Z
(b
= EI = constant) (4.15)
It can be verified that the generalized force vector in (4.15) represents the
"statically equivalent" forces and moments at nodes 1 and 2 due to the uniformly distributed load over the element. For any given function f, (4.12c) provides a straightforward way of computing the components of the generalized force vector {r}.
4.2.4 Assembly of Element Equations The assembly procedure for beam elements is the same as that used for bar elements, except that we must take into account the two degrees of freedom at each node. Recall that the assembly of elements is based on: (a) interelement continuity of the primary variables (deflection and slope) and (b) interelement equilibrium of the secondary variables (shear force and bending moment) at the nodes common to elements. To demonstrate the assembly procedure, we select a twoelement model (see Fig. 4.4). There are three global nodes and a total of six global generalized displacements and six generalized forces in the problem. The continuity of the primary variables implies the following relation' between the element degrees of freedom Ut and the global degrees of freedom U; (see Fig. 4.4): u~
= VI, ui = Uz , uj = ui = U3 u~ = u~ = V4 , u~ = Us, u~ = U6
(4.16)
In general, the equilibrium of the generalized forces at a node between two connecting elements Qe and Q! requires that .
Q3 + Q{ =' applied external point force
Q4 + Q~ = applied external bending moment
(4.17)
If no external applied forces are given, the sum should be equated to zero. In
equating the sums to the applied generalized forces (i.e., force or moment), the sign convention for the element force degrees of freedom (see Fig. 4.2)
152
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
t
VI '" 0,
V2
• .
~
Element 1
Ql
VS. Q~ = Fo
V6,
Ql '" Mo
U~ut
:+;lzd QI
:r
Qt = unknown = 0, Qi = unknown
4
Ql,+ Q1
w(x) 1.2
0.8
0.4
FIGURE 4.4 Assembly of two EulerBernoulli (or classical) beam finite elements, and the finite element solution as a linear combination of the nodal values and interpolation functions.
should be followed. Forces are taken positive acting upward and moments are taken positive acting clockwise. To impose the equilibrium of forces in (4.17), it is necessary to add the third and fourth equations (corresponding to the second node) of ge to the first and second equations (corresponding to the first node) of Qt. Consequently, the global stiffnesses K 33 , K 34 , K 43 , and K44 associated with global node 2 are the superpositions of the element stiffnesses: K 44 = K1 + Kb (4.18)
In general, the assembled stiffness matrix and force vector for beam elements
BENDING OF BEAMS
153
connected in series' have the following forms:
[K]=
Global· nodes 1
2
~
"""
3
...._..
K~I
K~2
K1I
Kh
Kjl
Kj2
Kj3 + Kil
K1+Ki2
Kl i
Kl2
Kl 3 + K~I
Kit +K~2
(4. 19a)
F~
F1
{F} =
Fj+Fi F~+F~
(4. 19b)
F~
F~
The connectivity matrix [B] (which will be used in computer implementation) for the twoelement mesh is
[B]=[~ ~] Since there are two primary degrees of freedom per node, repetition of a number in [B] indicates that the coefficients associated with both degrees of freedom will add up. For example, the repetition of the global node number 2 (which corresponds to global degrees of freedom 3 and 4) in rows 1 and 2 indicates that the global K))> K)4, K 43 , and K 44 have contributions from both elements 1 and 2. The assembled system of equations for a mesh of two elements with hI = h2 = ~L = h, and constant EI and f [hence, [W] and {r} are given by
154
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
(4.15)] is 0 3h 3h 0 6 2 z 0 3h 3h 2h 0 h 3h 3h 6 3h 3h 6+6 2E1 6 3 Z h 3h h 3h 3h 2hz+2h z 3h hZ 3h 0 3h 6 0 6 Z 0 s« 0 h 3h 2hz
U1
6
u, U3 U4 Us U6
Qt h Qi =fh 6+6 + ).?5 + Qi '" (4.20) 12 6
hh
Ql+Q~
6
Q~
h
Q~
The reader is cautioned that (4.20) does not' represent the assembled equations of any two beam elements; it is based on the assumption that hi = hz, (E1)1 = (£1)2, and (fh = (f)z. Equations (4.19) are more general.
4.2.5 Imposition of Boundary Conditions At this step of the analysis, we must specify the particular boundary conditions, i.e., geometric constraints and forces applied, of the particular problem to be analyzed. The type of essential (also known as geometric) boundary conditions for a specific beam problem depend on the nature of the geometric support. Table 4.1 contains a list of commonly used geometric supports for beams. For the sake of completeness, the boundary conditions on the axial displacement u are also included. The natural (also called force) boundary conditions involve the specification of generalized forces when the corresponding primary variables are not constrained. Here we consider a cantilever beam (i.e., a beam fixed at one end and free at the other) of length L, flexural rigidity £1 (= constant), and subjected to a uniformly distributed force fo, end force Po, and end moment Mo (see Fig. 4.4). First, we write the equilibrium conditions for the generalized forces. At global node 1, Qt and Qi (the shear force and the bending moment, respectively; i.e., the reactions at the fixed end) are not known. At global node 2, there are no externally applied shear forces and bending moment. Hence,
Qj+Qi=O,
Ql+Q~=o
(4.21a)
At global node 3, the shear force is given as 1'0, and the bending moment as Mo (note the sign convention for FO and Mo from Fig. 4.2):
Q~=_(bd2~:)1 dx
x=L
=Mo
(4.21b)
155
BENDING OF BEAMS
'fABLE 4.1
Types of commonly used support conditions for beams and frames
Type of support
zf E3~x'
Displacement boundary conditions
Force boundary conditions
None
All, as specified
u=o 11'=0
Moment is specified
FREE
zt
~x PINNED
q
~Rx
. u=o
Transverse force and moment are specified
(vertical)
px
11'=0
Horizontal force and bending moment are specified
u=O 11'=0 dw/dx=O
None specified
ROLLER (horizontal)
FIXED or CLAMPED
Next, we identify and impose the specified generalized displacements. Since the beam is clamped at global node 1, it follows that the deflection wand the slope dw/ dx are zero there:
ui == w~ = V
1
u~ ==
= 0,
01 = Vz = 0
(4.22)
Using (4.21) and (4.22) in (4.20), we obtain
6
3h
6
3h
0
0
2
Z
3h h 0 0 3h 2h , ·iz······6···:·6:..=.:3h :..=.6"3h·· 2EI h 3 3h h 2 o 4h Z 3h h 2
o
0
6
o
0
3h
3h hZ
6
3h
3h
2h Z
V 1=0
Qi
Vz=O
Q~
u;V4
foh 12
+
o o
V5
Po
V6
Mo
(4.23)
156
FINiTE ELEMEIIT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
4.2.6 Solution Equation (4.23) contains six algebraic equations in six unknowns: (QL Qi, U3 , U4 , Us, U6 ) . Because the algebraic equations for the unknown generalized displacements (U3 , U4 , Us, U6 ) do not contain the unknown generalized forces (QL QD, equations 3, 4, 5, and 6 of the system (4.23) can be solved independently; the known values of the displacements U, and Vz are used in equations 3 and 4. This provides us with the motivation to partition (shown by dashed lines) the matrix equation (4.23), which can be recast in the form (4.24) r
where {VI} denotes the column of known generalizeddisplacements, {U Z } the column of unknown generalized displacements, {pl} the column of unknown forces, and {pZ} the column of known forces. Equation (4.24) can be written, after carrying out the matrix multiplication, in the form [cf. (3.56)]
or (4.25a) (4.25b)
These are the condensed equations for the generalized displacements and forces, respectively. Since {U l } and {pZ} are known ({U l } = {On, we can use (4.25a) to solve for {V z}, and then use (4.25b) to compute the unknown reactions {Pi}:
6 3h]~l{ foh } 3h h 6 3h Po + YOh (4.26) z 3h 2h  + tzfoh2 2
0
Mo
Inverting the matrix (say, by Gaussian elimination) and performing the matrix
BENDING OF BEAMS
157
multiplication (i.e.;. solving the equations), we obtain Z
h {U } = 6EI Z
[
2h 3h 3h 6. 5h 2 9h
6
 3h
z 5h 9h 16h 2
12h
Fa + Yah
12h
12
 Mo + tzfoh 2
3h]{
foh
6
}
0
3
5h2Fa + 3hMo + ¥foh } h 9hFo  6Mo  7foh z = 6EI { 16h zFQ + 12hMo + 12foh3 12hFo  12Mo  8foh z
(4.27)
The reactions Ql and Qi can be obtained by substituting (4.27) into (4.25b). The Q'f obtained from the element (equilibrium) equations are more accurate than those obtained from the definitions (4.3), wherein the derivatives of ware obtained by differentiation of its finite element interpolation. The reactions from equilibrium (QDequil are
{
1}
F
=
{Q~} = 2EI 3
Qi
h
[63h
3h
h
2
0 0
i}_
°OJ{ Us
6}
foh { 12h
U6
. { (Fa+2foh) } = 2h(Fo + foh) + Mo
(4.28)
It can be verified that the reactions Q~ and Qi in (4.28) satisfy the static equilibrium equations of the beam:
Qi + (Fa + 2foh) = The reactions
0,
Qi
(2Foh + 2foh z + Mo)
=0
Qi and Qi can also be computed using the definitions Q 11 ==.!i dx
(E1 dZlzV) I, dx x=o
I
1 ( dZW) Qz== EI dx 2 x=o
(4.3): (4. 29a)
From (4.10b), we note that the second derivative of the Hermite cubic interpolation functions is linear over the element and the third derivative is constant over the element. Therefore, the reactions, i.e., bending moment and shear force, computed using the definition (4.3) are elementwise linear and constant, respectively. Further, at nodes connecting two elements, they yield discontinuous values because the second and third derivatives of ware not
158
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
made continuous across the interelement nodes. Substituting (4.10b) into (4.29a), we obtain the values
3ep1) Q}=EI ( U3 dJepl dX33+U4 ddx J4
I
r (12)
x 0
1
(6)] =(FQ+!foh)
=EI U3 h 3 +U4 h 2
=
e
(6 + U he2) = (Mo + 2Foh," + lifoh2)
Q2 = EI U3 h;
2.3
4
(4.29b)
e
'
which are in error by Ii = Yoh and 12 = rJOh 2 compared with those in (4.28).
4.2.7 Postprocessing of the Solution The finite element solution as a function of position x is given,by
e( )_{U34>~+U4ep~ for O~x~h wx2 2 2 2 U34>1 + U44> 2 + U54>3 + U6 4>4 for h ~x ~2h
(4.30a) .
where [see (4.10a)]
4>~=3Gf 2Gr epl=h[Gf GY]
4>1 =13(1 ~r  2(1~Y, 4>~ =h(1 ~)(2~r
(4.30b)
4>~=3(1~r+2(1~r 4>~=h[(l~r +(l~r] The exact solution of (4.1) subject to the boundary condition (4.21) can be obtained by direct integration, and is given by .
Elw(x) = fJ'ox  UFo +f oL)x + :!(Mo + FoL + YoL ) X EI£J(x) = tJoX 3 + HFo + foL )x 2  (Mo + FoL + VaLL,) 4
M(x) = Yox
2
2
3
(Fo + foL)x

2
}
for 0 ~ x ~ L
+ Mo + PoL + YoL
2
(4.31) The finite element solution for the data in Table
4.2
(4.30) and the exact solution (4.31) are compared
fo=24kNm 1 , 6kN
E = 200 x 10
m
Fo=60kN, 2
,
1 =29 X
L=3m, Mo=OkNm 106mm4 (El = nr')
5800kN
As expected, the finite element solution for wand (J coincides with the exact solution at the nodes. At other points, the difference between the finite element and exact solutions is less than 2%. This completes the finite element formulation and analysis of the fourthorder differential equation (4.1) governing the EulerBernoulli beam theory. Whenever the flexural rigidity b == El is a constant in each element, the
BENDING OF BEAMS
159
TABLE 4.2
Comparison of the finite element solution with the exact solution of the cantilever beam of Fig. 4.4 (2 elements; L = 3 m, EI = 5800kNm2,Jo=24kNm~t, lb=6(}kN, Mo=OkNm) w(m)
(J
=dw{dx
M{EI (kN m)
x [rn)
FEM
Exact
FEM
Exact
FEM
Exact
0.00 0.1875 0.3750 0.5625 0.7500 0.9375 1.1250 1.3125 1.5000 1.6875 1.8750 2.0625 2.2500 2.4375 2.6250 2.8025
O.OOOOt 0.0008 0.0033 0.0071 0.0124 0.0188 0.0263 0.0347 0.0439t 0.0539 0.0644 0.0754 0.0868 0.0986 0.1106 0.1228 0.1350t
0.ססOO
O.OOOOt 0.0088 0.0169 0.0244 0.0311 0.0372 0.0426 0.0472 0.0512t 0.0546 0.0575 0.0600 0.0620 0.0635 0.0645 0.0651 0.0652t
0.ססOO
0.0008 0.0033 0.0072 0.0124 0.0188 0.0263 0.0347 0.0439 0.0539 0.0644 0.0755 0.0869 0.0987 0.1107 0.1228 0.1350
0.0089 0.0171 0.0245 0.0311 0.0371 0.0425 0.0471 0.0512 0.0547 0.0576 0.0600 0.0620 0.0634 0.0644 0.0650 0.0652
0.0489t 0.0452 0.0415 0.0378 0.0341 0.D305 0.0268 0.0231 0.0194t 0.0169 0.0144 0.0118 0.0093 0.0068 0.0043 0.0017 O.0008t
0.0497 0.0455 0.0414 0.0375 0.0338 0.0301 0.0266 0.0234 0.0202 0.0171 0.0143 0.0115 0.0089 0.0065 0.0042 0.0020
3.ססOO
0.ססOO
t Nodal values; all others are computed by interpolation.
element stiffness matrix (4.15) can be used directly. The finite element solution for the generalized displacements at the nodes is exact for all problems for which E1 is constant and the applied transverse load f is a polynomial expression. Further, the solution is exact at all points if the distributed load is such that the exact solution is a cubic. The bending moment at any point in the element ge of the beam can be computed from the formula (4. 32a)
For beams of rectangular crosssection with height H (and width B) the maximum bending stress is (4.32b)
The minus sign is for the top and the plus sign for the bottom of the beam.
160
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
. .= 2400Ib re'
10,000lb
EI
=
O
I
(5.15)
where c is a constant. This estimate implies that the error goes to zero as the pth power of h as h is decreased (or the number of elements is increased). In other words, the logarithm of the error in the energy norm versus the logarithm of h is a straight line whose slope is k + 1  m. The greater the degree of the interpolation functions, the more rapid the rate of convergence. Note also that the error in the energy goes to zero at the rate of k + 1  m; the error in the ~ norm will decrease even more rapidly, namely, at the rate of k + 1, i.e., derivatives converge more slowly than the solution itself. Error estimates of the type in (5.15) are very useful because they give an idea of the accuracy of the approximate solution, whether or not we know the true solution. While the estimate gives an idea of how rapidly the finite element solution converges to the true solution, it does not tell us when to stop refining the mesh. This decision rests with the analysts, because only they know what a reasonable tolerance is for the problems they are solving. As an example of estimating the error in the approximation, i.e, (5.15), consider the linear (twonode) element for a secondorder equation (m = 1). We have for an element . u,,=ul(1s)+u2s
(5.16)
where s = x/h and x is the local coordinate. Since Uz can be viewed as a function of Ul via (5.16), one can expand in a Taylor series around node 1 to obtain
U2
(5.17)
l
204
FINITE ELEMENT ANALYSIS OF ONE·DlMENSIONAL PROBLEMS
where u'
~dulds.
Substituting this into (5.16), we obtain (5.18)
Expanding the true solution in a Taylor series about node 1, we obtain (5.19) Therefore, we have, from (5.18) and (5.19),
IUh  ul ~ 21 (s  s 2) max
0",, 0) or divergence (p < 0).
REFERENCES FOR ADDITIONAL READING Ciarlet, P. G.: The Finite Element Method for Elliptic Problems, NorthHolland, Amsterdam, 1978.
Davies, A. L.: The Finite Element Method, A FirstApproach, Clarendon Press, Oxford, 1980. Oden, J. T., and J. N. Reddy: An Introduction to the Mathematical Theory of Finite Elements, WileyInterscience, New York, 1982. Reddy, J. N.: Applied FUllctionalAnalysis and Variational Methods in Engineering, McGrawHili, New York, 1986; Krieger, Melbourne, FL, 1991. . Rektorys, K: Variational Methods in Mathematics, Science and Engineering, Reidel, Boston, 1977.
Strang, G., and G. J. Fix: An Analysis of the Finite Element Method, PrenticeHall, Englewood Cliffs, NJ, 1973. Wait, R., and A. R. Mitchell: Finite Element Analysis and Applications, John Wiley, New York, 1985.
CHAPTER
6 EIGENVALUE AND TIMEDEPENDENT PROBLEMS
6.1 EIGENVALUE PROBLEMS 6.1.1 Introduction Determination of the values of the parameter Asuch that the equation A(u) = AB(u)
(6.1)
where A and B denote linear differential operators, has nontrivial solutions u is called an eigenvalue problem. The values of A are called eigenvalues and the associated functions U are called eigenfunctions. For example, the equation d 2u  dx'1= AU,
d2
with A = dx 2'
B =1
which arises in connection with the axial oscillations of a bar or the transverse oscillations of a cable, constitutes an eigenvalue problem. Here A denotes the square of the frequency of vibration to, In general, the determination of the eigenvalues is of engineering as well as mathematical importance. In structural problems, the eigenvalues denote either natural frequencies or buckling loads. In fluid mechanics and heat transfer, eigenvalue problems arise in connection with the determination of the homogeneous parts of the solution. In these cases, eigenvalues often denote 209
210
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROllLEMS
amplitudes of the Fourier components making up the solution. Eigenvalues are also useful in determining the stability characteristics of temporal schemes, as discussed in Section 6.2. In this section, we develop finite element models of eigenvalue problems. In view of the close similarity between the equations of eigenvalue and boundary value problems, the steps involved in the construction of their finite element models are entirely analogous. Differential eigenvalue problems are reduced to algebraic eigenvalue problems by means of the finite element approximation. The methods of solution of algebraic eigenvalue problems are then used to solve for the eigenvalues and eigenvectors.
6.1.2 Formulation of Eigenvalue Problems Consider the parabolic partial differential equation
au a(kAau) = q(x, t)
pcA~
at ax
ax
(6.2)
which governs transient heat transfer in onedimensional systems (e.g., a plane wall). Here u denotes the temperature, k the thermal conductivity, p the density, A the crosssectional area, c the specific heat, and q the heat generation per unit length. The homogeneous solution (i.e., the solution when q = 0) of (6.2) is often sought in the form of a product of a function of x and a function of t (i.e., through the separation of variables technique):
u(x, t) = S(x)T(t)
(6.3)
Substitution of this assumed form of solution into the homogeneous form of (6.2) gives dT d ( kAdS) T=O peASdt
dx
dx
Separating the variables (assuming that peA and kA" are functions of x only), 1 dT
1 1d (
"Tdi= peASdx
dS)
kA dx
(6.4)
Note that the lefthand side of this equation is a function of t only while the righthand side is a function of x only. For two functions of two independent variables to always be equal, both must be equal to the same constant, say ,\:
. .!: dT =_1_!~ (kA dS) =_). T dt
peA Sdx
dx
(6.5a)
or dT dt
=}..T
'
d ( kAdS) ).pcAS=O dx dx
(6.5b)
EIGENVALUE AND TlME·DEPENDENT PROBLEMS
211
The negative sign In the constant A is based on the physical requirement that the solution SeX) be harmonic in x while T(/) must decay exponentially with increasing I. The solution of the first equation is T = 'ToeAt. When k, A, p, and c are constants, the solution of the second equation is
Sex) = B 1 sin VJ.x
+ B z cos YXx, }..
The constants A, To, B 1 , and B 2 are determined with the help of initial and boundary conditions. In view of the above discussion, the solution of (6.2) is of the form
u(x, I)
=
U(x)e M
(6.6)
This form is consistent with the solutions we derived above, with U(x) ~ Sex) 'To. Substituting (6.6) into the homogeneous form of (6.2), we obtain
d ( kA~ dU) eA'ApcAU(x)eA'=O dx dx or
d ( kA~ dU) ApcAU=O dx dx
(6.7)
We wish to determine). and nonzero Vex) such that (6.7) holds and the boundary conditions of the problem are met. Equation (6.7) describes an eigenvalue problem, A being the eigenvalue and U(x) the eigenfunction. The axial motion of a bar can be described by the hyperbolic equation
cru a( au)
p AZ   EA =q(x, I)
at
ax
ax
(6.8)
Here u denotes the axial displacement, E the modulus of elasticity, A the crosssectional area, p the density, and q the axial force per unit length. The natural axial oscillations of the bar are periodic, and they can be determined by assuming a solution of the form
uix, I) = U(x)e;wt, with i = v=I
(6.9)
where to denotes the frequency of natural axial motion (or vibration), and U(x) denotes the configuration of the bar, called the mode shape, during the vibration. For each value of oi, there is an associated mode shape. Substitution of (6.9) into the homogeneous form of (6.8) gives
[PAWZU 
~ (EA ~~) ]e;wl= 0
2U
FINITE ELEMENT ANALYSIS OF ONE·DIMENSIONAL PROBLEMS
or (6.10) where A = (1)2. Equation (6.10) is an eigenvalue problem, which involves determining the square of natural frequencies 1\. and mode shapes U. This equation also arises in the solution of (6.8) by means of separation of variables, as discussed in connection with the parabolic equation. Note that (6.7) and (6.10) are of the same form. Only the coefficients are different. Equations similar to (6.10) can be derived for the transverse vibrations of a beam using the EulerBernoulli or the Timoshenko beam theories. For the EulerBernoulli beam theory, we assume
\V(x, t) = W(x)eia>t where (I) is the frequency of natural transverse motion and W(x) is the mode shape of the transverse motion. Substitution of this form into the equation of motion of the EulerBernoulli beam theory
~w
pA at2
f? (
+ ax2
~W)
E1 ax 2 = 0
(6.11)
gives (6.12) where 1\.= w2 • For the Timoshenko beam theory, we assume
w(x, t) = W(x)eia>t,
W(x, t) = S(x)eia>t
and substitute into the equations of motion of the theory 2
a w2  a [GAK (aw pA +'P)] =0 at ax ax
a(aw)
(
(6.13)
aw ) =0 pOZw 12   £1 +GAK +'P at ax ax ax to obtain the eigenvalue problem

~ [ GAK( ~: + S) ] 
1\.pA W = 0
d(dS) £1 +GAK(dW) +S 1\.p1S=0 dx dx dx
(6.14)
The study of buckling of beamcolumns also leads to an eigenvalue problem. For example, the equation_ governing the equilibrium of a beam
EIGENVALUE AND TIME·DEPENDENT PROBLEMS
213
subjected to an axial force P, according to the EulerBernoulli beam theory, is
2
2
. d ( d 1V) d'w dx2 E~ dx 2 + P dx 2 = 0
(6.15)
which is an eigenvalue equation with A = P as the eigenvalue, which represents the buckling load. Often one is interested only in the smallest value of A, called the critical buckling load. For the Timoshenko beam theory, the buckling equations are d ~
(dW
)] +P=O d w
[ GAK +W dx dx
2
dx?
d(dW) EI +GAK (dW) +W =0 dx dx dx
(6.16)
~
This completes the formulation of eigenvalue problems associated with the model problems studied in this book. In the next section, we develop the finite element models of (6.7), (6.10), (6.12), and (6.14)(6.16).
6.1.3 Finite Element Models An examination of the eigenvalue equations derived in the previous section shows that they are a special case of the equations studied in Chapters 3 and 4. For example, consider
_~(adU) =f dx
dx
(6.17)
The eigenvalue equation associated with this is
 dx~ (a dU) = ACo U dx where a and transfer,
Co are
(6.18)
quantities that depend on the physical problem: for heat a = kA,
Co= pcA
where c is the specific heat, while, for a bar,
a = EA,
co=pA
Similarly, the eigenvalue equations associated with the transverse vibrations and buckling of beams are special cases of their static counterparts. Therefore, the finite element models of the eigenvalue equations can readily be developed. It is important to note that the spatial derivative operators of the static (i.e., nontimedependent) and eigenvalue equations are the same. The difference between (6.17) and (6.18) is that, in place of the source term f, we have ACOU in the eigenvalue equations. This difference is responsible for
214
FINITE ELEMENT ANALYSIS OF ONE·D1MENSIONAL PROBLEMS
another coefficient matrix, in addition to the usual coefficient matrix [K e ] , in the eigenvalue problems. The derivation of the finite element models of eigenvalue equations is presented next. Over a typical element Qe, we seek a finite element approximation of u in the form n
U=
2: uilJi/(x)
(6.19)
j=l
(6.20a) where w is the weight function, and Ql and Qn are the "usual secondary variables (Qi = 0, i =F 1, n) (6.20b)
Substitution of the finite element approximation into the weak form gives the finite element model of the eigenvalue equation (6.18): (6.21a)
where
«: _ fX 8 d1J!i d1/JJ d e:
a dx
XA
d
X
(6.21b)
X,
Equation (6.21a) represents the finite element model of the eigenvalue equations (6.18) and (6.10). The finite element model of (6.12) is (6.22a)
where {u e} and {Qe} are the columns of nodal generalized displacement and force degrees of freedom of the EulerBernoulli beam element:
L~(EI~~)l
W'l
(~: {u"} =
W2
(~:
{Qe} =
d2W) (El dx 2 1
L
[ ~ (EI~x~) 2W) ( EI ddx 2 2
(6.22b)
where the subscripts 1 and 2 refer to element nodes 1 and 2 (at x = XA and x = XB)' The matrices [Ke ] and [Me], known as the stiffness and mass matrices,
EIGENVALUE AND TIME·DEPENDENf PROBLEMS
are defined by Kij=
d2ifl7 d 2iflJ EJ2 d2 dx , dx X
XB
f
xA
Mij=
fX8
pA 0 (i.e., t > 0)
problem require
while the initial condition requires
Since the initial condition should be consistent with the boundary conditions, we take (U1)0 = 0.0. Using the boundary conditions, we can write for the oneelement model (h = 1.0) (6.59)
EIGENVALUE AND TIME·DEPENDENT PROBLEMS
235
which can be solved repeatedly for U2 at different times, s = 0, 1, .... Repeated use of (6.59) can cause the temporal approximation error to grow with time, depending on the value of a. As noted earlier, the forward difference scheme (a = 0) is a conditionally stable scheme. The critical time step is given by
At
er
where
Amax
== 2/}'max
is the maximum eigenvalue associated with (6.57c):
A{M]{U} + [K]{U} == {O} For the model at hand, this reduces to ~AhU2
+ h 1 U2 =
0, or A= 31h2= 3
Hence Mer= 0.6667. Thus, in order for the forward difference solution of (6.59) to be stable, the time step should be smaller than M« == 0.6667; otherwise, the solution will be unstable, as shown in Fig. 6.7. For unconditionally stable schemes (a ~ D, there is no restriction on the time step. However, to obtain a sufficiently accurate solution, the time step must be taken as a fraction of !:lt er • Of course, the accuracy of the solution also depends on the mesh size h. As this is decreased (i.e., the number of elements is increased), !:lt er decreases. Figure 6.8 shows plots of u(l, t) versus time for a = 0.5 and M = 0.05. Solutions predicted by meshes of one and two linear elements and the mesh of one quadratic element are compared with the exact solution. The convergence of the solution with increasing number of elements is clear. The finite element solutions obtained with different methods, time steps, and meshes are compared with the exact solution in Table 6.3.
1
I/(l, t)
o
1 ilt = 0.675, lL ex = 0.0 ~
2 +.,...,r.r..+ 6 8 4 2 o Time t
FIGURE 6.7 Transient solution of a parabolic equation according to the forward difference scheme (.t'l.t = 0.675, one linear element). The solution is unstable because the time step is larger than the critical time step.
236
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
1.2 +',.L_~_l.__'__...L_~__'__
__'__+
..._.__._. Analytical 1.0
0
1L
~2L
 lQ
0.8
tJ.t = 0.05
u(l, t)
" = 0.5
0.6
0.4 0.2 /
0.0 +.......r..,......,,...,,..+ 0.6 0.2 0.0 0.8 0.4 1.0 Time t FIGURE 6.8 Transient solution of a parabolic equation according to linear and quadratic finite elements. TABLE 6.3
A comparison of the finite element solutions obtained using various time approximation schemes and meshes with the analytical solution of a parabolic equation arising in conductive heat transfer 1L Time t a=O 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0040 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
«=1
Linear elements (« = 0.5) lL
2L
4L
8L
Quadratic elements ( n
(7.15)
In subparametric formulations, the geometry is represented by lowerorder elements than those used to approximate the dependent variables. An example of this category is provided by the EulerBernoulli beam element, where the Hermite cubic element is used to approximate the transverse deflection, while the geometry is approximated, when straight beams are analyzed, with linear interpolation functions. In isoparametric formulations (the most common in practice), the same element is used to approximate the geometry as well as the dependent unknowns: 'lJli(X) = ~i(~)' In the superparametric formulations, the geometry is represented with higherorder elements than those used to approximate the dependent variables. This formulation is seldom used in practice.
7.1.5 Numerical Integration The evaluation of integrals of the form
f
F(x) dx
(7.16)
by exact means is either difficult or impossible owing to the complicated form of the integrand F. Numerical integration is also required when the integrand is to be evaluated inexactly (as in the Timoshenko beam element) and when
252
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
the integrand depends on a quantity that is known only at discrete points (e.g., in nonlinear problems). The basic idea behind all numerical integration techniques is to find a function P(x) that is both a suitable approximation of F(x) 'and simple to integrate. The interpolating polynomials of degree n, denoted by PlI , which interpolate the integrand at n + 1 points of the interval [a, b], often produce a suitable approximation and possess the desired property of simple integrability. An illustration of the approximation of the function F(x) by the polynomial Pix) that exactly matches the function F(x) at the indicated base points is given in Fig. 7A(a). The exact value of (7.16) is given by the area under the solid curve, while the approximate value
f ...
Pix)dx
/
'"'"
F(x)
L.JL...._'_ _'_ _'_ _  '  _ _
a =
XJ
X
X3
(a)
''...L.x Xz
(b)
_F(x)
''''x Xz (c)
FIGURE 7.4 Numerical integration by the NewtonCotes quadrature: (a) ap· proximation of a function by P4(X); (b) the trapezoidal rule; (c) Simpson's rule.
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
253
is given by the area under the dashed curve. It should be noted that the difference (i.e., the error in the approximation) E = F(x)  P4(X) is not always of the same sign, and therefore the overall integration error may be small (because positive errors in one part cancel negative errors in other parts), even when P4 is not a good approximation of F. The commonly used integration methods can be classified into two basic groups: 1. The NewtonCotes formulae that employ values of the function at equally spaced base points. 2. The Gauss quadrature formula that employs unequally spaced base points. These are described here. THE NEWTONCOTES QUADRATURE. For r equally spaced base points, the NewtonCotes integration formula is given by
(b
r
)(1 F(x) dx = (b  a) ~l F(x[)w[
(7.17)
where WI are the weighting coefficients and x, are the base points, which are equally spaced. For r = 1, (7.17) gives the rectangle formula. For r = 2, it gives the familiar trapezoidal rule, in which the required area under the solid curve in Fig. 7.4(b) is approximated by the area under the dotted straight line [i.e., F(x) is approximated by Pl(X)]:
i~~~2 F(x) dx = !h[F(Xl) + F(X2)],
E = O(h
3)
(7.18)
where E denotes the error in the approximation and h is the spacing between two base points. The notation O(h), read as "order of h", is used to indicate the order of the error in terms of the spacing h. For r = 3, (7.17) gives the familiar Simpson's onethird rule (see Fig. 7.4c):
i~:~3 F(x) dx = !h[F(Xl) + 4F(X2) + F(X3)],
E = O(h5 )
(7.19)
The weighting coefficients for r = 1, 2, ... , 7 are given in Table 7.1. Note that I;~!fWI = 1. The base point location for r = 1 is Xl = a +!(b  a) = !(a + b). For r> 1, the base point locations are al
= a,
a2 = a
+ ~x,
""
Or
~a
+ (r
 1) ~x = b
where ~ = (b  a)/(r 1).A comment is in order on the use of the NewtonCotes integration formula (7.17). For mn r (i.e., when there is an even number of intervals or an odd number of base points), the formula is exact when F(x) is a polynomial of degree r + 1 or less; for odd r, the formula is exact for a polynomial of degree r or less. Conversely, a pthorder polynomial is integrated exactly by choosing r = p + 1 base points.
254
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
TABLE 7.1
Weighting coefficients for the NewtonCotes formula (7.17) r
WI
1 2 3
1
Wz
W4
W3
Ws
W6
w7
1 1
2
Z
t
'!
1
J.
8
5
.1
3
1
90
90
ik
2Q.
50
288
288
216
.zi.
272 840
90
Z8B
.31840
7
6
8 1£
8 3Z
90 .is,
6
1
6
4
8 32
840
840
.1
90 .J2.. 288
sa: 840
19 2118
ill 840
41
840
THE GAUSSLEGENDRE QUADRATURE. In the NewtonCotes quadrature, the base point locations have been specified. If the XI are not specified then
there will be 2r + 2 undetermined parameters, the weights )Vl and base points which define a polynomial of degree 2r + 1. The GaussLegendre quadrature is based on this idea. The base points XI and the weights WI are chosen so that the sum of the r + 1 appropriately weighted values of the function yields the integral exactly when F(x) is a polynomial of degree 2r + 1 or less. The GaussLegendre quadrature formula is given by (see Fig. 7.5)
Xl>
(7.20) where
WI
~
1
§
are the weight factors,
~~
0.57735
~ =
(a)
(b)
are the base points [roots of the Legendre
§= 1
0.57735
§= 1
1
0.;7~597
~=
~l
§= 0.774597
FIGURE 7.5 The twopoint (a) and threepoint (b) GaussLegendre quadratures.
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
255
TABLE 7.2
Weights and Gauss Legendre quadrature
Points ~i
points for the Gauss
T
Weights
0.000000000o
Onepoint formula
2.OOOOOOOOOO
±O.5773502692
Twopoint formula
1.OOOOOOOOOO
O. OOOOOOOOOO ±0.7745966692
Threepoint formula
0.8888888889 0.5555555555
±0.3399810435 ±0.8611363116
Fourpoint formula
0.6521451548 0.3478548451
O. ()()()(J{)(J()( ±0.5384693101 ±0.9061798459
Fivepoint formula
0.5688888889 0.4786286705 0.2369268850
±0.2386191861 ±0.6612093865 ±0.9324695142
Sixpoint formula
0.4679139346 03607615730 0.1713244924
polynomial Pr+l(;)], and
P is the
lVi
transformed integrand
PC;) =
F(x(;»J"(';)
(7.21)
The weight factors and Gauss points for the GaussLegendre quadrature (7.20) are given, for r = 1, ... , 6, in Table 7.2. The GaussLegendre quadrature is more frequently used than the NewtonCotes quadrature because it requires fewer base points (hence a saving in computation) to achieve the same accuracy. The error in the approximation is zero if the (2r + 2)th derivative of the integrand vanishes. In other words, a polynomial of degree p is integrated exactly by employing r =!(p + 1) Gauss points. When p + 1 is. odd, one should pick the nearest larger integer: r = [!(p
+ 1)]
(7.22)
In finite element formulations, we encounter integrals whose integrands F are functions of x, 1jJ;(x) and derivatives of 1jJ;(x) with respect to x. For the GaussLegendre quadrature, we must transform F(x) dx to p(.;) d; to use (7.20). For example, consider the integral
( ) d1jJ'f d1jJJ dX K ~. = fXB ax dx If
XA
dx
(7.23)
Using the chain rule of differentiation, we have (7.24)
256
FINITE ELEMENT ANALYSIS OF ONE·DIMENSIONAL PROBLEMS
Therefore, the integral can be written, with the help of (7.10), as I
Kij =
J
I
1 dljJ~ 1 dljJ:
a(x(~)) j d~' j d~J J d~
(7.25)
r
2: Fij(Sf)Wf
=
(7.26)
f~1
where ~e
F jj
1 dljJ1 d1jJj
_
aj d~

m
d{;'
j =
d,p~
~ xf d{;'
(7.27)
For the isoparametric formulation, we take 1jJf = ,pf. The transformation from S is not required for the NewtonCotes quadrature. "As noted earlier, the Jacobian matrix will be the same {Je = ~he) when the element is straight, even if the coordinate transformation is quadratic or cubic. However, when the element is curved, the Jacobian is a function of S for transformations other than linear. It is possible to determine the number of Gauss points required to evaluate the finite element matrices
x to
e
K jj =
fXB d1jJf d1jJj ~d d dx, XA x X
Mij=
fXB
1jJfljJJ dx,
Fi =
XA
1X8
1jJi dx
(7.28)
XA
exactly using linear, quadratic, and cubic interpolation functions. For linear: interpolation functions, the integrand of Kij is constant, requiring only onepoint Gauss quadrature. The integrand of the mass matrix Mij. is quadratic (p = 2), requiring [r = !(p + )1 =~] twopoint quadrature. The coefficients are evaluated exactly by onepoint quadrature. Similarly, for quadratic and cubic elements, we can estimate the number of Gauss points needed to evaluate Kij, Mij, and fi exactly. The results are summarized below. Note that, in estimating the quadrature points, it is assumed that the Jacobian is a constant, which holds true when the element is a straight line.
n
Number of Gauss quadrature points needed Element type
Kif
Mij
f~
Linear Quadratic Cubic
1 2 3
2 3 4
2
1 3
. If the matrices in (7.28) have variable coefficients or the elements are curved [and hence Je = j"(';)], the degree of the variation of the integrands changes and the number of Gauss points needed to evaluate the integral exactly changes. If the elements are straight, and the coefficients a = a(x) and c = c(x) together with f = f(x) are no more than linear in x, then the number
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
257
of Gauss points for evaluating the coefficients
Kij =
"'B dlJJ~ d d1/J~ dx, j . a d X X I
J
"'A
Mij =
j"'B clJJiVJJ dx
(7.29)
XA
remain the same as listed in the above table. However, the evaluation of fi requires one point more than before. Conversely, the twopoint quadrature for linear elements, threepoint quadrature for quadratic elements, and fourpoint quadrature for cubic elements would exactly evaluate Kij with a quadratic variation of a(x), Mij with linear variation of c(x), and fi with quadratic variation of f(x). The use of Gauss quadrature on (7.28) yields the following values (exact up to the fifth decimal place) when the element is straight and the isoparametric formulation is used: Quadratic (threepoint formula)
1 [
[K] "i:
e
h [ [M] = 1~
0.33333] 2.33333 2.66667 5.33333 2.66667 2.66667 0.33333 2.66667 2.33333 1.33333 0.66667 0.33333] 0.66667 0.66667 5.33333 0.33333 0.66667 1.333433
(7.30)
0' 166667} {F} = he 0.666667 { 0.166667 Cubic (fourpoint formula)
1.350 0.325] 3,700 4.725 1.350 1 4.725 10.800 7.425 [K]=~ he [ 1.350 7.425 10.800 4.725 0.325 1.350 4.725 . 3.700 0.761905 0.589286 0.214286 0.113095] 0.589286 3.857143 0.482143 0.214286 = [ ] M 10 0.214286 0.482143 3.857143 0.589286 [ 0.113095 0.214286 0.589286 0.761905
he
{F} =
(7.31)
he{~:~~} 0.125
In Section 7.2, we study the computer implementation of the steps involved in the finite element analysis of onedimensional problems. As a part
258
FINITE ELEMENT ANALYSIS OF ONE·DIMENSIONAL PROBLEMS
of the element calculations there, the computer implementation of the numerical integration ideas presented in this section will be studied. A model finite element program (FEM1DV2) for the solution of onedimensional problems is also described, and its application is demonstrated via several examples. Appendix 1 contains a source listing of program FEM1DV2.
7.2
COMPUTER IMPLEMENTATION
7.2.1 Introductory Comments Chapters 36 were devoted to the finite element formulations of two classes of boundary value, initial value, or eigenvalue problems in one dimension: I'
"
1. Secondorder differential equations (e.g., heat transfer, fluid mechanics, ID elasticity, bars, and the Timoshenko beam theory);' 2. Fourthorder differential equations governing the EulerBernoulli beam theory. The frame element, obtained by superposing the bar element and the beam element, was discussed in Chapter 4. By now, it should be clear to the reader that the steps involved in the finite element analysis of a general class of problems (e.g., single secondorder, single fourthorder, and a pair of secondorder equations) are systematic and can be implemented on a digital computer. Indeed, the success of the finite element method is largely due to the ease with which the analysis of a class of problems, without regard to a specific problem, can be implemented on a digital computer. For different geometries, boundary and initial conditions, and problem data, a specific problem from the general class can be solved by simply supplying the required input data to the program. For example, if we develop a general computer program to solve equations' of the form (7.32) I
then all physical problems described by (3.1) and (4.1) and their timedependent versions can be solved for any compatible boundary and initial conditions. The purpose of this section is to discuss the basic steps involved in the development of a computer program for second and fourthorder onedimensional differential equations studied in the preceding chapters. The ideas presented here are used in the development of the model program FEM1DV2 (a revised version of FEM1D from the first edition of this book), and they are meant to be illustrative of the steps used in a typical finite element program development. One can make use of the ideas presented here to develop a program of one's own.
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
259
PREPROCESSOR • Read the input data • Echo the input data • Write the input data and plot the mesh
PROCESSOR • Compute element coefficient matrices and column vectors • Assemble element equations • Impose boundary conditions • Solve equations or find eigenvalues
POSTPROCESSOR • Compute solution at points other than the nodes • Compute the gradient of the solution • Print/plot the results
7.2.2
FIGURE 7.6 The three main functional units of a finite element program.
General Outline
A finite element program consists of three basic parts (see Fig. 7.6):
1. Preprocessor 2. Processor 3. Postprocessor In the preprocessor part of the program, the input data of the problem are read in and/or generated. This includes the geometry (e.g., length of the domain and boundary conditions), the data of the problem (e.g., coefficients in the differential equation), finite element mesh information (e.g., element type, number of elements, element length, coordinates of the nodes, and connectivity matrix), and indicators for various options (e.g., print, no print, type of field problem analyzed, static analysis, eigenvalue analysis, transient analysis, and degree of interpolation).
260
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
In the processor part, all steps in the finite element method, except for postproeessing, discussed in the preceding chapters are performed. These include the following: 1. 2. 3. 4.
Generation of the element matrices using numerical integration. Assembly of element equations. Imposition of the boundary conditions. Solution of the algebraic equations for the nodal values of the primary variables.
In the postprocessor part of the program, the solution is computed by interpolation at points other than nodes, secondary variables that are deriv~ble from the solution are computed, and the output data are processed in a desired format for printout and/or plotting. The preprocessor and postprocessors can be a few Fortran statements to read and print pertinent information, simple subroutines (e.g., subroutines to generate mesh and compute the gradient of the solution), or complex programs linked to other units via disk and tape files. The processor, where typically large amounts of computing time are spent, can consist of several subroutines, each having a special purpose (e.g., a subroutine for the calculation of element matrices, a subroutine for the imposition of boundary conditions, and a subroutine for the solution of the equations). The degree of sophistication and the complexity of a finite element program depend on the general class of problems being programmed, the generality of the data in the equation, and the intended user of the program. It is always desirable to describe, through comment statements, all variables used in the computer program. A flow chart of the computer program FEM1DV2 is presented in Fig. 7.7. In the following sections, a discussion of the basic components of a typical finite element program is presented, and then the basic ideas are illustrated via FORTRAN statements (see Appendix 1) ..
7.2.3 Preprocessor The preprocessor unit consists of reading input data and generating finite element mesh, and printing the data and mesh information. The input data to a finite element program consist of element type IELEM (i.e., Lagrange element or Hermite element), number of elements used (NEM), specified boundary conditions on primary and secondary variables (number of boundary conditions, global node number and degree of freedom, and specified values of the degrees of freedom), the global coordinates of nodes Xl> and element properties (e.g., coefficients a., be' Ce, fe, etc.). If a uniform mesh is used, the length of the domain should be read in, and global coordinates of the nodes can be generated in the program. The preprocessor portion that deals with the generation of finite element mesh (when not supplied by the user) can be separated into a subroutine
r.,
:
I
~
/
I I
PREPROCESSOR I I I
1
I I
I
I I
I I I
I I
r        , I I
I
I I I
I I I I
POSTPROCESSOR
ASSMBL
I I I I
1 I I
I
I I
I
j
I I
I
l
I I
I J
I
'j
I
AXLBX:
Subroutine {o~ determining eigenvalues and eigenvectors of the equation [AJ(X} = L[B]{X}.
BONDRY:
Subroutine for imposing specified boundary conditions of essential type, natural type, and mixed type.
COEFNT:
Subroutine for computing clement matrices and Vectors for all problems except truss and frame problems. Subroutine to echo the input file.
Subroutine to generate the mesh (i.e , global coordinates 0/ the global nodes and the connectivity array).
SOLVER or AXLBX
,2.. 51 C!
n
;>
SHPlD
I .r......... ::
PSTPRC:
Subroutine to postprocess the solution for all' problems, except truss and frame problems.
REACTN:
Subroutine to calculate the reactions for truss and frame problems.
1 SHP1D:
Subroutine to compute the shape functions and their derivatives.
SOLVER:
Subroutine to solve Iinear algebraic equations.
TRSFRM:
Subroutine to compute element stiffness matrices and force vectors for truss .and frame problems.
r
I Z
~
I:;)
I
I
ECHO: MESH1D:
~.
I
I I
1
++
I
I
!
COEFNT or TRSFRM BONDRY
I I
Subroutine for the assembly of element equations: upperbanded form for static and transient problems: full matrix form for eigenvalue problems.
I
PROCESSOR
1
ECHO and MESHlD
ASSMBL:
PSTPRC Or REACTN
I I
n
o
~
";::
:;; e
m
zm
~
FIGURE 7.7 Flow chart of the computer program FEMIDV2 along with the subroutines.
~
~
262
FINITE ELEMENT ANALYSIS 01' ONE·DIMENSIONAL PROBLEMS
(MESH1D), depending on the convenience and complexity of the program. Mesh generation includes computation of the global coordinates X, and the connectivity array [B} = [NOD]. Recall that the connectivity matrix describes the relationship between element nodes to global nodes: NOD(I, J)
= global node number corresponding to the Jth node
of
element I This array is used in the assembly procedure as well as to transfer information from element to the global system and vice versa. For example, to extract the element nodal coordinates x7 of the element Q" from the global coordinates Xl> we can use the array NOD:
x7 = Xl>
1= NOD(n, i)
The arrays {ELX} and {GLX} are used in FEM1DV2 for {x7} and {Xl}, respectivel y.
7.2.4 Calculation of Element Matrices (Processor) The most significant part of a processor is where we generate' element matrices. The element matrices are computed in various subroutines (COEFNT and TRSFRM), depending on the type of problem being solved. These subroutines typically involve numerical evaluations of the element matrices [Ke] and [Me] (program variables ELK and ELM) and the element vector {Ie} (program variable ELF) for various field problems. The Gauss quadrature described in Section 7.1.5 is used to evaluate element matrices and vectors, and the arrays are assembled as soon as they are computed. Thus, a loop on the number of elements in the mesh (NEM) is used to compute element matrices and assemble them (subroutine ASSMBL). It is here that the connectivity array NOD plays a crucial role. By putting element matrices into global locations one at a time, we avoid the computation of all element matrices at once. Element matrices for different model equations (MODEL) and type of problem (NTYPE) are generated by assigning values as follows:
1. MODEL = 1, NTYPE = 0: all field problems described by the model equation (3.1), including radially symmetric heattransfertype problems. 2. MODEL = 1, NTYPE = 1: radially symmetric elastic disk problems (see Problem 3.33). 3. MODEL = 2, NTYPE = 0 (RIE) or 2 (CIE): the Timoshenko theory of beams. 4. MODEL = 2, NTYPE = 1 (RIE) or 3 (CIE): the Timoshenko theory for bending of circular plates. 5. MODEL = 3, NTYPE = 0: the EulerBernoulli theory of beams. 6. MODEL = 3, NlYPE = 1: the EulerBernoulli theory for bending of circular plates.
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
263
7. MODEL = 4, NTYPE = 0: the twonode truss element. 8. MODEL = 4, NTYPE = 1: the twonode EulerBernoulli frame element. 9. MODEL = 4, NTYPE = 2: the twonode, Timoshenko frame element.
The timedependent option is exercised through variable ITEM: ITEM = 0 ITEM = 1 ITEM = 2
static analysis firstorder time derivative (i.e., parabolic) problems secondorder time derivative (i.e., hyperbolic) problems
The element matrices are evaluated using the Gauss quadrature, except for MODEL = 4, where the explicit forms of element coefficients are programmed in the interest of computational efficiency. The element shape functions SF and their derivatives GDSF are evaluated at the Gauss points in subroutine SHP1D. The Gaussian weights and points associated with two, three, four, and fivepoint integration are stored in arrays GAUSWT and GAUSPT, respectively. The nth column of GAUSWT, for example, contains the weights corresponding to the npoint Gauss quadrature rule: GAUSPT(i, n)
= ith Gauss point corresponding to the zrpoint Gauss rule
The variable NGP is used to denote the number of Gauss points, which is selected to achieve good accuracy. As noted earlier, the linear, quadratic, and cubic interpolation functions require two, three, and four quadrature points, respectively, to evaluate the element coefficients exactly. Thus, if IELEM is the element type, I for linear IELEM = 2 for quadratic { 3 for cubic
(Lagrange elements)
n
[see (7.29)] exactly then NGP = IELEM + 1 would evaluate KJj , MJj , and when c(x) is linear, and a(x), b (x), and f(x) are quadratic functions. The Hermite cubic element is identified with IELEM = 0, for which case NGP is taken to be 4. The coefficients a(x) = AX, b(x) = BX, and c(x) = CX, together with f(x) = FX in the differential equation (7.32) are assumed to vary with x as follows: BX=BXO+BXl*X
(a=aO+alx) (b =bo+b1x)
CX=CXO+CXl*X
(C=CO+C1X)
FX= FXO + FXl *X + FX2*X*X
(f =/0+ hx + Jzx 2 )
AX=AXO+AXl*X
For radially symmetric elasticity problems, (AXO, AXl) [or (BXO, EXl) for circular plates] are used to input Young's modulus E and Poisson's ratio Y.
264
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
The Gauss quadrature formula (7.20) can be implemented on the computer as follows. Consider Kij of the form (7.33)
Let us use the following program variables for the quantities in (7.33): SF(I)
ex = c(x),
AX=a(x), NPE = n, After transforming x to
= 1jJ1,
d·l,e GDSF(I) =~'I'~/ dx ELX(I) =xi
the number of nodes in the element
s, n
X
=
L x~ljJ't
[= ~XA +1h e (1 + s)]
(7.34)
/=1
the coefficients Kij in (7.33) take the form [see (7,25)]
(7.35a) NGP
=
L
Fij(s/)JWi
(7.35b)
/=1
where F't; denotes the expression in the square brackets in (7.35a), .:J is the Jacobian, and (SI, WI) are the Gauss points and weights. Examination of (7 .35b) shows that there are three free indices: i, i, and 1. We take the Gausspoint loop on I as the outermost one.' Inside this loop, we evaluate F'ft at the Gauss point for each i and j, multiply with the Jacobian J = !he and the weights Wi, and sum:
s/
ELK(i, j) = ELK(i, j) + Fij(sl)JWi
(7.36)
To accomplish this, we must initialize all arrays that are being evaluated using the Gauss quadrature: DO 10 J DO 10 I 10 ELK(IrJ)
= =
=
1,NPE l,NPE 0.0
This initialization is made outside the doloop on number of Gauss points. The computation of coefficients F'ft in (7 .35b) requires evaluation of a, c, 1jJ/, and d1jJJds at the Gauss point SI' Hence, inside the loop on I, we call subroutine SHP1D to evaluate 1jJ;, d1jJJdx = (d1JJJd;)/,1. We now have all
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
265
quantities needed to compute Kij in (7.35b); DO 100 NI = l,NGP XI = GAUSPT(NI,NGP)
c c C
call subroutine SHPID to evaluate the interpolation functions (SF) and their global derivatives (GDSF) at the Gauss point XI
C
CALL SHPID(XI,NPE,SF,GDSF,GJ) CONST = GJ*GAUSWT(NI,NGP) C
x = 0.0 DO 20 I=l,NPE 20 X = X + SF(I)*ELX(I)
C
AX=AXO + AXl*X CX=CXO + CXl*X DO 30 J = I,NPE ELF(J) = ELF(J) + CONST*SF(J)*FX DO 30 I = 1,NPE 30 ELK(I,J)=ELK(I,J)+CONST* (AX*GDSF(I) *GDSF(J) +CX*SF(I) *SF (J))
In the same way, all other coefficients (e.g., Mij andfD can be evaluated. Recall that the element properties (i.e., Kij, Mij, and ff) are calculated by calling a suitable subroutine (COEFNT or TRSFRM) for the field problem being analyzed within a loop on number of elements. As soon as the element properties are available for a particular element, they are put into their proper locations (i.e., assembled) with the help of array NOD. The assembly is explained in the next section.
7.2.5
Assembly of Element Equations (Processor)
The assembly of element equations should be carried out as soon as they are computed, rather than waiting till element coefficients of all elements are computed. The latter requires storage of the element coefficients of each element. In the former case, we can perform the assembly in the same loop in which a subroutine is called to calculate element matrices. A feature of the finite element equations that enables us to save storage and computing time is the assembly of element matrices in upperbanded form. When element matrices are symmetric, as is the case in most problems of interest in this book, the resulting global (or assembled) matrix is also symmetric, with many zeros away from the main diagonal. Therefore, it is sufficient to store only the upper halfband of the assembled matrix. The halfbandwidth of a matrix is defined as follows. Let N; be the number of matrix elements between the diagonal element and the last nonzero element in the ith row, after which all elements in that row are zero; the halfbandwidth is the maximum of (N; + 1) X NDF, where NDF is the number of degrees of freedom per node:
b, = max [(N; + 1) X NDF] l~i";;;;n
266
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
where 11 is the number of rows in the matrix (or equations in the problem). Generalpurpose equation solvers are available for such banded systems of equations. The halfbandwidth NHBW of the assembled (i.e., global) finite element matrix can be determined in the finite element program itself. The local nature of the finite element interpolation functions (i.e., 1fJ'f are defined to be nonzero only over the element Qe) is responsible for the banded character of the assembled matrix. Recall from earlier discussions that if two global nodes do not belong to the same element then the corresponding entries in the global matrix are zeros: KlJ = 0 if global nodes 1 and J do not correspond to local nodes
of the same element "
This property enables us to determine the halfbandwidth 'NHBW of the assembled matrix: (
NHBW =
max
{abs [NOD(N, 1)  NOD(N, J)]
+ 1} x NDF
(7.37a)
I ",N"'NEM I~I,J~NPE
where NEM = number of elements in the mesh NPE = number of nodes per element
(7.37b)
NDF = number of degrees of freedom per element Clearly, for onedimensional problems with elements connected in series, the maximum difference between nodes of an element is equal to NPE  1. Hence, NHBW= [(NPE1) + 1] x NDF=NPExNDF
(7.38)
Of course, NHBW is always less than or equal to the total number of primary degrees of freedom, i.e., the number of equations, NEQ, in the finite element mesh of the problem. The logic for assembling the element matrices Kij into the upperbanded form of the global coefficients K lJ is that the assembly can be skipped whenever J < I and J > NHBW. The main diagonal, 1 = J, of the assembled square matrix (i.e., full storage form), becomes the first column of the assembled banded matrix (i.e., banded storage form), as shown in Fig. 7.8. The upper diagonals (parallel to the main diagonal) take the position of respective columns in the banded matrix. Thus, the banded matrix has dimension NEQ X NHBW, where NEQ denotes the total number of equations in the problem. . The element coefficients Kif and f? of a typical element Q" are to be assembled into the global coefficient matrix [K) and source vector {F}, respectively. If the ith node of the element is equal to the Ith global node, and the jth node of the element is equal to the Jth global node, we have KlJ = Kij,
Pi = Fi_
(for NDF = 1)
(7.39a)
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
Upperhalfbanded mode
Full storage mode ~ Halfbandwidth = NHBW X X
,X
"
,
x x x x
All blanks contain zeros
X X
'x X x ""'.
X X X
X X X X X X X
.
"'" All entries above x, x x '~this d.iagonal x x x x contain zeros
'x
, x xx'>.x, x "x, x x"' "
X X X
x x x x X X X
X
X X
X
X
X
,
X
X
X
x x x x
X
"x x x >. " "
Symmetric
267
X
X
X
X
X
X X X X
X
x,x X X, X X
X
X
X
X X X'~ X X·/ I // ( X .< J
~,
l:
.>, Main diagonal NEQ x NEQ
NEQ x NHBW
FIGURE 7.8 Finite element coefficient matrix storage in upperhalfbanded form.
The values of I and J can be obtained with the help of array NOD: 1= NOD(n, i),
J = NOD(n, j)
(7.39b)
Recall that it is possible that the same I and J may correspond to a pair of i and j of some other element Q"'. In that case, K~'j will be added to existing coefficients KIf during the assembly. For NDF> 1, the logic still holds, with the change K(NR)(NC)
= K(i+pl)(j+q~l)
= 1, 2, ' .. s NDF)
(7.40a)
NC = (J  1) X NDF + q
(7.40b)
(p, q
where NR
= (I 
1)
X
NDF + p,
and I and J are related to i and j by (7 .39b). A Fortran listing of the subroutine ASSMBL can be found in FEM1DV2 in Appendix 1.
7.2.6 Imposition of Boundary Conditions (Processor) Imposition of boundary conditions on the primary and secondary global degrees of freedom can be carried out through a subroutine (BONDRY), which remains unchanged for 2D or 3D problems. There are three types of boundary conditions for any problem:
1. Essential boundary conditions, i.e., boundary conditions on primary variables (Dirichlet boundary conditions).
268
FINITE ELEMENT ANALYSIS OF ONE·DIMENSIONAL PROBLEMS
2. Natural boundary conditions, i.e., boundary conditions on secondary variables (Neumann boundary conditions). 3. Mixed boundary conditions (Newton boundary conditions). The procedure for implementing the boundary conditions on the primary variables involves modifying the assembled coefficient matrix (GLK) and righthand column vector (GLF) by three operations: Step 1 Step 2
Step 3
moving the known products to the righthand column of the matrix equation. replacing the columns and rows of GLK corresponding to the known primary variable by zeros, and setting the coefficient on the main diagonal to unity. "replacing the corresponding component of the righthand column by the specified value of the variable.
Consider the following N algebraic equations in full matrix form: 13
K2ll1 K12 K K 22 [ K·3 1 K 32
K •• ']{Vl} K2J ... U2 K 33 ••• V3
·· . , .. . ~.
...
_
{Pi} IS. ~
.. .
where VJ and Fj are the global primary and secondary degrees of freedom, and KlJ are the assembled coefficients. Suppose that Us = is specified. Recall that when the primary degree of freedom at a node is known, the corresponding secondary degree of freedom is unknown, and vice versa. Set K ss = 1 and Fs = Os; further, set K SJ = K JS = 0 for 1= 1, 2, ... , N and 1=1= S. For S = 2, the modified equations are
o,
Ku
o
o
where
F;=F;Ki2 02 (i=1,3,4,5, ... ,n;i=f=2) Thus, in general, if Us =
Kss = 1,
o, is known, we have
Fs = Us; F; = F;  KiSOS; KSi = KiS = 0 where i = 1, 2, ... , S '1, S + 1, ... , n (i =1= S). This procedure is repeated for .every specified primary degree of freedom. It enables us to retain the original order of the matrix, and the specified boundary conditions on the primary degrees of freedom are printed as part of the solution. Of course, the logic should be implemented for a banded system of equations.
NUMERICAL INTEGRATION AND COMPlJfER IMPLEMENTATION
269
The specified secondary degrees of freedom (Qj) are implemented directly by adding their specified values to the computed values. Suppose that the point source corresponding to the Rth secondary degree of freedom is specified to be FR' Then
FR = fR+FR where fR is the contribution due to the distributed source f(x); fR is computed as a part of the element computations and assembled. Mixedtype boundary conditions are of the form du a dx + k(u  u) = 0 (u and k are known constants)
(7.41)
which contains both the primary variable u and the secondary variable a du/dx. Thus a dul dx at the node P is replaced by kp(up  up): Qp = kp( Up  Up)
This amounts to modifying Kpp by adding kp to its existing value, K pp oE K pp + k p
F;, oE F'p + k p Up All three types of boundary conditions are implemented in subroutine BONDRY for boundary, initial, and eigenvalue problems. The following are used in subroutine BONDRY (see Appendix 1): NSPV NSSV NNBC
VSPV VSSV VNBC
UREF ISPV
number of specified primary variables number of specified secondary variables number of Newton boundary conditions column of the specified values 0. of primary variables column of the specified values FR of secondary variables column of the specified values k p column of the specified values Op array of the global node and degree of freedom at the node that is specified [ISPV(I, 1) = global node of the lth boundary condition, ISPV(I, 2) = degree of freedom specified at the global node, ISPV(I, 1)]
(7.42)
Similar definitions are used for ISSV and INBC arrays.
7.2.7 Solution of Equations and Postprocessing Subroutine SOLVE is used to solve a banded system of equations, and the solution is returned in array GLF. The program performs the Gaussian elimination and backsubstitution to compute the solution. For a discussion of the Gaussian elimination used to solve a set of linear algebraic equations, the
270
FINITE ELEMENT ANALYSIS OF ONE·DlMENSIONAL PROBLEMS
reader is referred to the book by Carnhan, Luther, and Wilkes (1969). On most computing systems, a variety of equation solvers are available, and one can use any of the programs that suits the needs. See Appendix 1 for a listing of subroutine SOLVER. Postprocessing involves computation of the solution and its gradient at preselected points of the domain. Subroutine PSTPRC is used to evaluate the solution and its derivatives in any element: n
ue(x) = 2: uJ1JJj(x),
I t uj (d1JJj) I dx
( due) = dx
j=1
x
1=1
(7.43a)
x
for the Lagrange elements and
e
dmwm = ~ " uj (d"'tP~) ~d d m l \ (m = 1, 2,3) (7.....43b) 1=1 x 1=1 X x / for the Hermite cubic elements. The nodal values uJ of the element Qe are ( deduced from the global nodal values VI as follows: n we(x) = 2: ujtPj(x),
uJ= VI>
1= NOD(e, j),
For NDF> 1, I is given by I uj+p
when NDF = 1
= [NOD(e, j) 1] X NDF and
= V/+p (p = 1, 2, ... , NDF)
The values computed using the derivatives of the solution are often inaccurate because the derivatives of the approximate solution become increasingly inaccurate with increasing order of differentiation. For example, the shear force computed in the EulerBernoulli beam theory,
d (d
2
W)
V = dx b dx2 =
2tPe) d ( d b dxi
j~ uJ dx n
(7.44)
will be in considerable error compared with the true value of V. The accuracy increases, rather slowly, with mesh refinement and higherorder elements. When accurate values of the secondary variables are desired, it is recommended that they be computed from the element equations: n
Qi =
L KijuJ  fi
(i = 1,2, ... , n)
(7.45)
j=1
However, this requires recomputation or saving of element coefficients Kij and
/1. 7.3 APPLICATIONS OF THE COMPUTER PROGRAM FEMIDV2
7.3,.1 General Comments The computer program FEM1DV2, which embodies the ideas presented in the previous section, is intended to illustrate the use of the finite element models developed in Chapters 36 to a variety of onedimensional field problems, some of them not discussed in this book. FEM1DV2 is a modified version of
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
271
FEMID from the first edition of this book. It was developed as a learning computational tool for students of the course. In the interest of simplicity and ease of understanding, only the model equations discussed in this book and their immediate extensions are included in the program (see Appendix 1). Table 7.3 contains a summary of the definitions of coefficients of various model problems and their corresponding program variables. The table can be used as a ready reference to select proper values of AXO, AXl, and so on for different problems.
7.3.2 Illustrative Examples Here we revisit some of the example problems considered earlier to illustrate the use of FEMIDV2 in their solution. Only certain key observations concerning the input data are made, but complete listings of input files for each problem are given. In the interest of brevity, the output files for most problems are not included. Table 7.4 contains a description of the input variables to program FEMlDV2. Example 7.1 Steady heat transfer in a fin. The problem is governed by
dze z+ce.=o for O<x
z
'='
8 ""
,~. ~ ttl
"" iii
~
N
t;J
274
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
TABLE 7.4
Description of the input variables to the program FEMIDV2t
*
Data Card 1: TITLE  Title of the problem being solved (80 characters)
*
Data Card 2: MODEL  Model equation being solved (see below) NTYPE  Type of problem solved (see below) MODEL=l, MODEL=l, MODEL=l, MODEL=2, MODEL=2, MOOEL=2, MODEL=2, MODEL=3, MODEL=3, MODEL=4, MODEL=4, MODEL=4,
NTVPE=O: NTYPE=l: NTYPE>l: NTYPE=O: NTYPE=l: IlTYPE=2: IlTYPE>2: NTYPE=O: NTYPE>O: NTYPE=O: IlTYPE=l: NTYPE=2:
A A A A A A A A A A A A
problem of MODEL EQUATION 1 circular DISK (PLANE STRESS) circular DISK (PLANE STRAIN) Timoshenko BEAM (RIEl) problem Timoshenko PLATE (RIE) problem Timoshenko BEAM (CIE#/) problem Timoshenko PLATE (CIE) problem EulerBernoulli BEAM problem EulerBernoulli Circular plate plane TRUSS problem EulerBernoulli 'FRAME problem Timoshenko (CIE) FRAME problem
I  Reduced Integration Element (RIE) II  Consistent Interpolation Element (CIE)
ITEM  Indicator for transient analysis: ITEM=O, ITEM=l, ITEM=2, ITEM=3,
*
Steadystate solution Transient analysis of PARABOLIC equations Transient analysis of HYPERBOLIC equations Eigenvalue analysis
Data Card 3: IELEM  Type of finite element: IELEM=O, Hermite cubic finite element IELEM=l, Linear Lagrange finite element IELEH=2, Quadratic Lagrange finite element IlEM
*
 Number of elements in the mesh
Data Card 4: rCOIlT  Indicator for continuity of data for the problem: ICONT=l, Data (AX,BX,CX,FX and mesh) is continuous ICONT=Q, Data is element dependent NPRNT  Indicator for printing of element/global matrices: NPRNT=O, Not print element or global matrices but postprocess the solution and print NPRNT=l, Print Element 1 coefficient matrices only but postprocess the solution and print
t SKIP means that you omit the input data. It does not mean that you leave a blank. In the "free format" used here, variables of each data card are read from the same line; if the values are not found on the same line, the computer will look for them on the next line(s). Each data card should start with a new line. Note that the variable names included in the data sets of Tables 7.5~7.14 are only for the convenience of the readersthey are nOI read by the computer.
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
275
TABLE 7.4 (Continued)
Print Element 1 and global matrices but NOT postprocess the solution NPRNT>2, Not print element or global matrices and 'NOT postprocess the solution
NPRNT~2,
SKIP Cards 5 through 15 for TRUSS/FRAME problems (MODEL ~ 4), and read Cards 5 through 15 only if MODEL.NE.4. SKIP Cards 5 through 9 if data is discontinuous (ICONT ~ 0)
*
Data Card 5: DX(I)  Array of element lengths. DX(l) denotes the global coordinate of Node 1 of the mesh; DX(I) (I~2,NEMl) denotes the length of the (Il)st element. Here NEMl~NEM+l, and NEM~number of elements in the mesh. Cards 6 through 9 define the coefficients in the model equations. All coefficients are expressed in terms of GLOBAL coordinate x. See Table 7.3 for the meaning of the coefficients, especially for deformation of circular plates and Timoshenko elements. __
*
Data Card 6: AXO AXl
*
Data Card 7: BXO BXl
*
 Constant term of the coefficient AX  Linear term of the coefficient AX

Constant term of the coefficient BX  Linear term of the coefficient BX
Data Card 8: CXO CXl
 Constant term of the coefficient ex  Linear term of the coefficient CX
SKIP Card 9 for eigenvalue problems (Le. ITEM=3)
*
Data Card 9: FXO FXl FX2
 Constant term of the coefficient FX  Linear term of the coefficient FX  Quadratic term of the coefficient FX
SKIP Cards 10 through 15 if data is continuous (ICONT.ME.O). Cards 10 through 15 are read for each element (i.e., NEM times). All coefficients are expressed in terms of the GLOBAL coord _
*
Data Card 10: {GLX}  Global xcoordinates of the FIRST and LAST nodes of the element
*
Data Card 11: {NOD}  Connectivity of the element: NOD(N,I}=Global node number corresponding to the Ith node of 'Element N (I=l,NPE) NPE.denotes the Number of nodes Per Element
*
Data Card 12: {DCAX} Constant and linear terms of the coefficient AX
*
Data card 13: {DCBX} Constant and linear terms of the coefficient BX
276
FINITE ELEMENT ANALYSIS OF ONEDIMENSIONAL PROBLEMS
TABLE 7.4 (Continued)
*
Data Card 14: {DCCX} Constant and linear terms of the coefficient CX
*
Data Card 15: {DCFX} constant, linear and quadratic terms of FX READ Cards 16 through 21 for TRUSS/FRAME problems (MODEL SKIP Cards 16 through 21 if MODEL.NE.4
4).
* Data Card 16: NNM
 Number of nodes in the finite element mesh
SKIP Cards 17 through 19 for TRUSS problems (NTYPE = 0). Cards 19 through 19 are read for each element, i.e., NEM times
* Data Card 17: PR SE SL SA 81 CS SN

/
Poisson's ratio of the material# Young's modulus of the material Length of the element Crossectional area of the element Moment of inertia of the element cosine of the angle of orientation of the element Sine of the angle of orientation of the element Angle is measured counterclackwise from x axis
# PR is not used for the Euler Bernoulli element
*
Data Card 18: HF VF PF XB
*

Intensity of the horizontal distributed force Intensity of the transversely distributed force Point load on the element Oistance from node I, along the length of the element, to the point of load application, PF CMT  cosine of the angle of orientation of the load PF SNT  sine of the angle of orientation of the load PF Angle is measured counterclackwise from the element x axis Data Card 19: {NOD}  connectivity of the element: NOD(N,I) Global node number corresponding to the Ith node of Element N (I=1,NPE) READ Cards 20 and 21 only for TRUSS problems (NTYPE = 0). Cards 20 and 21 are read for each element;i.e. NEM times
*
Data Card 20: SE SL SA CS SN HF
*

Young's modulus of the material Length of the element Crossectional area of the element Cosine of the angle of orientation of the element Sine of the angle of orientation of the element Angle is measured counterclackwise from x axis  Intensity of the horizontal distributed force
Data Card 21: {NOD}  connectivity of the element: NOD(N,I) = Global node number corresponding to the 1th node of Element N (I=l,NPE)
NUMERICAL INTEGRATION AND COMPUTER IMPLEMENTATION
277
TABLE 7.4 (Continued)
* Data Card 22: NSPV
 Number of specified PRIMARY degrees of freedom
SKIP Card 23 if no primary variables is specified (NSPV=O). Repeat Card 23 NSPV times _ _~ ~ _ * Data Card 23:
(I = 1 to NSPV)
ISPV(I,l)  Node number at which the PV is specified ISPV(I,2)  specified local primary DOF at the node VSPV(I)  specified value of the primary variable (PV) (will not read for eigenvalue problems) SKIP. Card 24 for eigenvalue problems (i.e. when ITEM=3)
*
Data Card 24: NSSV
 Number of specified (nonzero) SECONDARY variables
SKIP Card 25 if no secondary variable is specified (NSSV=O). Repeat Card 25 NSSV times _ * Data Card 25:
(I = 1 to NSSV)
ISSV(I,l)  Node number at which the SV is specified ISSV(I,2)  specified local secondary DOF at the node VSSV(I)  specified value of the secondary variable (PV)
*
Data Card 26: NNBC
 Number of the Newton (mixed) boundary conditions
SKIP Card 27 if no mixed boundary condition is specified (NNBC=O). The mixed boundary condition is assumed to be of the form: SV+VNBC *(PVUREF) = O. Repeat Card 28 NNBC times
*
Data Card 27: INBC(I,l) INBC(I,2) VNBC(I) UREF(I)
(I = 1 to NNBC) 
Node number at which the mixed B.C. is specified Local DOF of the PV and SV at the node Value of the coefficient of the PV in the B.C. Reference value of the PV
SKIP Card 28 if ITEM=O (read value problems
*
only for timedependent and eigen
_
Data Card 28: CTO  Constant part of CT = CTO + CTl*X CTI  Linear part of CT = CTO + CTl*X SKIP remaining cards if steadystate or eigenvalue analysis is to be performed (ITEM=O or ITEM=3)
* Data Card 29: DT  Time increment (uniform) ALFA  Parameter in the timeapproximation scheme BETA  Par~meter in the time approximation scheme
*
Data Card 30: INCOND Indicator INCOND=O, 1NCOND>O, NTIME  Number of INTVL  Time step
for initial conditions: Homogeneous (zero) initial conditions Nonhomogeneous initial conditions time steps for Which solution is sought intervals at which solution is to printed
278
ANITE ELEMENT ANALYSIS OF ONE·DIMENSIONAL PROBLEMS
TABLE 7.4 (Continued) SKIP Cards 31 and 32 if initial conditions are zero (INCOND=O) ___
*
Data Card 31: {GUO}  Array of initial values of the primary variables SKIP Card 32 for parabolic equations (ITEM=l)
*
Data Card
32~
{GU1}  Array of initial values of the first timederivatives of the primary variables
For Set 2 boundary conditions, we have one essential boundary condition and one mixed boundary condition: ;' "VI = 0,
Qi+
t
(Us 
V~) = 0
Hence, NSPV = 1, NSSV = 0, and NNBC = 1, and the values in the mixed boundary conditions are input as
VNBC(l) = fJ/k
= 1.28,
UREF(1) = V~= 0.0,
INBC(l, 1) = 5,
INBC(I, 2)
= 1.
The output file for Set 1 boundary conditions is presented in Table 7.5(a). The input files are echoed in the output. The input file for Set 2 boundary conditions is given in Table 7.5(b). Example 7.2 (Example 3.2) Radially symmetric heat transfer in a solid cylinder. The governing equation of the problem is given by [see (3.78)] 
~ (21I'kr ~~ =
2:rqo for O', 6t = 10 6 s, and assume that the beam is of square crosssection 0.5 in by 0.5 in. Using five EulerBernoulli elements in the halfbeam, obtain the finite element solution and compare with the series solution at midspan for the case c = !L. 7.53. Repeat Problem 7.52 for c = ~L and eight elements in the full span. 7.54. Repeat Problem 7.52 for P(t) = Po at midspan and eight elements in the full span. 7.55. Consider a cantilevered beam with a point load at the free end. Using the load and data of Problem 7.52, find the finite element solution for the transverse deflection using EulerBernoulli beam elements. 7.56. Repeat Problem 7.52 for a clamped beam with the load at the midspan. 7.577.61. Solve Problems 7.527.56 using Timoshenko beam elements. Use v = 0.3.
REFERENCES FOR ADDITIONAL READING Fluid mechanics Bird, R. B., W. E. Stewart, and E. N. Lightfoot; Transport Phenomena, John Wiley, New York, 1960. Duncan, W. J., A. S. Thorn, and A. D. Young: Mechanics of Fluids, 2d ed., Elsevier, New York, 1970. Eskinazi, S.: Principles of Fluid Mechanics, Allyn and Bacon, Boston, 1962. Harr, M. E.: Ground Water and Seepage, McGrawHili, NeVI York, 1962. Nadai, A.: Theory of Flow and Fracture of Solids, vol. II, McGrawHill, New York, 1963. Schlichting, H.: BoundaryLayer Theory (translated by J. Kestin), 7th ed., McGrawHili, New York, 1979. Shames, I. H.: Mechanics of Fluids, McGrawHili, New York, 1962. Vallentine, H. R.: Applied Hydrodynamics, Butterworths, London, 1959. Verrujit, A.: Theory of Groundwater FlolV, Gordon and Breach, New York, 1970.
Heat transfer Carslaw, H. S., and J. C. Jaeger: Conduction of Heat ill Solids, Clarendon Press, Oxford, 1959. Holeman, J. P.: Heal Transfer, 6th ed., McGrawHili, New York, 1986. Kreith, F.: Principles of Heat Transfer, 3rd ed., Harper & Row, New York, 1973. Myers, G. G.: Analytical Methods in Conduction Heat Transfer, McGrawHill, New York, 1972. OZisik, M. N.: Heat Transfer: A Basic Approach, McGrawHili, New York, 1985.
292
FINITE ELEMENT ANALYSIS OF ONE· DIMENSIONAL PROBLEMS
Plane elasticity Budynas, R. G.: Advanced Strength and Applied Stress Analysis, McGrawHili, New York, 1977. Harris, C. M., and C. E. Crede: Shock and Vibration Handbook, vol. I, McGrawHili, New York, 1961. Ugural, A. C., and S. K. Fenster: Advanced Strength and Applied Elasticity, Elsevier, New York, 1975. Volterra, E., and J. H. Gaines: Advanced Strength of Materials, PrenticeHall, Englewood Cliffs, NJ, 1971.
Time approximations Bathe, K. J.: Finite Element Procedures in Engineering Analysis, PrenticeHall, Englewood Cliffs,
NJ,1982. Clough, R. W., and J. Penzien: Dynamics of Structures, McGrawHill, New York, 1975. Hughes, T. J. R.: The Finite Element Method, PrenticeHall, EnglewoodCliffs, NJ, 1987. Numerical integration
/
Carnahan, B., H. A. Luther, and J. O. Wilker: Applied Numerical Methods, John WHey, New York, 1969. ( Loxan, A. N., N. Davids, and A. Levenson: "Table of the Zeros of the Legendre Polynomials of Order 116 and the Weight Coefficients for Gauss' Quadrature Formula," Bulletin of the American Mathematical Society, vel, 48, pp. 739743, 1942. Silvester, P.: "NewtonCotes Quadrature Formulae for NDimensional Simplexes," Proceedings of 2d Canadian Congress of Applied Mechanics, Waterloo, Ontario, 1969. University of . Waterloo, Waterloo, Ontario, 1969. Stroud, A. H., and D. Secrest: Gaussian Quadrature Formulas, PrenticeHall, Englewood Cliffs,
NJ,1966. Zienkiewicz, O. C., and R. L. Taylor: The Finite Element Method, vols. 1 and 2, McGrawHili, 1989 and 1991.
PART
3 FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
293
CHAPTER
8 SINGLEVARIABLE PROBLEMS
8.1 INTRODUCTION The finite element analysis of twodimensional problems involves the same basic steps as those described for onedimensional problems in Chapter 3. The analysis is somewhat complicated by the fact that twodimensional problems are described by partial differential equations over geometrically complex regions. The boundary I' of a twodimensional domain Q is, in general, a curve. Therefore, the finite elements are simple twodimensional geometric shapes that can be used to approximate a given twodimensional domain as well as the solution over it. Thus, in twodimensional problems, we not only seek an approximate solution to a given problem on a domain but we also approximate the domain by a suitable finite element mesh. Consequently, in the finite element analysis of twodimensional problems, we shall have approximation errors due to the approximation of the solution as well as discretization errors due to the approximation of the domain. The finite element mesh (discretization) consists of simple twodimensional elements, such as triangles, rectangles, and/or quadrilaterals, that allow unique derivation of the interpolation functions. The elements are connected to each other at nodal points on the boundaries of the elements. The ability to represent domains with irregular geometries by a collection of finite elements makes the 295
TABLE 8.1
Some examples of the Poisson equation V· (kVu) = f Natural boundary condition: k au/ an + f3(u  uoo) = q. Essential boundary condition: u = it Primary variable Field or application
u
Material constant k
1. Heat transfer
Temperature T
Conductivity k
2. Irrotational flow of an ideal fluid
Stream function 'l/J
Density p
f Heat source Q
Secondary variables
au au q, ax' ay Heat flow q Icomes from conduction k aT/an and convection h(T  T~)l Velocities:
Mass production (normally zero)
a'l/J ax
(J
a1J! ay
 =  v , :;;::u
Velocity potential 3. Groundwater flow
Source variable
tP
Piezometric head 4>
Density p
Mass production a (normall zero)
~=U
Permeability K
Recharge Q (or pumping, Q)
Seepage q = k
a4> ax
a4> =v 'ay
~:
Velocities:
a4> ax
a4J ay
u=k~, v=k~
4. Torsion of members with constant crosssection
Stress·function 'II
4J
5. Electrostatics
Scalar potential
6. Magnetostatics
Magnetic potential
7. Transverse deflection of elastic.membranes
Transverse deflection u
k=l
1=2
G = shear modulus
o ':' angle of twist per GO aw = a umt length ay xz
Dielectric constant e Charge density
4J
a'll GO ax = ayz
p
Permeability /.I.
Charge density p
Tension Tin membrane
Transversely distributed load
Displacement flux density D n Magnetic flux density s,
t!'
Normal force q
SINGLE·VARIABLE PROBLEMS
297
method a valuable practical tool for the solution of boundary, initial, and eigenvalue problems arising in various fields of engineering. The objective of this chapter is to extend the basic steps discussed earlier for onedimensional problems to twodimensional boundary value problems involving a single dependent variable. Once again, we describe the basic steps of the finite element analysis with a model secondorder partial differential equation, namely the Poisson equation, governing a single variable. This equation arises in a number of fields, including electrostatics, heat transfer, fluid mechanics, and solid mechanics (see Table 8.1).
8.2 BOUNDARY VALUE PROBLEMS 8.2.1 The Model Equation Consider the problem of finding the solution differential equation
II
of the secondorder partial
(8.1) for given data aij (i, j = 1, 2), aoo and f, and specified boundary conditions. The form of the boundary conditions will be apparent from the weak formulation. As a special case, one can obtain the Poisson equation from (8.1) by setting all = an = a and a12 = a21 = aoo = 0:
V'(aVu)=f in
Q
(8.2)
where V is the gradient operator. If i and J denote the unit vectors directed along the x and y axes, respectively, the gradient operator can be expressed as (see Section 2.2.2) A
a
A
a
V=i+j
ax
ay
and (8.2) in the cartesian coordinate system takes the form
au)
~(a ~(a ax ax oy
au) =f ay
(8.3)
In the following, we shall develop the finite element model of (8.1). The major steps are as follows:
1. 2. 3. 4. 5.
Discretization of the domain into a set of finite elements. Weak (or weightedintegral) formulation of the governing differential equation. Derivation of finite element interpolation functions. Development of the finite element model using the weak form. Assembly of finite elements to obtain the global system of algebraic equations.
298 6. 7. 8.
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
Imposition of boundary conditions. Solution of equations. Postcomputation of solution and quantities of interest.
Steps 6 and 7 remain unchanged from onedimensional finite element analysis because at the end of Step 5 we have a set of algebraic equations whose form is independent of the dimension of the domain or nature of the problem. In the following sections, we discuss each step in detail.
8.2.2 Finite Element Discretization In two dimensions, there is more than one simple geometric shape that can be used as a finite element (see Fig. 8.1). As we shall see shortly,«, the interpolation functions depend not only on the number of nodef in the element but also on the shape of the element. The shape of the element must be such that its geometry is uniquely defined by a set of points, which serve as the element nodes in the development of the interpolation functions. As will be discussed later in this section, a triangle is the simplest geometric shape, followed by a rectangle. The representation of a given region by a set of elements (i.e., discretization or mesh generation) is an important step in finite element analysis. The choice of element type, number of elements, and density of elements depends on the geometry of the domain, the problem to be analyzed,
r
y
Discretization error
(a)
y
~: D'ds
1 x
r (b)
FIGURE 8.1 Finite element discretization of an irregular domain: (a) discretization of a domain by triangular and quadrilateral elements; (b) a typical triangular element (boundary P", the unit normal fi on the boundary of the element).
SINGLEVARIABLE PROBLEMS
299
and the degree of accuracy desired. Of course, there are no specific formulae for obtaining this information. In general, the analyst is guided by his or her technical background, insight into the physics of the problem being modeled (e.g., a qualitative understanding, of the solution), and experience with finite element modeling. The general rules of mesh generation for finite element formulations include the following: 1. 2.
3.
4.
Select elements that characterize the governing equations of the problem. The number, shape, and type (i.e., linear or quadratic) of elements should be such that the geometry of the domain is represented as accurately as desired. The density of elements should be such that regions of large gradients of the solution are adequately modeled (i.e., more or higherorder elements should be used in regions of large gradients). Mesh refinements should vary gradually from highdensity regions to lowdensity regions. If transition elements are used, they should be used away from critical regions (i.e., regions of large gradients). Transition elements are those that connect lowerorder elements to higherorder elements (e.g., linear to quadratic).
Additional discussion of finite element meshes and types of elements (linear and higherorder Lagrange elements, transition elements, etc.) will be presented in Chapter 9.
8.2.3 Weak Form In the development of the weak form, we need only consider an arbitrarily typical element. We assume that Qe is such an element, whether triangular or quadrilateral, of the finite element mesh, and we develop the finite element model of (8.1) over ~Y. Various twodimensional elements will be discussed in the sequel. Following the threestep procedure presented in Chapters 2 and 3, we develop the weak form of (8.1) over the typical element Qe. The first step is to multiply (8.1) with a weight function iv, which is assumed to be oncedifferentiable with respect to x and y, and then integrate the resulting equation over the element domain Qe: 0=
L, w[
:x (Ft) 
~ (F2) + aoou  f] dx dy
(8Aa)
where (8Ab)
In the second step, we distribute the differentiation equal1y between u and w. To achieve this, we integrate the first two terms in (8Aa) by parts. First we
300
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
note the identities or
aFl aw
a
ax ax
ax
w~=F1~(wF1)
(8.5a) (8.5h)
or
Next, we recall the component form of the gradient (or divergence) theorem in (2. 17b):
(8.6a) I'
Ln' vy: (WF;) dx dy =,[Try wF;n
y
ds
(8.6b)
where n, and n y are the components (i.e., direction cosines) of the unit normal vector ii = n) +
fd = cos at + sin aj
(8.7)
on the boundary I", and ds is the arclength of an infinitesimal line element along the boundary (see Fig. 8.th). Using (8.5) and (8.6) in (8Aa), we obtain
(8.8) From an inspection of the boundary term in this equation, we note that the specification of u constitutes the essential boundary condition, and hence u is the primary variable. The specification of the coefficient of the weight function in the boundary expression
(8.9) constitutes the natural boundary condition; thus, q" is the secondary variable of the formulation. The function q" = q,,(s) denotes the projection of the vector a1.Vu along the unit normal Ii. By definition, q" is positive outward from the surface as we move counterclockwise along the boundary I", The secondary variable q" is of physical interest in most problems. For example, in the case of heat transfer through an anisotropic medium, a'l are the conductivities of the medium, and q" is the heat flux normal to the boundary of the element (into the element).
SINGLEVARIABLE PROBLEMS
301
The third and last step of the formulation is to use the definition (8.9) in (8.8) and write the weak form of (8.1) as
(8.10a)
or B(w, u) = I(w)
where the bilinear form B(·, .) and linear form
(8.l0b)
10 are
I(W)=l wfdxdy+i. wqc ds Q
0 2a
(8.93a)
where
c = u1 uo,
b
= 3u1 + 4U2 
u~,
a = 2(ui 2U2 + u~)
(8.93b)
Equation (8.93a) is to be applied on any three lines in the element until three different values h > So > 0 are found. The computational problem considered here [i.e., (8.65)] has several physical interpretations (see Table 8.1). The problem can be viewed as one of finding the temperature u in a unit square with uniform internal heat generation, where the sides x = 0 and y = 0 are insulated and the other two sides are at zero temperature (see Section 8.4.1). Another interpretation of the equation is that it defines the torsion of a cylindrical bar of 2 in square crosssection (see Section 8.4.3). In this case, U denotes the stress function 'II,
FIGURE 8.15 Isolines for triangular and quadrilateral elements (linear elements).
334
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
and the components of the gradient of the solution are the stresses (which are of primary interest):
a
YZ
aw ax
=G(J
where G is the shear modulus and () is the angle of twist per unit length of the bar. A third interpretation of (8.65) is provided by groundwater (seepage) and potential flow problems. In this case, u is the piezometric head rjJ, stream function 1./J or velocity potential rjJ (see Section 8.4.2). The x and y components of the velocity for the groundwater flow are defined as
where all and an are the permeabilities of the soil along the x and y directions, respectively. ( Examples of each of these field problems will be considered in Section 8.4.
8.3 SOME COMMENTS ON MESH GENERATION AND IMPOSITION OF BOUNDARY CONDITIONS 8.3.1 Discretization of a Domain The representation of a given domain by a collection of finite elements requires engineering judgement on the part of the finite element practitioner. The number, type (e.g., linear or quadratic), shape (e.g., triangular or rectangular), and density (i.e., mesh refinement) of elements used in a given problem depend on a number of considerations. The first is to discretize the domain as closely as possible with elements that are admissible. As we shall see later, one can use one set of elements for the approximation of a domain and another set for the solution. In discretizing a domain, consideration must be given to an accurate representation of the domain, point sources, distributed sources with discontinuities (i.e., sudden change in the intensity of the source), and material and geometric discontinuities, including areentrant corner. The discretization should include, for example, nodes at point sources (so that the point source is accurately lumped at the node), reentrant corners, and element interfaces where abrupt changes in geometry and material properties occur. A second consideration, which requires some engineering judgement, is to discretize the body or portions of the body into sufficiently small elements so that steep gradients of the solution can be accurately calculated. The engineering judgement should come from both a qualitative understanding of the behavior of the solution and an estimate of the computational costs involved in the mesh refinement (i.e., reducing the size of the elements). For example, consider inviscid flow around a cylinder in a
J
SINGLEVARIABLE PROBLEMS
335
Vo FIGURE 8.16 Flow of an inviscid fluid around a cylinder (streamlines).
channel. The flow entering the channel at the left goes around the cylinder and exits the channel at the right (see Fig. 8.16). Since the section at the cylinder is smaller than the inlet section, it is expected that the flow accelerates in the vicinity of the cylinder. On the other hand, the velocity field far from the cylinder (e.g., at the inlet) is essentially uniform. Such knowledge of the qualitative behavior of the flow allows us to employ a coarse mesh (i.e., elements that are relatively large in size) at sites sufficiently far from the cylinder, and a fine one at closer distances to the cylinder (see Fig. 8.17). Another purpose of using a refined mesh near the cylinder is to accurately represent the curved boundary of the domain there. In general, a refined mesh is required in places where acute changes in geometry, boundary conditions, loading, material properties or solution occur. A mesh refinement should meet three conditions: (1) all previous meshes should be contained in the refined mesh; (2) every point in the body can be included within an arbitrarily small element at any stage of the mesh refinement; and (3) the same order of approximation for the solution may be retained through all stages of the refinement process. The last requirement eliminates comparison of two different approximations in two different meshes. When a mesh is refined, care should be taken to avoid elements with very large aspect ratios (i.e., the ratio of the smallest to the largest side of the element)
Yt i,,,,.
.. x FIGURE 8.17 Finite element mesh considerations for inviscid flow around a cylinder. A typical mesh for a quadrant of the domain.
336
FINITE ELEMENT ANALYSIS OF TWO·D1MENSIONAL PROBLEMS
or small angles. Recall from the element matrices (8.41) and (8.68) that the coefficient matrices depend on the ratios of a to band b to a. If the value of a/b or b/a is very large, the resulting coefficient matrices are illconditioned (i.e., numerically not invertible). Although the safe lower and upper limits on b/a are believed to be 0.1 and 10, respectively, the actual values are much more extreme (say, 1000), and they depend on the nature of the physical phenomenon being modeled. For example, in the inviscid flow problem discussed above, large aspect ratios are allowed at the entrance of the channel. The words "coarse" and "fine" are relative. In any given problem, one begins with a finite element mesh that is believed to be adequate (based on experience and engineering judgement) to solve the problem at hand. Then, as a second choice, one selects a mesh that consists of a larger number of elements (and includes the first one as a subset) to solve the problem once again. If there is a significant difference between the two solutions, on~ sees the benefit of mesh refinement, and further refinement may be warranted. If the difference is negligibly small, further refinements are ,not necessary. Such numerical experiments with mesh refinements are not always feasible in practice, mostly because of the computational costs involved. In cases where computational cost is the prime concern, one must depend on one's judgement concerning what is a reasonably good mesh, which is often dictated by the geometry and qualitative understanding of the variations of the solution and its gradient. Since most practical problems are approximated in their engineering formulations, one should not be overly concerned with the numerical accuracy of the solution. A feel for the relative proportions and directions of various errors introduced into the analysis helps the finite element practitioner to make a decision on when to stop refining a mesh. In summary, scientific (or engineering) knowledge and experience with a given class of problems is an essential part of any approximate analysis.
8.3.2
Generation of Finite Element Data
An important part of finite element modeling is mesh generation, which involves numbering the nodes and elements, and the generation of nodal coordinates and the connectivity matrix. While the task of generating such data is quite simple, the type of the data has an effect on the computational efficiency as well as on accuracy. More specifically, the numbering of the nodes directly affects the bandwidth of the final assembled equations, which in turn increases the storage requirement and computational cost if equation solvers with the Gauss elimination procedure are used. The elements can be numbered arbitrarily, because this has no effect on the halfbandwidth. In a generalpurpose program with a preprocessor, options to minimize the baridwidth are included. The saving of computational cost due to a smaller bandwidth in the solution of equations can be substantial, especially in problems where a large number of nodes and degrees of freedom per node are involved. While element numbering does not affect the halfbandwidth, it may
SINGLE·VARIABLE PROBLEMS
337
affect the computer time required to assemble the global coefficient matrixusually, a very small percentage of the time required to solve the equations. The accuracy of the .finite element solution can also depend on the choice of the finite element mesh. For instance, if the selected mesh violates the symmetry of the problem, the resulting solution will be less accurate than one obtained using a mesh that agrees with the physical symmetry of the problem. Geometrically, a triangular element has fewer (or no) lines of symmetry compared with a rectangular element, and therefore one should use meshes of triangular elements with care (e.g., one should select a mesh that does not violate the mathematical symmetry present in the problem). The effect of the finite element meshes shown in Fig. 8.18 on the solution of the Poisson equation in Example 8.3 is investigated. The finite element solutions obtained by the three meshes are compared with the series solution in Table 8.3. Clearly, the solution obtained using mesh 3 is less accurate. This is to be expected, because mesh 3 is symmetric about the diagonal line connecting node 3 to node 7, whereas the mathematical symmetry is about the diagonal line connecting node 1 to node 9 (see Fig. 8.18). Mesh 1 is the most
7~
fY
8
0/0)/
7t
9
8
9
®
5
4
Y
%
6
4
@ (1) 3
L
0""""'"" x
2
3
%
x
(b)
(a)
4 %%1,"'i'l
(c)
FIGURE 8.18 Various types of triangularelement meshes for the domain of Example 8.3: (a) mesh 1; (b) mesh 2; (c) mesh 3.
338
ANITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
TABLE 8.3
Comparison of the finite element solutions obtained using various linear triangularelement meshest with the series solution of the problem in Example 8.3 Finite element solution Node
Mesh 1
Mesh 2
Mesh 3
Series solution
1 2 4 5
0.31250 0.22917 0.22917 0.17708
0.29167 0.20833 0.20833 0.18750
0.25000 0.20833 0.20833 0.16667
0.29469 0.22934 0.22934 0.18114
t See Fig. 8.18 for the finite element meshes.
/'
desirable of the three, because it does not violate the mathematical symmetry of the problem. Next, the effect of mesh refinement with rectangular elements is investigated. Four different meshes of rectangular elements are shown in Fig. 8.19. Each mesh contains the previous mesh as a subset. The mesh shown in Fig. 8.19(c) is nonuniform; it is obtained by subdividing the first two rows and columns of elements of the mesh shown in Fig. 8.19( b). The finite element solutions obtained by these meshes are compared in Table 8.4. The numerical
(a)
(b)
(c)
(d)
FIGURE 8.19 Mesh refinement; the meshes in (a), (b), and (d) are uniform; the mesh in (c) is nonuniform: (a) 2 x 2 mesh; (b) 4 x 4 mesh; (c) 6 x 6 mesh; (d) 8 x 8 mesh.
SINGLE·VARIABLE PROBLEMS
339
TABLE 8.4
Convergence of the finite element solution refinementt) of the problem in Example 8.3 Location
(with
mesh
Finite element solution
x
r
2X2
4x4
6x6
8X8
Series solution
0.0 0.125 0.250 0.375 0.50 0.625 0.750 0.875 0.125 0.250 0.375 0.50 0.625 0.750 0.875
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.125 0.250 0.375 0.50 0.625 0.750 0.875
0.31071
0.29839
0.23220
0.29641 0.29248 0.28055 0.26022 0.23081
0.14137
0.14064
0.29560 0.29167 0.27975 0.24943 0.23005 0.19067 0.14014 0.07709 0.28781 0.26498 0.22873 0.18179 0.12813 0.07332 0.02561
0.29469 0.29077 0.27888 0.25863 0.22934 0.19009 0.13973 0.07687 0.28692 0.26415 0.22799 0.18114 0.12757 0.07282 0.02510
0.28239 0.24107
0.18381
0.28862 0.26580 0.22960 0.18282
0.07506
0.07481
0.26752 0.19286
t See Fig. 8.19 for the finiteelement meshes.
convergence of the finite element solution of the refined meshes to the series solution is apparent from the results presented.
8.3.3 Imposition of Boundary Conditions In most problems of interest, one encounters situations where the portion of the boundary on which natural boundary conditions are specified has points in common with the portion of the boundary on which the essential boundary conditions are specified. In other words, at a few nodal points of the mesh, both the primary and secondary degrees of freedom may be specified. Such points are called singular points. Obviously, one cannot impose boundary conditions on both the primary and secondary variables at the same point. As a general rule, one should impose the essential boundary condition (i.e., the boundary condition on the primary variables) at the singular points and disregard the natural boundary condition (i.e., the boundary condition on the secondary variables), because the essential boundary conditions are often maintained more strictly than the natural ones. Of course, if the true situation in a problem is that the natural boundary conditions are imposed and the essential boundary conditions are a result of this then consideration must be given to the former. Another type of singularity that one encounters in the solution of boundary value problems is the specification of two different values of a primary variable at the same boundary point. An example of such a case is provided by the problem in Fig. 8.20, where u is specified to be zero on the
340
FINITE ELEMENT ANALYSIS OF TWO·DlMENSIONAL PROBLEMS
US=Y 5
I 6
11
16
   
 
21
6
4 x 4 mesh (mesh 1)
(a)
11
16
21
(b)
,   
1 10 ....
73 (c)
73
1 10 ....
8 x 8 mesh (mesh 2)
Cd)
FIGURE 8.20 Effect of specifying (either of the) two values of a primary variable at a boundary node [node 5 in (a) and (b) and node 9 in (c) and (d»).
boundary defined by the line x = 0 (for any y), and to be unity on the boundary defined by the line y = 1 (for any x). Consequently, at x = 0 and y = 1, u has two different values. In the finite element analysis, one must make a choice between the two values, or a weighted average of the two values can be used. When a choice is made between two values, often the larger one is imposed. In any case, the true boundary condition is replaced by an approximate condition. The closeness of the approximate boundary condition to the true one depends on the size of the element containing the point (see Fig. 8.20). A mesh refinement in the vicinity of the singular point often yields an acceptable solution.
8.4 APPLICATIONS 8.4.1 Heat Transfer In Chapter 3, Section 3.3.1, heat transfer (by conduction and convection) in onedimensional (axial and radially symmetric) systems was considered. Here
SINGLE·VARIABLE PROBLEMS
341
we consider heat transfer in twodimensional plane and axisymmetric systems. The derivation of twodimensional heat transfer equations in plane and axisymmetric geometries follows the same procedure as in one dimension, but considers heat transfer in the two directions. Details of such derivations can be found in textbooks on heat transfer (see the references at the end of the chapter). Here we record the governing equations for various cases, construct their finite element models, and present typical applications. For heat conduction in plane or axisymmetric geometries, the finite element models developed in Sections 8.2 and 8.3 are immediately applicable with the following interpretation of the variables:
u = T == temperature (in "C) qn == negative of heat flux (in W m 2 OC 1) all, a22
== conductivities (in W m"
I)
of an orthotropic medium whose principal material axes coincide with the (x, y) axes OC
1== internal heat generation (in W m ?
(8.94)
°C 1)
aoo=O
For convective heat transfer, i.e., when heat is transferred from one medium to the surrounding medium (often a fluid) by convection, the finite element model developed earlier requires some modification. This is because, in twodimensional problems, the convective boundary is a curve as opposed to a point in onedimensional problems. Therefore, the contributions of the convection (or Newtontype) boundary condition to the coefficient matrix and source vector are to be computed by evaluating boundary integrals involving the interpolation functions of elements with convective boundaries. The model to be presented allows the computation of the additional contributions to the coefficient matrix and source vector whenever the element has the convection boundary condition. PLANE SYSTEMS. The governing equation for steadystate heat transfer in plane systems is a special case of (8.1), and is given by
 ~ax (k
x
all ~(k aT\ = I(x, y) a~) ay YayJ
(8.95)
where T is the temperature (in "C), k; and k y are the thermal conductivities (in W m" °C 1) along the x and y directions, respectively, and I is the internal heat generation per unit volume (in W m 3 ) . For a convective boundary, the natural boundary condition is a balance of energy transfer across the boundary due to conduction and/or convection (i.e., Newton's law of cooling): (8.96)
where
fJ is the convective conductance (or the convective heat transfer
342
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
coefficient) (in W m2 "C 1) , Too is the (ambient) temperature of the surrounding fluid medium, and qn is the specified heat flow. The first two terms account for heat transfer by conduction, the third for heat transfer by convection, while the term on the righthand side accounts for the specified heat flux, if any. It is the presence of the term f3(T  Too) that requires some modification of (8.10). The weak form of (8.95) can be obtained from (8.8). The boundary integral should be modified to account for the convective heat transfer term in (8.96). Instead of replacing the coefficient of w in the boundary integral by qn, we use (8.96) (qn is replaced by qn> which is obtained on the element boundary):
0=
i( >.I'
aw sr awa T!) dx dy> Tr w (aT aT) kx+k w kxn +k n ds ax ax Y ay ay rax x y oy y
"
= B(w, T) lew)
(8.97a)
where w is the weight function, and R(·, .) and 1(·) are the bilinear and linear forms aWaT owaj R(w, T)= kx+k dxdy+ f3wTds y >.I' ax ax oy ay r(8. 97b) lew) = w] dx dy + 1. f3wToo ds + 1. wq" ds
i
1( 1
Jr
.I'
Jr.e
The finite element model of (8.58) is obtained by substituting the finite element approximation of the form ~~,/: l ~
,~ &
n
T
= 2:
Tj1jJj(x, y)
(8.98)
;~l
for T and
1jJf for w into (8.97): n
2: (Kij + Hf;)Tj = Ff + r; ;=1
(8.99a)
where
ic:I; = L (kx a1jJi a1jJi a1jJ;) d d a a1jJ; a + kY::l a x y >.I'
Ff =
X
X
c,jy
Y
L f1jJi dx dy + Jr< 1. q~1jJi ds == If + Qf >.I'
(8.99b)
SINGLEVARIABLE PROBLEMS
343
Note that by setting the heat transfer coefficient f3 equal to zero, we obtain the heat conduction model with no account taken of convection. The additional coefficients Hij and Pf due to the convective boundary conditions can be computed by evaluating boundary integrals. These coefficients must be computed only for those elements and boundaries that are subject to the convective boundary condition. The computation of the coefficients for the linear triangular and rectangular elements is presented in the following paragraphs. The coefficients H'f; and P'f for a linear triangular element are defined by "h Ih~1 Hij = f3~2 0 tp'f1jJ'j ds + f3~3 0 1jJ'f1jJ'j ds + f3jl 0 1jJ'ftp'j ds (8.100)
l
.
P'f = f3~2T~2
Ihh
l"h 0
1jJ'f ds
+ f3~ Tz; Jr"~) 1jJ'f ds + f3jl T~l 1"3 0 o
1
1jJ'f ds
where f3ij is the film coefficient (assumed to be constant) for the side connecting nodes i and j of the element g", T% is the ambient temperature on that side, and hij is the length of the side. For a rectangular element, the expressions in (8.100) must be modified to account for four line integrals on four sides of the element.
3
/~
f"\
/"
Boundary of the discretized domain
"
"
....
l
\
l
FIGURE 8.21 Triangular and quadrilateral elements, with node numbers and local coordinates for the evaluation of the boundary integrals. Also shown are the boundary approximation and flux representation using linear and quadratic elements.
344
FINITE ELEMENT ANALYSTS OF TWODIMENSIONAL PROBLEMS
The boundary integrals are line integrals involving the interpolation functions. The local coordinate s is taken along the side, with its origin at the first node of the side (see Fig. 8.21). As noted earlier, the interpolation functions on any given side are the onedimensional interpolation functions. Therefore, the evaluation of integrals is made easy. Indeed, the integrals hfj
L wiwi ds, o
have been evaluated in Chapter 3 in connection with mass matrix coefficients and source vector coefficients for linear and quadratic elements. We summarize the results here. For a linear triangular element, the matrices [He] and {P"] are given by
[2 10] +
[He] = f312 h12 1 2 0 6 0 0 0 e
e
f323 h 23 6 e
e
[0 00] + 0 2 1
0 1 2
[~
"
0
f331 h3; 0 0 6 1 0 e
e
(8.101a)
{P"l =
Ph~!.'hu
Ii} +Ph~'2hh G}
+p;,T
2;!h;,
G}
(8.101b)
For a quadratic triangular element, we have
4 2 1 2 16 2 1 2 4 [He] = f3hhh 0 0 0 30 0 0 0 0 0 0
+ f3S1 h Sl 30
4 0 0 0 1 2
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0
1 4 f32~ T I3h e {pe} = 13 '" 13 1 0 6 0 0
0 0 0
0
0 0 0 + f335 h3S 0 30 0 0
0 0 0 0 0 0 0 0 0 4 2 1 2 2 16 0 2 4 0 1 O· 0 0 0
0 0 0 0 0 0
2 0
0 0 4 2 2 16
0 0
e f3e35 T 3sh ee 35 + 6
(8.102a)
0 0 1 4 1 0
+
W51 T 51heSl cc
6
1 0 0 0 1 4
(8.102b)
SINGLEVARIABLE PROBLEMS
345
For a linear rectangular element, the matrix [He] is of the form 0] {He] = f:3hhh ·12 21 00 00] + f:3~3hb [0 0 20 10 a 6 [ 0000 60120 0000 0000
+ f:3~4h~4
0 a a 0] + [
0 0 0 0 60021 0012
f:3~lh41
[2 001]
6
0 0 0 0 000 0 1002
and {pel is given by
{P'}
~Ph~~hh
{;} +
PbT2~hb {J} + P~~~h~
(8.103a)
mP:,~~'h:,{~} +
(8.103b)
Similar expressions hold for a quadratic rectangular element. AXISYMMETRIC SYSTEMS. For symmetric heat transfer about the z axis
(i.e., independent of the circumferential coordinate), the governing equation is given by
 [!~ T\J = fer, z) ror (rk oT\ a;}+~ az (k a'a;} r
z
(8.104)
where r is the radial coordinate and z is the axial coordinate. The temperature gradient vector is defined by
, (oL
aL)
q = r k, or i + k, az j
and the normal derivative of T (i.e., the negative of the heat flux) across the surface is qn
=
r(kr
oT n; + k; aT n z )
ar
oz
(8.105)
where n, and nz are the direction cosines of the unit normal ft, ft = n)
+nzJ
The weak form of (8.104) is given by
L w{[!~ a1\J  f}rdrdz ror (krr a1\ a;) +~(kz Bz 0; J OW et oW O T ) J: = 2n k, or or + k, oz oz wf r dr dz  2JtYr< wq; ds ( L
0= 2n
Q'
Q'
(8.106)
346
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
where the factors of 2n come from the integration with respect to the circumferential coordinate over (0,2n), and qn is given by (8.105). The convective boundary condition is of the form
q" + rf3(T  To» = qll Substituting for q" = rf3(T  To»
O=2n
oWoT
+ fill
in (8.106), we obtain
oIVoT)
rdrdz2n L (k,;+kzw[ or or OZ OZ Q'
i
w[rf3(TToo)+fi,,]ds
P'
(8.107) The finite element model of this equation is "
(8.108a)
where K':. = 2n ( (k o1J1f 01J1; + k otJIf o1J1i)r dr dz IJ lQ U2 , U3 , Us, U6 , and U7 have the same form as before. The equations for the unknown temperatures (i.e., the condensed equations for the unknown primary variables) are given by
k 
6
4 1 0 1
2 0
1 8 1
0 1
2 2 2
0
2 2 2 2
2 2
2
16
2
0
2
16
1
2
8
2 2
0 8
0
0
fl ! I U2
V3 Us V6 U7
k
=6
0 0 To+V3To 0 To+ To 2To+ V3 V3 To+ To
(8.117)
Their solution is UI = O.6128To,
Us = O.7030To,
U2 = 0.5307To, U6 = O.6088To,
U:,=O.30641Q
U7 = O.3515To
(8.118)
The value of the heat at node 4 is given by
We note that the results obtained using the 3 X 2 mesh of rectangular elements is not as accurate as that obtained with 3 X 2 mesh of triangular elements. This is due to the fact that there are only half as many elements in the former case as in t,he latter. Table 8.5 gives a comparison of the finite element solutions with the analytical solution (8.114a) for two different meshes of linear triangular and rectangular elements, and Fig. 8.23 shows plots of T(x, 0) and qxCx, O)/To, computed using various meshes of
350
ANITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
TABLE 8.5
Comparison of the nodal temperatures T(x, y)/To, obtained using various finite element meshes, t with the analytical solution of (8.109) (Example 8.4) Triangles
Rectangles
x
y
3x2
6x4
3X2
6x4
Analytical solution (8. 114a)
0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0
0.6362
0.6278 0.6064 0.5437 0.4439 0.3139 0.1625 0.7148 0.6904 0.6190 0.5054 0.3574 0.1850
0.6128
0.6219 0.6007 0.5386 0.4398 0.3110 0.1610 0.7102 0.6860 0.6150 0.5022 0.3551 0.1838
0.6249 0.6036 0.5412 0.4419 0.3124 0.1617 0.7125 0.6882/ 0.6171 0.5p38 0.3563 0.1844
1.5 2.0 2.5
0,5510 0.3181 0.7214 0.6248 0.3607
0.5307 0.3064 0.7030 0.6088 0.3515
t See Fig. 8.22 for the geometry and meshes.
triangular and rectangular elements. Note that Qf are heats (in W) whereas q" is the flux (in W m") in the x direction (qx = k aT/ ax). Next we consider an example involving convective heat transfer.
Example 8.5. Consider heat transfer in a rectangular region of dimensions a by b, subject to the boundary conditions shown in Fig. 8.24. We wish to write the finite element algebraic equations for the unknown nodal temperatures and heats. For illustrative purposes, a 4 x 2 mesh of rectangular elements is chosen. We assume that the medium is orthotropic, with conductivities k, and k, in the x and y directions, respectively. It is assumed that there is no internal heat generation. The heat transfer in the region is governed by the equation
 ~ax (k
x
aT\  ~ (k aT\ = 0 8~)
8y
Yay)
The finite element model of the equation is given by [Ke + He]{ue} = {Qe} + {P'} ({fe} = {O})
(8.119)
where uf denotes the temperature at node i of the element Qe. The element matrices are
[Ke]=kxr[~ ~1 :2 2:]+!s:[ 1~ 2~ 6 1 6r
1 1
[He] =
2
{3~iJ
2
[n ! ~] o
0 0 0
2
1
(e = 4, 8)
1
2]
2 1 2 1
1 2
40
Analytical ••·A Triangles  0 Rectangles
60 T(x,O)
To
I
I
A
30
 r
= 100
!
I
Analytical 3 X 2 mesh
To = 100 , k=~tr+__ _..I ..eo.
1·j2~1 l_Ll
qx (x, 0)
'40
I
20
I
Rectangular elements
1O are governed by the Laplace equation, the boundary conditions on them are different in a flow problem, as should be evident from the definitions (8.127) and (8.132). In this section, we consider applications of the finite element method to potential flows, i.e., the solution of (8.128) and (8.133). We consider two examples of fluid flow. The first deals with a groundwater flow problem and the second with the flow around a cylindrical body. In discussing these problems, emphasis is placed on certain modeling aspects, data generation, and postprocessing of solutions. Evaluation of element matrices and assembly is amply illustrated in previous examples, and will not be discussed, since it takes substantial space to write the assembled equations even for the crude meshes used in these examples. Example 8.6 Groundwater flow or seepage. The governing differential equation for a homogeneous aquifer (i.e., one where material properties do not vary withrposition) of unit depth, with flow in the (x, y) plane, is given by
_.i. (a a4» ax 11 ax
.». (a ay
22
a¢) =[ ay
(8.135)
where alI and an are the coefficients of permeability (in m dayf) along the x and y directions, respectively, 4> is the piezometric head (in rn), measured from a reference level (usually the bottom of the aquifer), and [is the rate of pumping (in m3 dayl m"). We know from the previous discussions that the natural and essential boundary conditions associated with (8.135) are Natural
(8.136a)
356
FINITE ELEMENT ANALYSIS OF nVO,DlMENSIONAL PROBLEMS
Essential
(8. 136b)
where f t and f 2 are the portions of the boundaryT of Q such that I', + f 2 = r. Here we consider the following specific problem: find the lines of constant potential 4J (equipotential lines) in a 3000 m X 1500 m rectangular aquifer Q (see Fig. 8.25) bounded on the long sides by an impermeable material (i.e., acpj an = 0) and on the short sides by a constant head of 200 m (o/a = 200 m). In the way of sources, suppose that a river is passing through the aquifer, infiltrating it at a rate of 0.24 m3 day " m2 , and that two pumps are located at (1000, 670) and (1900, 900), pumping at rates Ql = 1200 m3 day"! m t and Q2 = 2400 m3 day" m I , respectively. Pumping is considered as a negative point source. A mesh of 64 triangular elements and 45 nodes is used to model the domain (see Fig. 8.26a). The river forms the interelement boundary bet''een elements (33,35,37,39) and (26,28,30,32). In the mesh selected, neither pump is located at a node. This is done intentionally for the purpose of illustrating the calculation of the generalized forces due to a point source within an element. If the pumps are located at a node then the rate of pumping Qa is input as the specified secondary variable of the node. When a source (or sink) is located at a point other than a node, we must calculate its contribution to the nodes. Similarly, the source components due to the distributed line source (i.e., the river) should be computed. First, consider the line source. We can view the river as a line source of constant intensity, 0.24'm 3 day" m 2 • Since the length of the river is equally divided by nodes 2125 (into four parts), we can compute the contribution of the infiltration of the river at each of the nodes 2125 by evaluating the integrals (see Fig. 8.26b): node 25: node 24:
a 8C'l II
..,;: '0
., '" .c::
...., '" s,
"
E
'"
~
ax
~=o = Uo
ay
D
ct>=o
C
B
~=o ay
FIGURE 8.31 Computational domain and boundary conditions for the velocity potential formulation of inviscid flow around a cylinder (see Fig. 8.28).
SINGLEVARIABLE PROBLEMS
363
boundary conditions there. Along AE, the velocity III = a4>/ax is specified to be Uo. On the surface of the cylinder, the normal velocity, Iln =  84>/an is zero. Thus all boundary conditions, so far, are of the flux type. On the boundary CD, we must know either 4> or 84>/an = a4>/ax. It is clear that 84>/8x = III is not known on CD. Therefore, we assume that 4> is known, and we set it equal to 4>0 = constant. The constant 4>0 is arbitrary and does not contribute to the velocity field (because 84>/ax = III and a4>/8y = lh are independent of 4>0)' It should be noted that specification of 4> at least one point of the mesh is necessary to determine the constant part in the solution for 4> (i.e., eliminate the rigid body motion), We take 4> = 4>0 = 0 on CD. The mathematical boundary conditions of the problem must be translated into finite element data. The boundary conditions on the primary variables are from the boundary CD. We have U2 1 = Un
= U23 = U24 = U25 = 0.0
The only nonzero boundary conditions on the secondary variables come from the boundary AE. There we must evaluate the boundary integral
( aa¢ 1/Ji ds = u;
Jr
at all nodes at which the solution is unknown [because [(x, t) is a known function of time, and the sum of QJ at these nodes is zero]. It should be recalled from Section 6.2.3 that, for different values of a, we obtain different wellknown approximation schemes, with given order of
SINGLEVARIABLE PROBLEMS
375
accuracy:
1' «= {1, 2,
2
3,
the the the the
forward difference scheme (conditionally stable); OeM) CrankNicolson scheme (unconditionally stable); OeM? Galerkin scheme (unconditionally stable); O(M)2 backward difference scheme (unconditionally stable); O(~t)
For the forward difference scheme, the stability requirement is (8.164)
where Amax is the largest eigenvalue of the finite element equations (8.161). We next consider examples of eigenvalue and timedependent problems. Example 8.9 Eigenvalue analysis. Consider the differential equation
au _ (cfu + cfu) =[ at ax ay 2
(8. 165a)
2
in a square region, subject to the boundary conditions
au ax (0, y, t) = 0,
au (x, 0, t) = 0,

ay
u(x, 1, t)
=
0,
u(I, y, t)
=
0
(8. 165b)
and initial condition (8. 165c)
u(x, y, 0) =0
As a first choice, we take a 1 X 1 mesh of triangular (2 elements) and rectangular elements. Alternatively, for the choice of triangles, we can use the diagonal symmetry and model the domain with one triangular element. The triangular element matrices are
[K'] =
~[
1 1 1 2 o 1
:}
[Me]
=~ [: ~ ~] 24 1 1 2
The eigenvalue problem becomes
The boundary conditions require U2 = UJ = O. Hence, we have
(nA + ~)Ul = 0,
or
A= 6
The eigenfunction becomes
U(x, y) =!JII(X, y) = 1 x which is defined over the octant of the domain. For a quadrant of the domain, by symmetry, it becomes Utx, y) = (1 x)(I y).
376
FINITE ELEMENT ANALYSIS OF nVO·DIMENSIONAL PROBLEMS
For one rectangular element, we have
[KO]
ll
1 2 4 1 2 1 4 1 ' 2 1 4
=!6 [: 2 ~1
4 2 1 A
(
2
4
2
2l
4 2 1 1 2 4 2 1 [MO] = 36 1 2 4 2 [
2 124
2l1 +61 [ 14  14
36 1 2 4 2 [ 212 4
2 1
12
Using the boundary conditions Uz = U3 = U4 = 0, we obtain (k)..x4+~)Ul=0,
or
}'=6
The eigenfunction over the quadrant of the domain is given by U(x, Y) = 1J!,(x, y) = (1 x)(l y)
For this problem, the oneelement mesh of triangles in an octant of the domain gives the same solution as the oneelement mesh of rectangular elements in a quadrant of the domain. Table 8.9 gives eigenvalues obtained with various meshes of triangular and rectangular elements, along with the analytical solution of the problem. It is clear that the convergence of the minimum eigenvalue obtained using the finite element method to the analytical value is rapid compared with the convergence of the higher eigenvalues; i.e., the errors in the higher eigenvalues are always larger than that in the minimum eigenvalue. Also, the minimum eigenvalue converges faster with mesh refinements. The mesh used must be such that the required eigenmodes are represented accurately. Example 8.10. We wish to solve the transient heat conduction equation (8.1600)
TABLE 8.9
Comparison of finite element solutions for eigenvalues, obtained using various meshes, with the analytical solution (Example 8.9) Triangles
Rectangles Analyticalt
lXl
A
Al (An) 6.000 ~ (A'13) A3(~I)
A4 (A33) A4 (AI5) A5 (A51)
txt
4x4
8x8
lxl
tx2
4x4
8X8
Am rt
5.415 32.000 38.200 76.390
5.068 27.250 28.920 58.220 85.350 86.790
4.969 25.340 25.730 48.080 69.780 69.830
6.000
5.193 34.290 34.290 63.380
4.999 27.370 27.370 49.740 84.570 84.570
4.951 25.330 25.330 45.710 69.260 69.260
4.935 24.674 24.674 44.413 64.152 64.152
t The analytical solution is
I
Am n =1.11"\111 2 + 112)
(m, 11 = 1,3, 5, ... ).
SINGLEVARIABLE PROBLEMS
377
subject to the boundary conditions, for t "" 0,
aT
aT ay (x, 0, t) = 0
ax (0, y, t) = 0,
(8.l66b)
T(x, 1, t) =0
T(l, y, t) = 0,
and the initial condition
T(x, y, 0) = 0 for all (x, y) in
(8.166c)
Q
We choose a 4 x 4 mesh of linear triangular elements (see Fig. 8.36) to model the domain, and investigate the stability and accuracy of the CrankNicolson method (i.e., a = 0.5) and the forward difference scheme (a = 0.0) for the temporal approximation. Since the CrankNicolson method is unconditionally stable, one can choose any value of !!t. However, for large values of !!t, the solution may not be accurate. The forward difference scheme is conditionally stable; it is stable if !!t < Men where
2
2
M or = Am•• = 386.4 = 0.00518
where the maximum eigenvalue of (8.166a) for the 4 x 4 mesh of triangles is 386.4. The element equations are given by (8.163), with [Me], [Ke], and {fe} defined by (8.159c), wherein c = 1, all = 1, an = 1, ao = 0, and f = 1. The boundary conditions of the problem for the 4 x 4 mesh are given by Us = Uw = U I5 = U20
= U2 1 = U» = Un =
U24 = U25
= 0.0
Beginning with the initial conditions U, = 0 (i = I, 2, ... ,25), we solve the assembled set of equations associated with (8.163). The forward difference scheme would be unstable for !!t > 0.00518. To illustrate this point, the equations are solved using a = 0, M = 0.01 and a = 0.5, M = 0.01. The CrankNicolson method gives a stable and accurate solution, while the forward difference scheme yields an unstable solution (i.e., the solution error grows unboundedly with time), as can be seen from Fig. 8.37. For !!t = 0.005, the forward difference scheme yields a stable solution.
r
T=O
5
1.0.:
Insulated, aT =0 aX
aT
  'iPT= 1
T= 0
at
aT Insulated,  ay
=0
x " 1.0
r
10
15
20
25
4
24
3
23
2
22
I 6
11
16
FIGURE 8.36 Domain, boundary conditions, and finite element mesh for the transient heat conduction problem of Example 8.10.
378
FINITE ELEMEl'IT ANALYSIS OF TWODIMENSIONAL PROBLEMS
0.14 0.12

c0
0.10
···O",Ci.
~
=
0.5, t:.t
=
0.01
= 0.0, t:.t = 0.005
,, , \
= 0.0, lit = 0.01
4 x 4 mesh of linear triangles Mer = 0.00518
E::;'
.,...
Ci.
  Ci.
0.08
~
B ~ ....
8
.,E
f4
0.119
0.06 0.04
/
0.02 0.00 0.00
0.04
0.02
0.06
0.10
0.08
Time t FIGURE 8.37 Stability of the transient solutions of the heat conduction problem in Example 8.10 analyzed using a 4 x 4 mesh of linear triangular elements and the CrankNicolson (a = 0.5) and forward difference (a = 0.0) time integration schemes.
TeO,O.t)
0.3j:;:.,..,..COO""()__i
0.2
2.0 Time (b)
0.0
0.25
0.5
x(aty=O)or
0.75
1.0
y(atx=O)
(0)
FIGURE 8.38 Variation of the temperature as a function of position x and time t for the transient heat conduction problem of Example 8.10.
SINGLE·VARIABLE PROBLEMS
379
TABLE 8.10'
Comparison of finite difference and finite element solutions with the exact solution of the heat conduction problem in Example 8.10
Node 1
2 3 4 5 7
8 9 10 13 14 15 19 20 25
Exact solution
Finile difference solution Error
Finite element solution (steady)
0.2947 0.2789 0.2293 0.1397 0.0000 0.2642 0.2178 0.1333 0.0000 0.1811 0.1127 0.0000 0.0728 0.0000 0.0000
0.2911 0.2755 0.2266 0.1381 0.0000 0.2609 0.2151 0.1317 0.0000 0.1787 0.1110 0.0000 0.0711 0.0000 0.0000
0.3013 0.2805 0.2292 0.1392 0.0000 0.2645 0.2172 0.1327 0.0000 0.1801 0.1117 0.0000 0.0715 0.0000 0.0000
0.0036 0.0034 0.0027 0.0016 0.0000 0.0033 0.0027 0.0016 0.0000 0.0024 0.0017 0.0000 0.0017 0.0000 0.0000
Error 0.0066 0.0016 0.0001 0.0005 0.0000 0.0003 0.0006 . 0.0006 0.0000 0.0010 0.0010 0.0000 0.0013 0.0000
0.0000
Finite element solution (unsteady) at t= 1.0 0.2993 0.2786 0.2278 0.1385 0.0000 0.2628 0.2159 0.1320 0.0000 0.1791 0.1111 0.0000 0.0712 0.0000 0.0000
The CrankNicolson method gives a stable and accurate solution even for t!>.t = 0.05. The temperature T(x, 0, t) is plotted versus x for various values of the time in Fig. 8.38(a). The steady state is reached at time t = 1.0. The temperature T(O, 0, t)
predicted by the CrankNicolson method is plotted versus time in Fig. 8.38(b), which indicates the evolution of the temperature from zero to the steady state. Table 8.10 gives a comparison of the transient solution at t = 1.0 with the steadystate finite element, the finite difference, and the analytical solutions. Table 8.11 gives the finite element solutions for temperature predicted by 4 x 4 meshes of triangles and rectangles and various values of t!>.t and IX = 0.5.
8.5.3 Hyperbolic Equations The transverse motion of a membrane, for example, is governed by a partial differential equation of the form
c
8"u a ( au) a ( 22 au) at2  ax all ax  ay a ay . + aou = f(x, y, t)
(8. 167a)
where u(x, y, t) denotes the transverse deflection, c is the material density of the membrane, a11 and a22 are the tensions in the x and y directions of the membrane, ao is the modulus of the elastic foundation on which the membrane is stretched (often ao = 0, i.e., there is no foundation), and f(x, y, t) is the transversely distributed force. Equation (8.167a) is known as the wave equation, and is classified mathematically as a hyperbolic equation. The function u must be determined such that it satisfies (8.167a) in a region Q,
380
FINITE ELEMENT ANALYSIS OF nVO·D1MENSIONAL PROBLEMS
TABLE 8.11
Comparison of the transient solutions of (8.166) obtained using a mesh of triangular and rectangular elements (a' = 0.5) (Example 8.10) Temperature along the line y =0: Tix, 0, I) Time t
X
10
Element] x=O.O
x =0.25 x=O.5
x =0.75
0.1
T1 Rl 1"2 R2
0.9758 0.9684 0.9928 0.9841
0.9610 0.9556 0.9798 0.9718
0.9063 0.8956 0.9168 0.9020
0.7104 0.6887 0.6415 0.6323
0.2
T1 Rl 1"2 R2
1.8003 1.7723 1.7979 1.7681
1.7238 1.7216 1.7060 1.6990
1.4891 1.4829 1.4644 1.4626
0.9321 0.9367 0.9462 0.9469
T1 R1 T2 R2
2.3130 2.2747 2.2829 2.2479
2.1671 2.1650 2.1448 2.1432
1.7961 1.8084 1.7943 1.8018
1.1469 1.1499 1.1249 1.1319
T1 R1 1"2 R2
2.9960 2.9648 2.9925 2.9621
2.7871 2.8053 2.7862 2.8037
2.2804 2.3090 2.2776 2.3065
1.3843 1.4059 1.3849 1.4053
0.3
1.0
I'
1
tTl, triangular element mesh with tit =0.1; 1"2, triangular clement mesh with tJ.1 = 0.05; RI, rectangular clement mesh with tit = 0.1; Rl, rectangular element mesh with tJ.1 = O. os.
together with the following boundary and initial conditions:
u=
u(x, y,
a
or qn =
0) = uo(x, y),
q" on au
I' (t;;o 0)
(8. 167b)
ai(x, y, 0) = vo(x, y)
(8. 167c)
where a and qn are specified boundary values of u and qn [see (8.155b)], and Uo and Vo are specified initial values of u and its time derivative, respectively. The weak form of (8.167a,b) over a typical element Qe is similar to that of (8.154) [see (8.157)], except that here we have the second time derivative of u:
0=
1[ ( a2uat Q'
av au av au] dx dy>v c 2+aou!) +aU;+a22;ox ax oy ay
i
qnvds
P'
(8.168)
where v = vex, y) is the weight function. The semidiscrete finite element model of (8.167a) is obtained by substituting the finite element approximation (8.158) for u and v = Wj into
SINGLEVARIABLE PROBLeMS
381
(8.168): (8. 169a)
or, in matrix form, (8. 169b)
The coefficients Mij, Kij, and
Ii are
the same as those in (8.159c).
EIGENVALUE ANALYSIS. The problem of finding uit) = ~e~io)l (i = v=I) such that (8.169) holds for homogeneous boundary and initial conditions and I:=: 0 is called an eigenvalue problem of (8.167). We obtain
I (_w2[Me]+[Ke]){ue}:=:{Qe} I
(8.170)
The eigenvalues (j)2 and eigenfunctions L:j=l U/li'j(x, y) are determined from the assembled equations associated with (8.170), after imposing the homogeneous boundary conditions. For a membrane problem, (j) denotes the frequency of natural vibration. The number of eigenvalues of the discrete system (8.170) of the problem is equal to the number of unknown nodal values of U in the mesh. Example 8.11. Consider the free vibrations of a rectangular membrane of homogeneous material, of dimensions a by b (in ft), material density p (in slugs ft2), and fixed on all its edges, i.e., u:::: 0 on T. Although the problem has symmetry about the horizontal centerline and vertical centerlines of the domain (see Fig. 8.39), use of any symmetry in the finite element analysis will eliminate the unsymmetric modes of vibration of the membrane. For example, if we consider a quadrant of the domain in the finite element analysis, the frequencies W mn (m, 11 :f= 1,3,5, ... ) and associated eigenfunctions will be missed in the results [i.e., we can only obtain Wm" (m,lI = 1,3,5, ... )]. By considering the full domain, the first N frequencies allowed by the mesh can be computed, where N is the number of unknown nodal values in the mesh. 
u=o
T 2 It
u= o
1
."
~=o iJy ,q" =0 (a)
(b)
FIGURE 8.39 Geometry, computational domain, finite element mesh, and boundary conditions for the transient analysis of a rectangular membrane with initial deflection: (a) actual geometry; (b) computational domain, finite element mesh of rectangular elements, and boundary conditions (4 x 4 mesh of linear elements or 2 x 2 mesh of ninenode quadratic elements).
382
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
TABLE 8.12
Comparison of natural frequencies computed using various meshes of linear triangular and rectangular elements with the analytical solution of a rectangular membrane fixed on aU its sides (au = an = 12.5, P = T = 2.5) Triangular (linear)
Rectangular (linear)
Wm"
2X2
4x4
8x8
2X2
4X4
8X8
WlI
5.0000
4.2266 5.9083 8.2392 8.3578 10.0618 12.1021 13.2011 14.6942 15.8117
4.0025 5.2068 6.8788 7.5271 8.4565 8.8856 9.9280 11.1193 11.4425
4.3303
4.0285 5.2899 7.2522 7.9527 8.6603 9.9805 12.7157 13.1700 14.0734
3.9522 3.9270 5.0478 4.9673 6.6026 6.3321 7.4200 7.2410 8.0571 7.8540 8.5145 7.85~0 9.1117 8.7810 10.5797 9.4574 10.7280 8.9346
wll W 31
W 12 W 22 W41 W 32 W S1 W 42
Analytical
If only the first eigenvalue W n is of interest or only symmetric frequencies are required, one can use a quadrant of the domain in the analysis. Indeed the results of Example 8.9 are applicable here, with Am n = w~. The results presented in Table 8.9 can be interpreted as the squares of the symmetric natural frequencies of a square a = b = 2) membrane with p = 1 and all = an == T = 1. The exact natural frequencies of a rectangular membrane of dimensions a by b, with tensions all = an = T and density p are
T\ Irl. (m2 n2)1rl.
Wm"
= n"pJ (
7 + b2
(m, n = 1, 2, ... )
To obtain all frequencies, the full domain must be modeled. Table 8.12 contains the first nine frequencies of a rectangular membrane of 4 ft by 2 ft, tension T = 12. 51b fC I, and density p = 2.5 slugs ft 2 , computed using various meshes of linear triangular and rectangular elements. The convergence of the finite element results to the analytical solution is clear. The mesh of linear rectangular element yields more accurate results than the mesh of linear triangular elements.
TRANSIENT ANALYSIS. The hyperbolic equation (8.169b) can be reduced to
a system of algebraic equations by approximating the second time derivative. As discussed in Section 6.2, the Newmark time integration scheme is the most commonly used method, and therefore it is used here. Since the mathematical form of (8.169b) is exactly the same as that in (6.44), the results in (6.47) and (6.49) are immediately applicable to the former. For the sake of convenience, the results are repeated here. The Newmark scheme is
{u ls+l = {u L + M{u}s + ~(At)2{iiL+y
+ M{ii}s+Cl' {ii}sH = (1 8){ii}s + B{ii}s+l {u }s+l = {u}s
(8. 171a)
SINGLEVARIABLE PROBLEMS
383
where, for example, £1'1 z,
y=t
the constantaverage acceleration method (stable)
£1'_1  2,
y=~,
the linear acceleration method (conditionally stable)
aI 2,
.y= 0,
the central difference method (conditionally stable)
(8.171b) The stability criterion is
I ~t 6.t cr ' Figure 8.40 also shows a comparison of the finite element solutions with the analytical solution (8.177). The finite element solutions are in good agreement with the analytical one.
8.6 SUMMARY A stepbystep procedure for finite element formulation of secondorder differential equations in two dimensions with a single dependent variable has been presented. The Poisson equation in two dimensions has been used to illustrate the steps involved. These include weak formulation of the equation, development of the finite element model, derivation of the interpolation functions for linear triangular and rectangular elements, evaluation of element matrices and vectors, assembly of element equations, solution of equations, and postcomputation of the gradient of the solution. A number of illustrative problems of heat transfer (conduction and convection), fluid mechanics and solid mechanics have been discussed. Finally, the eigenvalue and time
SINGLEVARIABLE PROBLEMS
385
0.5,....',    Analytical •••••0 ••• RL4 ,.1, RL2 _._
2 m3 )
I·
a.p
day"?
1
,p
= 200 m
= 1 m day"!
"'" Stream3 (0.5 m m"? dayI)
!
I'
0
(9.45)
From this it is clear that the functions 1; = 1;(x, y) and TJ = TJ(x, y) must be continuous, differentiable,' and invertible. Moreover, the transformation should be algebraically simple so that the Jacobian matrix can be easily evaluated. Transformations of the form (9.34) satisfy these requirements and the requirement that no spurious gaps between elements or overlapping of elements occur. We consider an example to illustrate the invertibility requirements.
426
FlNITE ELEMEl'IT ANALYSIS OF TWODIMENSIONAL PROBLEMS
~ 4
3
A
.'" I fiR ..~2
t t
'l'/  1
=0
1/11 = ~(1  t)(1  7)),
0
everywhere in the element Qe
(9.71)
where [je] is the Jacobian matrix in (9.44b). Geometrically, the Jacobian je represents the ratio of an area element in the real element to the correspond. ing area element in the mas~er element: dA ==dxdy = je d~d7] If je is zero then a nonzero area element in the real element is mapped into
zero area in the master element, which is unacceptable. Also, if :r < 0, a righthanded coordinate system is mapped into a lefthanded coordinate system.
440
FINITE ELEMElIT ANALYSIS OF TWO·DlMENSIONAL PROBLEMS
In general, the Jacobian is a function of £ and 77, implying that the real element is nonuniformly mapped into the master element, i.e., the element is distorted. Excessive distortion of elements is not good, because a nonzero area element can be mapped into a zero or nearly zero area. To ensure :r > 0 and keep within the extreme limits of acceptable distortion, certain geometric shapes of elements must be avoided. For example, the interior angle at each vertex of a triangular element should not be equal to either 00 or 1800 • Indeed, in practice the angle should reasonably be larger than 0 and less than 180 to avoid numerical illconditioning of element matrices. Although the acceptable range depends on the problem, the 0
0
(a)
(b)
(c)
FIGURE 9.18 Finite elements with unacceptable vertex angles; (a) linear quadrilateral elements; (b) linear triangular elements; (c) quadratic quadrilateral and triangular elements. The angles marked are either too small compared with 00 or too large compared with 1800.
INTERPOLATION FUNCTIONS.) NU!liERICAL INTEGRATION AND MODELING
7
1
I
IE
I
..( hl
441
2
1 (a)
4
7
3
T
l~
6
1~ lh] ~
h]
2
i
(b)
FIGURE 9.19 Range of acceptable locations of the midside nodes for quadratic elements: (a) eightnode quadratic element and sixnode quadratic triangular element; (b) the quarterpoint quadrilateral element.
range 15°165° can be used as a guide. Figure 9.18 shows elements with unacceptable vertex angles for straight and curvedsided elements. For higherorder Lagrange elements (also called CO elements), the location of the interior nodes contributes to the element distortion, and therefore they are constrained to lie within a certain distance from the vertex nodes. For example, in the case of a quadratic element, the midside node should be at a distance not less than onequarter of the length of the side from the vertex nodes (see Fig. 9.19). When the midside node is located exactly at a distance of onequarter of the side length from a vertex, the element exhibits special properties (see Problem 9.20). Such elements, called quarterpoint elements, are used in fracture mechanics problems to represent an inverse squareroot singularity in the gradient of the solution at the nearest vertex n o d e . , , . II'})
9.3.3 Mesh Generation
· .•'O)f/ \':.
Generation of a finite element mesh for a given problem should follow the guidelines listed below:
442
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
1. The mesh should represent the geometry of the computational domain and load accurately.' 2. It should be such that large gradients in the solution are adequately represented. 3. It should not contain elements with unacceptable geometries, especially in regions of large gradients. Within the above guidelines, the mesh used can be coarse (i.e., have few elements) or refined (i.e., have many elements), and may consist of one or more orders and types of elements (e.g., linear and quadratic, triangular and quadrilateral). A judicious choice of element order and type could save computational cost while giving accurate results. It should be noted that the choice of elements and mesh is problemdependent. What works well for one problem may not for another. An analyst with physical insight into the process being simulated can make a better choice of elements and mesh for the problem at hand. One should start with a coarse mesh that meets the three requirements listed above, exploit symmetries available in the problem, and evaluate the results thus obtained in the light of physical understanding and approximate analytical and/or experimental information. These results can be used to guide subsequent mesh refinements and analyses. Generation of meshes of a single element type is easy because elements of the same degree are compatible with each other (see Fig. 9.20). Mesh refinements involve several options. The mesh can be refined by subdividing existing elements into two or more elements of the same type (see Fig. 9.21a). This is caned huersion mesh refinement. Alternatively, existing elements can be replaced by elements of higher order (see Fig. 9.21b). This is called pversion mesh refinement. There is also It, puersion mesh refinement, in which elements are subdivided into two or more elements in some places and replaced with higherorder elements in other places. Generally, local mesh refinements should be such that elements of very small size are not placed adjacent to those of very large size (see Fig. 9.22). Combining elements of different kinds naturally arises in solid and structural mechanics problems. For example, plate bending elements (2D) can be connected to a beam element (lD). If the plate element is based on classical plate theory (see Chapter 12), the beam element should be based on the EulerBernoulli beam theory so that they have the same degrees of freedom at the connecting node. When a plane elasticity element (see Chapter 10) is connected to a beam element, which is not compatible with the former in terms of the degrees of freedom at the nodes, one must construct a special element that makes the transition from the 2D plane elasticity element to the 1D beam element (see Problem 10.8). Such an element is called a transition element. Combining elements of different order, say a linear to a quadratic element, may be necessary to accomplish local mesh refinements. There are two ways to do this. One way is to use a transition element that has a different number of nodes on different sides (see Fig. 9.23a). The other way is to
INTERPOLATION FUNCTIONS. NUMERICAL INTEGRATION AND MODELING
443
(a)
(b) FIGURE 9.20 Connecting elements of the same order. The C' elements of the same order ensure the C' continuity along the element interfaces: (a) linear elements; (b) quadratic elements.
!
1r
....!, , r
t
r
{J
,I ~
, I
I
;
I
r
I
I
I
(a)
•
(b)
FIGURE 9.21 The hversion (a) and pversion (b) mesh refinements.
,.: r .I
r
I
I
:
I
f fII
r r
I
I
iiI
:
r
I
tI I
'Tr
!


444
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
(b)
(a)
/
,
•
• (c)
(d)
FIGURE 9.22 Finite element mesh refinements. Meshes shown in (a) and (c) are acceptable, and those shown in (b) and (d) are unacceptable. The typical element of the mesh is shown in each case.
/
I
/
Linear element
/
Transition element
/
/
Quadratic element
(a)
/ I
/
I
Lin ear elements
/
I
Constraint condition (b)
!.
I Quadratic element
FIGURE 9.23 Combining different order elements: (a) use of a transition element that has three sides linear and one side quadratic; (b) use of a linear constraint equation to connect a linear side to a quadratic side.
INTERPOLATION FUNCTIONS, NUMERICAL JNfEGRATION AND MODELING
445
impose a conditionthat constrains the midside node to have the same value as that experienced at the node by the lowerorder element (see Fig. 9.23b). However, such combinations do not enforce interelement continuity of the solution. Figure 9.24 snows element connections that do not satisfy CO continuity along the connecting 'sides. Use of transition elements and constraint conditions in local mesh refinements is a common practice. Figure 9.25 shows a few examples of such refinements. Vex)
T
T
ut = III ...l.
'_'_,IL
Vex)

,
3
U(x)
T
ul = uI
,
ul,=
uj
lS==~~=~s 3
(c)
FIGURE 9.24 Various types of incompatible connections of finite elements. In all cases the interelement continuity of the function is violated along the connecting side.
446
FlNITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
9.3.4
Load Representation
Computation of the nodal contributions of a distributed boundary source was discussed in Chapter 8 [see (8.45)]. The accuracy of the result depends on the element and mesh used to represent the domain. For example, in heat transfer problems, the boundary source is the heat flux across (i.e., normal to) the boundary. Use of linear elements, for example, to represent the boundary will change the actual distribution (see Fig. 9.26). Of course, hversion or pversion mesh refinements will improve the representation of the boundary flux. Another situation where boundary forces are subject to different approximations is the force due to contact between two bodies. For example, a solid plate in contact with a circular disk generates a reactive force that can be represented either as a point load or a locally distributed force. Representation of the contact force between deformable bodies as a point load is. an approximation of the true distribution. A sine distribution might be a more realistic representation of the actual force (see Fig. 9.27).
~/ 1\/ \/
/
cPu ax oy x P ot2 oaxy oay OZv ox + oy + J;, :::: P ot2
(10.3)
where Ix and J;, denote the body forces per unit volume along the x and y directions; respectively, and p is the density of the material. STRAINDISPLACEMENT RELATIONS
2E
xy
au oU oy ax
::::+
(lOA)
STRESSSTRAIN (OR CONSTITUTIVE) RELATIONS
{
ax} :::: [Ctl Oy C 12 0xy
0
12 C
O]{ EX}
C22
0
Ey
0
C66
2Exy
(10.5)
where Cij (Cji:::: Ci) are the elasticity (material) constants for an orthotropic medium with the material principal directions coinciding with the coordinate axes (x, y, z) used to describe the problem. The cil can be expressed in terms
458
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
of the engineering constants (E 1 , E z , V12, and G IZ) for an orthotropic material, for plane stress and plane strain given below (VIZEZ = VZIEl): Plane stress
(10.6)
Plane strain
(VZ3
ell =
c lZ =
= V12)
E I(1 1
VIZ)
,
V12  2v12v21
v 12 E z
1 VIZ 
(to.7) ,
2v 1ZvZl
For an isotropic material, E
G=~
2(1 + v)
(10.8)
BOUNDARY CONDITIONS
Natural (10.9)
Essential
u = u,
v
=D
on I',
(10.10)
where (n.", ny) are the components (or direction cosines) of the unit normal vector il. on the boundary I', r 1 and r z are (disjoint) portions of the boundary (I', and r z do not overlap except for a small number of discrete pointssingular points), tx and ty are specified boundary stresses (or tractions), and u and D are specified displacements. Only one element of each pair, («, tx) and (u, ty ) , may be specified at a boundary point. Equations (10.3) can be expressed in terms of just the displacements u and v by substituting (10.4) into (10.5), and the result into (10.3):
(10.11)
the boundary stress components (or tractions) can also be expressed in terms
PLANE ELASTICITY
459
of the displacements:
(Cll au ax + Cl2 au) ay n, + C66 (au ay + au) ax ny au au) (au au) ty = C66 ( ay + ax nx + Cl ax + C2l oy ny
(. =
(10.12)
l
This completes the review of the governing equations of a plane elastic body undergoing small deformations (i.e., with strains small compared with unity) in the (x y) plane. The material of the body will be assumed in the present study to be linearly orthotropic. J
10.3 WEAK FORMULATIONS 10.3.1 Preliminary Comments Here we study two different ways of constructing the weak form and associated finite element model of the plane elasticity equations (10.11) and (10.12). The first uses the principle of virtual displacements (or the total potential energy principle), expressed in terms of matrices relating displacements to strains, and strains to stresses. This approach is used in most finite element texts on solid mechanics. The second follows a procedure consistent with the previous sections and employs the weak formulation of (10.11) and (10.12) to construct the finite element model. Of course, both methods give, mathematically, the same finite element model, but differ in their algebraic forms.
10.3.2 Principle of Virtual Displacements in Matrix Form First, we rewrite (10.3)(10.5) in matrix form:
_a/ox
[ o
0
a/ay] .. __ { ax} a + {f,x} = p {~} y
a/ay a/ax
«,
(IO.13a)
h
or [T*]{a}
+ {f} =p{~}
(lO.13b)
and
(10. 14a) or
(10.14b)
460
FINITE ElEMENT ANALYSIS OF TWO·DlMENSIONAL PROBLEMS
and {a}
= [C]{E}
(10.15)
where [T] and [T*] are matrices of differential operators. Note that [T] is the transpose of [T*]:
[T] = [T*f =
[a~ax a/nay]
(10.16)
a/ay a/ax Next, we use the (dynamic version of the) principle of virtual displacements [see Reddy (1984)] applied to a plane elastic body:
0=
J.v (O;j fJE;i + pili fJu;) dV _1r~ f; bu, dV _1J.se ti fJUi;fs
(10."17)
where V e denotes the volume of element e [Ve = ge X (~h, !h)], S" is the boundary of ge, he is the thickness of the element, 0 denotes the variational operator, a;; and Eij are the components of stress and strain tensors, respectively, and f; and ii are the components of the body force and boundary stress vectors, respectively: all = ax>
/1 = ix,
012 = 0xy,
/2 =h'
022 = 0y,
t 1 = tx,
etc.
t2 = ty
The first term in (10.17) corresponds to the virtual strain energy stored in the body, the second represents the virtual work done by the body forces, and the third represents the virtual work done by the surface tractions. For plane stress problems with thickness he' we assume that all quantities are independent of the thickness coordinate z. Hence, 0= he
1
(ax
be; + Oy
OEy
+ 20xy fJExy + pii fJu + pv ov) dx dy
Q
K 34 = K54 + Ki2 + K12'
K ll = Kit
etc.
(10.38)
Note ,that K l 2 , for example, denotes the coupling stiffness coefficient between the first and second global degrees of freedom, both of which are at global node 1. Similar arguments apply for the assembly of the force vector. With regard to the boundary conditions, the primary degrees of freedom are the displacements and the secondary degrees of freedom are the forces. At
PLANE ELASTICITY
469
any node on the boundary of the domain of a plane elasticity problem, one has the following four distinct possibilities: Case 1 Case 2 Case 3 Case 4
u and v are specified (and t; and ty are unknown) u and ty are specified (and tx and v are unknown)
tx and v are specified (and u and ty are unknown) tx and ty are specified (and u and v are unknown)
(10.39)
In general, only one of the quantities of each of the pairs (u, tJ and (v, ty) must be specified at a point on the boundary. For timedependent problems, the initial displacement and velocity must be specified: (10.40)
10.7 EXAMPLES Here we consider a couple of examples of plane elasticity problems to illustrate the load computation and imposition of boundary conditions. The stresses are evaluated at reduced Gauss points of the elements [see Barlow (1976, 1989)]. These examples are actually analyzed using the program FEM2DV2, which is discussed in Chapter 13. Example 10.1. Consider a thin elastic plate subjected to a uniformly distributed edge load, as shown in Fig. 10.5(a). We wish to determine the equilibrium (i.e., static) solution. First we consider a twoelement discretization of the plate by triangular elements, and then we perform all the algebra to obtain the nodal displacements. The assembly of element matrices for elements with two degrees of freedom (DOF) is described in Section 10.6. For the finite element mesh at hand, the correspondence between the global and local nodes and stiffness is as follows:
Nodal correspondence Global (DOF)
Local (DOF)
Stiffness correspondence Global
Local
K ll
1 of element 2 (1,2)
K 12
K1I+K~1 Kh+Kh Kl 2+ K~2
2 (3, 4)
2 of element 1 (3,4)
K 33 K 44 K 34
1 K 44 K~4
3 (5, 6)
3 of element 1 (5,6) 2 of element 2 (3,4)
4 (7,8)
3 of element 2 (5, 6)
1 of element 1 (1,2)
Kn
1 (1,2)
Kj3
K 55 K 66 K 56
K~5+K~3
Kn
K~5
K ss K7 8
K~+K1 K~6+K~ K 266
K~6
(10.41)
470
FINITE ELEMENT ANALYSIS OF lWO·D1MENSIONAL PROBLEMS
CD I
I
2
/l a = 120 in  ...j2 (a)
y
y
160 in
160 in 0::0
/
o+x 120 in
(b)
120 in
(c)
FIGURE 10.5 Geometry and finite element mesh of a plane elasticity problem by the CST elements: (a) a plane elasticity problem; (b) element I; (c) elemenl2.
If two global nodes correspond to two (local) nodes of the same element then the
corresponding stiffness coefficient is nonzero; otherwise it is zero. The specified (primary) degrees of freedom are (10.42)
The known secondary degrees of freedom (forces) are already included in {PI} and {p2}. Equilibrium of forces requires that (note that Po is the force per unit length) (10.43)
The first two rows and columns and the last two rows and columns of the assembled [K] can be deleted (since the specified boundary conditions are homogeneous) to obtain the condensed form of the matrix equation:
(10.44)
or (using a = 120 in, b
= 160 in, h = 0.036 in, v = 0.25,
E = 30 X 106lb in ", and Po =
PLANE ELASTICITY
93.0 36.0 16.2 72.0 21.6 36.0 16.2 21.6 43.2 0.0
104
114.4
·93.0
14.4]{ U2} Y2 = {800' 0.00}
43.2 0.0 72.0
U3
V;
471
(10.45)
800.0 0.0
Inverting the matrix, we obtain
U2} V; = 10{
U3 V3
6
l4.07 2.8.6534 2.34
3
=
1O4 {
0.17 0.59
0.17 1.6 3.63 1.6 4.72 0.99
1~:~~}
10.10 1.09
0'59]{800'0.00} 8.00 6.88 0.0
4.72 0.99
in
(10.46)
Note that the solution is not symmetric (i.e., U2 =I U3 , Y2 =I  V3) about the horizontal centerline. The exact solution should be symmetric about the centerline. This is due to the asymmetry in the finite element discretization. Table 10.1 gives the finite element solutions (obtained using FEM2DV2) for the displacements at the points (120,0) and (120, 160) of isotropic and orthotropic plates
TABLE 10.1
Finite element results for a thin plate (plane stress assumption) using various meshes of triangular and rectangular elements and material propertiest Mesh
Material
U2
V2
(x 10 4)
U3 (x 10 4 )
V3 (X
11.291 10.853
1.9637 2.3256
10.113 10.853
2.3256
1.6662 2.6758
10.650 10.728
1.579 2.6758
2.157 2.002
10.821 10.778
2.157
(x 10 4)
[zJ
Isotropic: E = 30 x 111' Ib in 2 v = 0.25 G = E/2(l + v)
0
Orthotropic: E I = 31 x uti Ib in 2 E 2 = 2.7 x 111' Ib in 2 10.767 G 12 = 0.75 x Hf Ib in 2 10.728 Vl2 = 0.28
1
1
~1.08
2
2
~ 1
10 4)
Orthotropic: . same as above
10.821 10.778
HlO2
2
For each mesh, the firstrow corresponds to triangular elements and the second row to one rectangular element.
472
FINITE ELEMENT ANALYSIS OF TWO·DlMENSIONAL PROBLEMS
TABLE 10.2
Deflections and stresses in an isotropic plate subjected to uniform edge load (Example 10.1) u(120,0) (x 10 4 )
u(120,0) (x 10 4 )
1x t
11.291
1.9637
2x2
11.372
2.1745
4x4
11.284
2.1255
1x 1
10.853
2.3256
2x2
11.078
2.0212
4x4
11.150
2.0094
Mesh
Triangular
Rectangular
a..
ay
a xy
285.9 (80, 53.33)t 279.7 (80,53.33) 280.0 (80,53.33)
14.40 (80,53.33) 69.36 (20,53.33) 69.59 (10,26.67)
10.80 (80,53.33) 23.20 (40,26.67) 35.93 (20,13.33)
277.8 (60,80) 277.8 (30,40) 280.4 (75,60)
25.84 (60,80) 37.46 (30,40) , 49.74 (15,60)
0.0 / (60,80) 13.23 (30,40) 27.73 (15,20)
t Location of the stress. y
'f
= 150 psi
r '1 !.
a = 10 in
»
~xDEC~2i" .j
b (a)
75
,/".111 12 __ 13 ___14 __ ~. i i 6E::I==:E~~t=:»>
10
1
~
2
3
4
(b)
FIGURE 10.6 Finiteelement meshes for an endloaded cantilever beam.
5
= 1 in
"c
PLANE ELASTICITY
473
for the meshes shown. Table 10.2 gives the deflections and stresses obtained with refined meshes of triangular and rectangular elements. Example 10.2. Consider the' cantilever beam (E = 30 x 106 psi, v = 0.25, a = 10 in, b = c = 1 in) shown in Fig. 10.6(a). 'We wish to determine, using the elasticity equations, the maximum deflection and bending stress in the beam when it is subjected to a uniformly distributed shear stress l' = 150 psi. The boundary conditions of the problem are uta, y) = 0,
t,
(" = 0,
= ty = 0 f y = 1' ty
=
v(a, c)
=
0
= 0, 2c for any x x = 0 for any y
at y at
0 at x = a and for any y
(10.47)
*' c
We shall solve the problem using the plane stress assumption. The elastic coefficients clj for the plane stress case are defined (assuming that steel is isotropic) by E Cll =CI2=1_
y'
C12 =
Ev 1 v2>
C 66
E = 2(1 + v) (= G)
(10.48)
Three different finite element meshes, increasingly refined, are shown in Fig. 1O.6(b). These meshes are those consisting of linear rectangular elements. Equivalent triangular element meshes are obtained by joining node 1 to node 3 of each rectangular element (see the dashed lines). Equivalent meshes of ninenode quadratic Lagrange elements are obtained by considering a 2 x 2 mesh of linear Lagrange elements as a quadratic element. For the finite element model, the boundary conditions on the primary and secondary variables, e.g., for the ISnode mesh, are given by U5 = UIO = Vw = U15 = 0.0 F{= 75.0,
F~=
150.0,
F{1 = 75.0
(10.49)
and all other forces are zero on the boundary. Table 10.3 gives a comparison of the finite element solutions with the elasticity solutions for the tip deflection (i.e., the deflection at the center node of the left end) and bending stress 0", obtained using twodimensional elasticity theory [see Reddy (1984)J. The linear triangular element mesh has the slowest convergence compared with the linear and quadratic rectangular elements.
The final example in this chapter deals with free vibration and transient analysis of a plane elasticity problem. We consider the cantilever beam of Example 10.2. Example 10.3. Consider the cantilever beam shown in Fig. 10.6a. We wish to determine the natural frequencies and transient response using the plane elasticity elements. We use the finite element meshes of linear triangular and rectangular elements shown in Fig. 1O.6(b) and their nodal equivalent meshes of quadratic elements to analyze the problem. Table lOA gives a comparison of the first 10 natural frequencies obtained with various meshes. The convergence of the natural frequencies with mesh refinement is clear.
474
FINITE ELEMENT ANALYSIS OF TWO·DlMENSIONAL PROBLEMS
TABLE 10.3
Comparison of the finite element solution with the elasticity solution for a cantilever beam subjected to a uniform shear load at the free end Tip deflection, u
X
10 2
Normal stress ax Linear rectangles
Number of nodes
Linear triangles
Linear rectangles
Quadratic rectangles
Linear triangles
15
0.16112
0.31335
0.50310
27
0.26623
0.43884
0.51288
51
0.316459
0.48779
0.51374
1196 1209 (8.75,1.5):1: (15)t 2270 1793 (31) (9.375, 1.5) 2829 2056 (63) (9.6875,1.5) 2876 2180 (9.583, 1.667) (9.6875,1.5)
Elasticityj
0.51875 (0.0. 1.0)
t Element number.
Quadratic rectangles 2196 (8.943. 1.577) 2439 (9.471, 1.577) 2526 (9.736,1.577) 2528 "(9.736, 1.577)
,
:j:Quadrature points. 3 2 2 § From Reddy (1984), p. 53: u(O, 1) = (PL /3£/)[1 + 3(1 + v)/ L ]; ax = Px(1  y)/l, 1=3'
For transient analysis with the linear acceleration scheme (a = 0.5, Y = time step /1t is given by the inequality At for example, will be replaced by the ith interpolation function used in the approximation of u. Similarly, W 3 will be replaced by the ith interpolation function used in the approximation of P. When different interpolations are used for (u, v) and P, this interpretation becomes necessary. The second step of the weak formulation for the present case needs some comments. The integration by parts to distribute differentiation equally among the variables (u, v, P) and weight functions (WI> IV2, IV3) helps relax the continuity requirements on the finite element approximation functions used for u, v, and P. However, in any problem, such trading of differentiability is subject to the restriction that the resulting boundary expressions must be physically meaningful. Otherwise, the secondary variables of the formulation may not be the quantities the physical problem admits as boundary conditions. Examination of the boundary stress components t, and ty shows that the pressure term is a part of them [see (11.3c)]. For example, tx is the x component of the total boundary stress, which is the sum of the viscous boundary stress, 2jl(au/ax)nx + jl(aujay + av/ox)ny , and the hydrostatic boundary stress,  Pn.. Therefore, it is necessary that each term, except the body force term, must be integrated by parts in the two momentum equations. By trading the differentiation from Ponto' WI and W2 in the two momentum equations, we gain both physically meaningful natural boundary conditions and the symmetry of the finite element equations, as we shall see shortly. No integration by parts is used in the continuity equation, because no relaxation of continuity on u and v can be accomplished; further, the resulting boundary conditions are not physically justifiable.
486
FINITE ELEMENT ANALYSIS OF TWODIMENSIONAL PROBLEMS
Keeping the above comments in mind, we carry out the remaining two steps of the weak formulation. The result is 0=
au
OWl au
aWl (au OV)
at ax ax ay +ay ax L [pWI+2jl+f.lQ'
 aWl ax Pwdx ] dxdy  1. Wl tx ds
X
0=
av awz (au au) au  + 2awz flat ax oy+ax oy ay L [PWz+fl
(11.7)
Q'
 awz ay P  wzh ] dx dy  1JI" wzty ds 0= 
au av) L (ax+oy dx dy W3
Q'
Note that there is no boundary integral involving W3, because no integration by parts is used. This implies that P is not a primary variable; it is part of the secondary variables (tx and ty ) . This in turn requires that P not be made continuous across interelement boundaries. If P by itself is not specified in a problem (but Lx and ty are specified) then P is arbitrarily set equal to a value at some node to determine the constant part of the pressure. Thus, P can be determined only to within an arbitrary constant. The minus sign in the third statement is inserted because P ~  W3, which makes the resulting finite element model symmetric. An examination of the weak form reveals that u and v are the primary variables that should be made continuous at interelement boundaries, while P is a nodal variable that is not t~be made continuous across the interelement boundaries. Therefore, the Lagrange family of finite elements can be used for (u, v, P). The weak form shows that the minimum continuity requirements on (u, v, P) are
(u, v) linear in x and y P constant
Thus, there are different continuity requirements on the interpolation of the velocity field and pressure: n
U=
2: Uj1J.!'j, j=1
n
m
V = 2: Vj1/JJ, j=1
P = ~ Ej} r .,T , I
1 .}
I
I
I
i
'"t
_ ..... _...J_ ...........
I
;
I
I
Ix =
0
ty = 0
r fI _ _ x
U=lx=O
i o."
e
(b)
FIGURE 11.3 Boundary conditions and finite element meshes for the problem of an incompressible viscous fluid squeezed between plates: (a) geometry of the flow domain; (b) finite element mesh of a quadrant.
8c '"
§ ;g"" ~
63 r
rrJ
E 6 '"
~
496
FINITE ELEMENT ANALYSIS OF TWO· DIMENSIONAL PROBLEMS
large compared with both their width and the distance between them, we have a case of plane flow (in the plane formed by the plate width and the distance between them). Assurnining that a state of plane flow exists, we determine the velocity and pressure fields. Let "Va be the velocity with which the two plates are moving toward each other (i.e., squeezing out the fluid), and let 2b and 2a denote, respectively, the distance between and the width of the plates (see Fig. 11.3a). An approximate analytical solution to this twodimensional problem is provided by Nadai (1963): 3 VoX ( y2) u=v;1 b2
v = _ Voy 2b
J
(3 _b"y2)
3~Voa2
PO=2j;3
ax == 2tt au  P = Po2 (x 2  3y 2  a2 + 2 b2)
ax
av
a == 2~ 
ay 
y
r xy
(ll.37)
a
2 2 2 2) P = Po( 2 x + y  a  2b
a
av)
== ~(aU + ay ax
= 
2PoXY 2 a
Owing to the biaxial symmetry of the problem, it suffices to model only a quadrant of the domain. A 5 X 3 nonuniform mesh of ninenode quadratic elements is used in the velocitypressure model, and a 10 x 6 mesh of fournode linear elements and a 5 x 3 mesh of quadratic elements is used in the penalty model (see Fig. 11.3b). The nonuniform mesh, with smaller elements near the free surface (i.e., at x = a), is used to accurately approximate the singularity in the shear stress at the point (a, b) = (6, 2). The mesh used for the penalty model has exactly the same number of nodes as the mesh used for the velocitypressure model. There are no specified TABLE 11.2
Comparison of the horizontal velocity u(x, 0) obtained using the penalty modelt and the mixed (velocitypressure) model with the analytical solution
x
4·node
9·node 4·node
Mixed Mixed model model 9node 4·node 9·node 8·node 9·node
1 2 3 4 4.5 5.0 5.25 5.? 5.75 6.0
0.0303 0.0677 0.1213 0.2040 0.2611 0.3297 0.3674 0.4060 0.4438 0.4793
0.0310 0.0691 0.1233 0.2061 0.2631 0.3310 0.3684 0.4064 0.4443 0.4797
0.6513 1.3062 1.9769 2.6730 3.0463 3.3956 3.5732 3.6874 3.7924 3.7862
Y= 1.0
Y= 100
0.6563 1.3165 1.9911 2.6960 3.0718 3.4347 3.6120 3.7388 3.8316 3.8362
y=lOS
0.7576 1.5135 2.2756 3.0541 3.4648 3.8517 4.0441 4.1712 4.2654 4.2549
0.7505 1.4992 2.2557 3.0238 3.4307 3.8029 3.9944 4.1085 4.2160 4.1937
0.7496 1.5038 2.2563 3.0213 3.4331 3.8249 4.0074 4.1450 4.2188 4.2659
0.7497 1.5031 2.2561 3.0203 3.4292 3.8165 3.9893 4.1204 4.2058 4.2364
Analytical solution 0.7500 1.5000 2.2500 3.0000 3.3750 3.7500 3.9375 4.1250 4.3125 4.5000
t The threepoint Gauss rule for nonpenalty terms and the twopoint Gauss rule for penalty terms are used for quadratic elements.
•
FLOWS OF VISCOUS INCOMPRESSIBLE FLUIDS
497
nonzero secondary variables in the problem. The velocities u(x, 0) obtained using the two models are compared with the analytical solution in Table 11.2. The influence of the penalty parameter on the accuracy of the solution is clear from the results. Next, a 12 X 8 mesh of linear elements and a 6 X 4 mesh of quadratic elements are used to evaluate the relative accuracies of the rectangular and triangular elements. The 12 X 8 mesh of linear triangular elements with full integration of [K 1] and [r] [see (11.35)] and selective integration (i.c., full integration of [K 1] and reduced integration of [K 2 ]) both give the same results for the velocity field. However, in both cases, erroneous results for pressure and stresses are obtained. The 6 X 4 mesh of quadratic triangular elements with full and selective integrations gives the same velocity fields, while the stresses and pressure are predicted to be the same at the same quadrature points. Both the 12 X 8 mesh of linear rectangular elements and the 6 X 4 mesh of ninenode rectangular elements give good results for velocities, pressure, and stresses. Figure 11.4 shows plots of the velocity u(x, y) for x = 4 and 6, and Fig. 11.5 shows plots of pressure Pix, y), for y = constant, computed using the quadratic triangular and ninenode rectangular elements. The pressure in the penalty model was computed using (11.34) with the 2 X 2 Gauss rule for the quadratic rectangular element and the onepoint formula for the quadratic triangular element. If the pressure in the penalty model were computed using the full quadrature rule for rectangular elements, we should obtain erroneous values. The linear triangular element with full as well as reduced integrations gives unstable pressures, while the quadratic triangular element with one or twopoint rules yields good results. In general, the same quadrature rule as that used for the evaluation of the penalty terms in the coefficient matrix must be used to evaluate the pressure, and one should avoid using the linear triangular element.
The next example deals with the finite element analysis of a lubrication problem. Example 11.3 Flow of a viscous lubricant in a slider bearing. The slider (or slipper) bearing consists of a short sliding pad moving at a velocity u = Un relative to a stationary pad inclined at a small angle with respect to the stationary pad, and the small gap between the two pads is filled with a lubricant (see Fig. 11.6a). Since the ends of the bearing are generally open, the pressure Pothere is atmospheric. If the upper pad is parallel to the base plate, the pressure everywhere in the gap must be atmospheric y
y
2
2
01.,,..0
o
1
2
Velocity u(4, y)
2
3
4
5
Velocity u( 6, y)
FIGURE 11.4 Horizontal velocity distributions at x = 3 and 6 for the problem of a viscous incompressible fluid squeezed between parallel plates (the solid lines represent the analytical solutions and the circles finite element results with all elements for the meshes used). .
~
00
:!l
:z
~
m
81 0 .0....
8
00,
"t:J...
""
;:;
~ 6
o
\,
""'" ~,,o.. ....... ~ \ '0\
o
~
,
Analytical o FEM (R9) (y '" 0.1057)\
~
OJ
~
2
o
  Analytical 0 FEM (R9) (y
!
o
I
=;
1.8943)
I
~
'f.
~
I
2
3 Distance x
(a)
\ 
~
\
~ 2
~
Analytical FEM (T6) (y '" ~ )
o
O.l
.\ , i ' 6 5 4
I
~ \,
d I:l..
I
,
C;;
o"rl
,
>:::4
I
\\~' \,: ,
~
::l
\.
~
:
~
'',
o'\.. '