Preface Finite element procedures are now an important and frequently indispensable part of engineering analysis and des
249 73 44MB
English Pages xiv,1037 [1050] Year 1996
Table of contents :
Preface xiii
CHAPTER ONE —
An Introduction to the Use of Finite Element Procedures 1
1.1 Introduction 1
1.2 Physical Problems, Mathematical Models, and the Finite Element Solution 2
1.3 Finite Element Analysis as an Integral Part of ComputerAided Design 11
1.4 A Proposal on How to Study Finite Element Methods 14
CHAPTER TWO
Vectors, Matrices, and Tensors 17
2.1 Introduction 17
22 Introduction to Matrices 18
23 Vector Spaces 34
24 Definition of Tensors 40 ⋅
2.5 The Symmetric Eigenproblem Av = Av 51
2.6 The Rayleigh Quotient and the Minimax Characterization
of Eigenvalues 60
2.7 Vector and Matrix Norms 66
2.8 Exercises 72
CHAPTER TRHEE
Some Basic Concepts of Engineering Analysis and an Introduction 77
to the Finite Element Method 77
3.1 Introduction 77
3.2 Solution of DiscreteSystem Mathematical Models 78
3.2.1 SteadyState Problems, 78
3.2.2 Propagation Problems, 87
3.2.3 FEigenvalue Problems, 90
3.2.4 On the Nature of Solutions, 96
3.2.5 Exercises, 101
33 Solution of ContinuousSystem Mathematical Models 105
3.3.1 Differential Formulation, 105
3.3.2 Variational Formulations, 110
3.3.3 Weighted Residual Methods; Ritz Method, 116
3.3.4 An Overview: The Differential and Galerkin Formulations, the Principle of
Virtual Displacements, and an Introduction to the Finite Element Solution, 124
3.3.5 Finite Difference Differential and Energy Methods, 129
3.3.6 Exercises, 138
3.4 Imposition of Constraints 143
3.4.1 An Introduction to Lagrange Multiplier and Penalty Methods, 143
3.4.2 Exercises, 146
CHAPTER FOUR
Formulation of the Finite Element Method—Linear Analysis in Solid and Structural Mechanics 148
4.1 Introduction 148
4.2 Formulation of the DisplacementBased Finite Element Method 149
4.2.1 General Derivation of Finite Element Equilibrium Equations, 153
4.2.2 Imposition of Displacement Boundary Conditions, 187
4.2.3 Generalized Coordinate Models for Specific Problems, 193
4.2.4 Lumping of Structure Properties and Loads, 212
4.2.5 Exercises, 214
4.3 Convergence of Analysis Results 225
4.3.1 The Model Problem and a Definition of Convergence, 225
4.3.2 Criteria for Monotonic Convergence, 229
4.3.3 The Monotonically Convergent Finite Element Solution: A Ritz Solution, 234
4.3.4 Properties of the Finite Element Solution, 236
4.3.5 Rate of Convergence, 244
4.3.6 Calculation of Stresses and the Assessment of Error, 254
4.3.7 Exercises, 259
4.4 Incompatible and Mixed Finite Element Models 261
4.4.1 Incompatible DisplacementBased Models, 262
4.4.2 Mixed Formulations, 268
4.4.3 Mixed Interpolation— Displacement/Pressure Formulations for
Incompressible Analysis, 276
4.4.4 Exercises, 296
4.5 The InfSup Condition for Analysis of Incompressible Media and Structural
Problems 300
4.5.1 The InfSup Condition Derived from Convergence Considerations, 301
4.5.2 The InfSup Condition Derived from the Matrix Equations, 312
4.5.3 The Constant (Physical) Pressure Mode, 315
4.5.4 Spurious Pressure Modes—The Case of Total Incompressibility, 316
4.5.5 Spurious Pressure Modes—The Case of Near Incompressibility, 318
4.5.6 The InfSup Test, 322
4.5.7 An Application to Structural Elements: The Isoparametric Beam Elements, 330
4.5.8 Exercises, 335
CHAPTER FIVE
Formulationand Calculation of Isoparametric Finite Element Matrices 338
5.1 Introduction 338
5.2 Isoparametric Derivation of Bar Element Stiffness Matrix 339
53 Formulation of Continuum Elements 341
5.3.1 Quadrilateral Elements, 342
5.3.2 Triangular Elements, 363
5.3.3 Convergence Considerations, 376
5.3.4 Element Matrices in Global Coordinate System, 386
5.3.5 Displacement/Pressure Based Elements for Incompressible Media, 388
5.3.6 Exercises, 389
5.4 Formulation of Structural Elements 397
5.4.1 Beam and Axisymmetric Shell Elements, 399
5.4.2 Plate and General Shell Elements, 420
5.4.3 Exercises, 450
55 Numerical Integration 455
5.5.1 Interpolation Using a Polynomial, 456
5.5.2 The NewtonCotes Formulas (OneDimensional Integration), 457
5.5.3 The Gauss Formulas (OneDimensional Integration), 461
5.5.4 Integrations in Two and Three Dimensions, 464
5.5.5 Appropriate Order of Numerical Integration, 465
5.5.6 Reduced and Selective Integration, 476
5.5.7 Exercises, 478
5.6 Computer Program Implementation of Isoparametric Finite Elements 480
CHAPTER SIX
Finite Element Nonlinear Analysis in Solid and Structural Mechanics 485
6.1 Introduction to Nonlinear Analysis 485
6.2 Formulation of the Continuum Mechanics Incremental Equations of
Motion 497
6.2.1 The Basic Problem, 498
6.2.2 The Deformation Gradient, Strain, and Stress Tensors, 502
6.2.3 Continuum Mechanics Incremental Total and Updated Lagrangian
Formulations, Materially NonlinearOnly Analysis, 522
6.2.4 Exercises, 529
6.3 DisplacementBased Isoparametric Continuum Finite Elements 538
6.3.1 Linearization of the Principle of Virtual Work with Respect to Finite Element
Variables, 538
6.3.2 General Matrix Equations of DisplacementBased Continuum Elements, 540
6.3.3 Truss and Cable Elements, 543
6.3.4 TwoDimensional Axisymmetric, Plane Strain, and Plane Stress Elements, 549
6.3.5 ThreeDimensional Solid Elements, 555
6.3.6 Exercises, 557
6.4 Displacement/Pressure Formulations for Large Deformations 561
6.4.1 Total Lagrangian Formulation, 561
6.4.2 Updated Lagrangian Formulation, 565
6.4.3 Exercises, 566
6.5 Structural Elements 568
6.5.1 Beam and Axisymmetric Shell Elements, 568
6.5.2 Plate and General Shell Elements, 575
6.5.3 Exercises, 578
6.6 Use of Constitutive Relations 581
6.6.1 Elastic Material Behavior—Generalization of Hooke’s Law, 583
6.6.2 Rubberlike Material Behavior, 592
6.6.3 Inelastic Material Behavior; Elastoplasticity, Creep, and Viscoplasticity, 595
6.6.4 Large Strain Elastoplasticity, 612
6.6.5 Exercises, 617
6.7 Contact Conditions 622
6.7.1 Continuum Mechanics Equations, 622
6.7.2 A Solution Approach for Contact Problems: The Constraint Function Method, 626
6.7.3 Exercises, 628
6.8 Some Practical Considerations 628
6.8.1 The General Approach to Nonlinear Analysis, 629
6.8.2 Collapse and Buckling Analyses, 630
6.8.3 The Effects of Element Distortions, 636
6.8.4 The Effects of Order of Numerical Integration, 637
6.8.5 Exercises, 640
CHAPTER SEVEN
Finite Element Analysis of Heat Transfer, Field Problems, and Incompressible Fluid Flows 642
7.1 Introduction 642
7.2 Heat Transfer Analysis 642
7.2.1 Governing Heat Transfer Equations, 642
7.2.2 Incremental Equations, 646
7.2.3 Finite Element Discretization of Heat Transfer Equations, 651
7.2.4 Exercises, 659
7.3 Analysis of Field Problems 661
7.3.1 Seepage, 662
7.3.2 Incompressible Inviscid Flow, 663
7.3.3 Torsion, 664
7.3.4 Acoustic Fluid, 666
7.3.5 Exercises, 670
1.4 Analysis of Viscous Incompressible Fluid Flows 671
7.4.1 Continuum Mechanics Equations, 675
7.4.2 Finite Element Governing Equations, 677
7.4.3 High Reynolds and High Peclet Number Flows, 682
7.4.4 Exercises, 691
CHAPTER EIGHT
Solution of Equilibrium Equations in Static Analysis 695
8.1 Introduction 695
8.2 Direct Solutions Using Algorithms Based on Gauss Elimination 696
8.2.1 Introduction to Gauss Elimination, 697
8.2.2 The LDLT Solution, 705
8.2.3 Computer Implementation of Gauss Elimination—The Active Column Solution, 708
8.2.4 Cholesky Factorization, Static Condensation, Substructures, and Frontal
Solution, 717
8.2.5 Positive Definiteness, Positive Semidefiniteness, and the Sturm Sequence
Property, 726
8.2.6 Solution Errors, 734
8.2.7 Exercises, 741
83 Iterative Solution Methods 745
8.3.1 The GaussSeidel Method, 747
8.3.2 Conjugate Gradient Method with Preconditioning, 749
8.3.3 Exercises, 752
8.4 Solution of Nonlinear Equations 754
8.4.1 NewtonRaphson Schemes, 755
8.4.2 The BFGS Method, 759
8.4.3 LoadDisplacementConstraint Methods, 761
84.4 Convergence Criteria, 764
8.4.5 Exercises, 765
CHAPTER NINE
Solution of Equilibrium Equations in Dynamic Analysis 768
9.1 Introduction 768
9.2 Direct Integration Methods 769
9.2.1 The Central Difference Method, 770
9.2.2 The Houbolt Method, 774
9.2.3 The Wilson 0 Method, 777
9.2.4 The Newmark Method, 780
9.2.5 The Coupling of Different Integration Operators, 782
9.2.6 Exercises, 784
9.3 Mode Superposition 785
9.3.1 Change of Basis to Modal Generalized Displacements, 785
9.3.2 Analysis with Damping Neglected, 789
9.3.3 Analysis with Damping Included, 796
9.3.4 Exercises, 801
9.4 Analysis of Direct Integration Methods 801
9.4.1 Direct Integration Approximation and Load Operators, 803
9.4.2 Stability Analysis, 806
9.4.3 Accuracy Analysis, 810
9.4.4 Some Practical Considerations, 813
9.4.5 Exercises, 822
9.5 Solution of Nonlinear Equations in Dynamic Analysis 824
9.5.1 Explicit Integration, 824
9.5.2 Implicit Integration, 826
9.5.3 Solution Using Mode Superposition, 828
9.5.4 Exercises, 829
9.6 Solution of Nonstructural Problems; Heat Transfer and Fluid Flows 830
9.6.I The aMethod of Time Integration, 830
9.6.2 Exercises, 836
CHAPTER TEN
Preliminaries to the Solution of Eigenproblems 838
10.1 Introduction 838
10.2 Fundamental Facts Used in the Solution of Eigensystems 840
10.2.1 Properties of the Eigenvectors, 841
10.2.2 The Characteristic Polynomials of the Eigenproblem K& = AMd and of Its
Associated Constraint Problems, 845
10.2.3 Shifting, 851
10.2.4 Effect of Zero Mass, 852
10.2.5 Transformation of the Generalized Eigenproblem Kb = AMd 1o a Standard
Form, 854
10.2.6 Exercises, 860
10.3 Approximate Solution Techniques 861
10.3.1 Static Condensation, 861
10.3.2 RayleighRitz Analysis, 868
10.3.3 Component Mode Synthesis, 875
10.3.4 Exercises, 879
10.4 Solution Errors 880 .
10.4.1 Error Bounds, 880
10.4.2 Exercises, 886
CHAPTER ELEVEN o ——
Solution Methods for Eigenproblems 887
11.1 Introduction 887
11.2 Vector Iteration Methods 889
11.2.1 Inverse Iteration, 890
11.2.2 Forward Iteration, 897
11.2.3 Shifting in Vector Iteration, 899
11.2.4 Rayleigh Quotient Iteration, 904
11.2.5 Matrix Deflation and GramSchmidt Orthogonalization, 906
11.2.6 Some Practical Considerations Concerning Vector Iterations, 909
11.2.7 Exercises, 910
11.3 Transformation Methods 911
11.3.1 The Jacobi Method, 912
11.3.2 The Generalized Jacobi Method, 919
11.3.3 The HouseholderQRInverse Iteration Solution, 927
11.3.4 Exercises, 937
11.4 Polynomial Iterations and Sturm Sequence Techniques  938
11.4.1 Explicit Polynomial Iteration, 938
11.4.2 Implicit Polynomial Iteration, 939
11.4.3 Iteration Based on the Sturm Sequence Property, 943
11.4.4 Exercises, 945
11.5 The Lanczos Iteration Method 945
11.5.1 The Lanczos Transformation, 946
11.5.2 Iteration with Lanczos Transformations, 950
11.5.3 Ezxercises, 953
11.6 The Subspace Iteration Method 954
11.6.1 Preliminary Considerations, 955
11.6.2 Subspace lteration, 958
11.6.3 Starting Iteration Vectors, 960
11.6.4 Convergence, 963
11.6.5 Implementation of the Subspace Iteration Method, 964
11.6.6 Exercises, 978
CHAPTER TWELVE
implementation of the Finite Element Method 979
12.1 Introduction 979
12.2 Computer Program Organization for Calculation of System Matrices 980
12.2.1 Nodal Point and Element Information Readin, 981
12.2.2 Calculation of Element Stiffness, Mass, and Equivalent Nodal Loads, 983
12.2.3 Assemblage of Matrices, 983
12.3 Calculation of Element Stresses 987
124 Example Program STAP 988
12.4.1 Data Input to Computer Program STAP, 988
12.4.2 Listing of Program STAP, 995
12.5 Exercises and Projects 1009
References 1013
Index 1029
Finite Element Procedures KlausJiirgen Bathe Professor of Mechanical Engineering Massachusetts Institute of Technology
PRENTICE HALL, Upper Saddle River, New Jersey 07458
Library of Congress CataloginginPublication Data BATHE, KLAUSJURGEN.
Finite element procedures / KlausJiirgen Bathe.
p.
cm
Revision (.)fi Finite element procedures in engineering analysis.
1982. Includes bibliographical references and index. ISBN 0133014584 1. Finite element method.
2. Engineering mathematics.
1. Bathe,
KlausJiirgen. Finite element procedures in engineering analysis. 1I. Title. TA347.F5B36 1996 620'.001'51535—dc20 951231 CIP
Acquisitions editor: William Stenquist Production editor: Merrill Peterson Cover designer: Marianne Frasco Buyer: William Scazzero Editorial assistant: Meg Weist
The authcr and publisher of this book have used their best efforts in preparing
this book. These efforts include the development, research, and testing of the theories and educational computer programs given in the book. The author and publisher make no warranty of any kind, expressed or implied, with regard to any text and the programs contained in this book. The author and publisher
shall not be liable in any event for incidental or consequential damages in connection with, or arising out of, the furnishing, performance, or use of this text and these programs. © 1996 by PrenticeHall, Inc.
Simon & Schuster / A Viacom Company Upper Saddle River, New Jersey 07458 All rights reserved, No part of this book may be reproduced, in any form or by any means, without permission in writing from the publisher. Printed in the United States of America
10
9
8
7
6
5
4
ISBNO01330L4584 PRENTICEHALL INTERNATIONAL (UK) LIMITED, London PRENTICEHALL OF AUSTRALIA PTY. LIMITED, Sydney PRENTICEHALL CANADA INC., Toronto PRENTICEHALL HISPANOAMERICANA, S.A., Mexico PRENTICEHALL OF INDIA PRIVATE LIMITED, New Delhi PRENTICEHALL OF JAPAN, INC., Tokyo
SIMON & SCHUSTER AsIA PTE. LTD., Singapore EDITORA PRENTICEHALL DO BRASIL, LTDA., Rio de Janeiro
To my students
. . . Progress in design of new structures seems to be unlimited. Last
sentence
of article:
“The
Use
of the
Electronic
in Structural Analysis,” by K. J. Bathe Computer (undergraduate student), published in Impact, Journal of the University of Cape Town Engineering Society, pp. 57—
61, 1967.
Contents
Preface
xiii
cHAPTERONE — An Introduction to the Use of Finite Element Procedures
1
1.1
Introduction
1
1.2
Physical Problems, Mathematical Models, and the Finite Element Solution
1.3
Finite Element Analysis as an Integral Part of ComputerAided Design
1.4
A Proposal on How to Study Finite Element Methods
Introduction
22
Introduction to Matrices
23
Vector Spaces
24
Definition of Tensors
2.5
17
17 18
34 40
⋅
The Symmetric Eigenproblem Av = Av
51
2.6
The Rayleigh Quotient and the Minimax Characterization of Eigenvalues 60
2.7
Vector and Matrix Norms
2.8
Exercises
72
11
14
CHAPTERTWO o Vectors, Matrices, and Tensors
2.1
2
66
vi
Contents
CHAPTER THREE Some Basic Concepts of Engineering Analysis and an Introduction to the Finite Element Method 3.1
Introduction
3.2
Solution of DiscreteSystem Mathematical Models 3.2.1 3.2.2 3.2.3 3.2.4 3.2.5
33
3.4
77
77
78
SteadyState Problems, 78 Propagation Problems, 87 FEigenvalue Problems, 90 On the Nature of Solutions, 96 Exercises, 101
Solution of ContinuousSystem Mathematical Models
105
3.3.1 3.3.2 3.3.3 3.3.4
Differential Formulation, 105 Variational Formulations, 110 Weighted Residual Methods; Ritz Method, 116 An Overview: The Differential and Galerkin Formulations, the Principle of Virtual Displacements, and an Introduction to the Finite Element Solution, 124
3.3.5 3.3.6
Finite Difference Differential and Energy Methods, 129 Exercises, 138
Imposition of Constraints
143
3.4.1
An Introduction to Lagrange Multiplier and Penalty Methods, 143
3.4.2
Exercises, 146
CHAPTER FOUR Formulation of the Finite Element Method—Linear Analysis in Solid and Structural Mechanics 4.1
Introduction
42
Formulation of the DisplacementBased Finite Element Method
4.3
4.4
148
4.2.1 4.2.2 4.2.3 4.2.4
General Derivation of Finite Element Equilibrium Equations, 153 Imposition of Displacement Boundary Conditions, 187 Generalized Coordinate Models for Specific Problems, 193 Lumping of Structure Properties and Loads, 212
4.2.5
Exercises, 214
Convergence of Analysis Results
149
225
4.3.1 4.3.2
The Model Problem and a Definition of Convergence, 225 Criteria for Monotonic Convergence, 229
4.3.3 4.3.4 4.3.5 4.3.6 4.3.7
The Monotonically Convergent Finite Element Solution: A Ritz Solution, 234 Properties of the Finite Element Solution, 236 Rate of Convergence, 244 Calculation of Stresses and the Assessment of Error, 254 Exercises, 259
Incompatible and Mixed Finite Element Models 4.4.1 4.4.2 4.4.3 4.4.4
261
Incompatible DisplacementBased Models, 262 Mixed Formulations, 268 Mixed Interpolation— Displacement/Pressure Formulations for Incompressible Analysis, 276 Exercises, 296
148
Contents
45
vii
The InfSup Condition for Analysis of Incompressible Media and Structural Problems 300 4.5.1 4.5.2 4.5.3 4.5.4 4.5.5 4.5.6
The InfSup Condition Derived from Convergence Considerations, 301 The InfSup Condition Derived from the Matrix Equations, 312 The Constant (Physical) Pressure Mode, 315 Spurious Pressure Modes—The Case of Total Incompressibility, 316 Spurious Pressure Modes—The Case of Near Incompressibility, 318 The InfSup Test, 322
4.5.7 4.5.8
An Application to Structural Elements: The Isoparametric Beam Elements, 330 Exercises, 335
00O CHAPTERFIVE Formulationand Calculation of Isoparametric Finite Element Matrices 5.1
Introduction
5.2
Isoparametric Derivation of Bar Element Stiffness Matrix
53
Formulation of Continuum Elements
5.4
55
5.6
338 339
341
5.3.1
Quadrilateral Elements, 342
5.3.2
Triangular Elements, 363
5.3.3 5.3.4 5.3.5 5.3.6
Convergence Considerations, 376 Element Matrices in Global Coordinate System, 386 Displacement/Pressure Based Elements for Incompressible Media, 388 Exercises, 389
Formulation of Structural Elements
397
5.4.1
Beam and Axisymmetric Shell Elements, 399
5.4.2 5.4.3
Plate and General Shell Elements, 420 Exercises, 450
Numerical Integration
455
5.5.1
Interpolation Using a Polynomial, 456
5.5.2
The NewtonCotes Formulas (OneDimensional Integration), 457
5.5.3
The Gauss Formulas (OneDimensional Integration), 461
5.5.4 5.5.5 5.5.6
Integrations in Two and Three Dimensions, 464 Appropriate Order of Numerical Integration, 465 Reduced and Selective Integration, 476
5.5.7
Exercises, 478
Computer Program Implementation of Isoparametric Finite Elements
00O CHAPTER swe Finite Element Nonlinear Analysis in Solid and Structural Mechanics 6.1
Introduction to Nonlinear Analysis
6.2
Formulation of the Continuum Mechanics Incremental Equations of Motion 497 6.2.1 6.2.2
338
485
The Basic Problem, 498 The Deformation Gradient, Strain, and Stress Tensors, 502
480
485
Contents
viii
6.2.3 6.2.4
6.3
DisplacementBased Isoparametric Continuum Finite Elements 6.3.1 6.3.2 6.3.3 6.3.4
6.3.5 6.3.6
6.4
6.6
6.8
568
Beam and Axisymmetric Shell Elements, 568 Plate and General Shell Elements, 575 Exercises, 578
622
Continuum Mechanics Equations, 622 A Solution Approach for Contact Problems: The Constraint Function Method, 626 Exercises, 628
Some Practical Considerations 6.8.1 6.8.2 6.8.3 6.8.4 6.8.5
581
Elastic Material Behavior—Generalization of Hooke’s Law, 583 Rubberlike Material Behavior, 592 Inelastic Material Behavior; Elastoplasticity, Creep, and Viscoplasticity, 595 Large Strain Elastoplasticity, 612 Exercises, 617
Contact Conditions 6.7.1 6.7.2 6.7.3
628
The General Approach to Nonlinear Analysis, 629 Collapse and Buckling Analyses, 630 The Effects of Element Distortions, 636 The Effects of Order of Numerical Integration, 637 Exercises, 640
CHAPTER SEVEN Finite Element Analysis of Heat Transfer, Field Problems, and Incompressible Fluid Flows 7.1
Introduction
7.2
Heat Transfer Analysis 7.2.1 7.2.2 7.2.3 7.2.4
561
Total Lagrangian Formulation, 561 Updated Lagrangian Formulation, 565 Exercises, 566
Use of Constitutive Relations 6.6.1 6.6.2 6.6.3 6.6.4 6.6.5
6.7
Variables, 538 General Matrix Equations of DisplacementBased Continuum Elements, 540 Truss and Cable Elements, 543 TwoDimensional Axisymmetric, Plane Strain, and Plane Stress Elements, 549 ThreeDimensional Solid Elements, 555 Exercises, 557
Structural Elements 6.5.1 6.5.2 6.5.3
538
Linearization of the Principle of Virtual Work with Respect to Finite Element
Displacement/Pressure Formulations for Large Deformations 6.4.1 6.4.2 6.4.3
6.5
Continuum Mechanics Incremental Total and Updated Lagrangian Formulations, Materially NonlinearOnly Analysis, 522 Exercises, 529
642
642
Governing Heat Transfer Equations, 642 Incremental Equations, 646 Finite Element Discretization of Heat Transfer Equations, 651 Exercises, 659
642
ix
Contents
7.3
Analysis of Field Problems 7.3.1 7.3.2 7.3.3 7.3.4 7.3.5
1.4
661
Seepage, 662 Incompressible Inviscid Flow, 663 Torsion, 664 Acoustic Fluid, 666 Exercises, 670
Analysis of Viscous Incompressible Fluid Flows
671
7.4.1 7.4.2
Continuum Mechanics Equations, 675 Finite Element Governing Equations, 677
7.4.3 7.4.4
High Reynolds and High Peclet Number Flows, 682 Exercises, 691
_____ CHAPTER EMGHT — Solution of Equilibrium Equations in Static Analysis
8.1
Introduction
8.2
Direct Solutions Using Algorithms Based on Gauss Elimination 8.2.1 8.2.2 8.2.3 8.2.4
83
8.4
717 726
Solution Errors, 734 Exercises, 741
Iterative Solution Methods 8.3.1 8.3.2 8.3.3
696
Positive Definiteness, Positive Semidefiniteness, and the Sturm Sequence Property,
8.2.6 8.2.7
695
Introduction to Gauss Elimination, 697 The LDLT Solution, 705 Computer Implementation of Gauss Elimination—The Active Column Solution, 708 Cholesky Factorization, Static Condensation, Substructures, and Frontal Solution,
8.2.5
695
745
The GaussSeidel Method, 747 Conjugate Gradient Method with Preconditioning, 749 Exercises, 752
Solution of Nonlinear Equations
754
8.4.1
NewtonRaphson Schemes, 755
8.4.2
The BFGS Method, 759
8.4.3 84.4
LoadDisplacementConstraint Methods, 761 Convergence Criteria, 764
8.4.5
Exercises, 765
CHAPTER NINE Solution of Equilibrium Equations in Dynamic Analysis 9.1
Introduction
9.2
Direct Integration Methods
768
769
9.2.1
The Central Difference Method, 770
9.2.2
The Houbolt Method,
9.2.3
The Wilson 0 Method, 777
774
768
x
Contents
9.2.4 9.2.5 9.2.6
9.3
Mode Superposition 9.3.1 9.3.2 9.3.3 9.3.4
9.4
9.5
9.6
The Newmark Method, 780 The Coupling of Different Integration Operators, 782 Exercises, 784
785
Change of Basis to Modal Generalized Displacements, Analysis with Damping Neglected, 789 Analysis with Damping Included, 796 Exercises, 801
Analysis of Direct Integration Methods
785
801
9.4.1 9.4.2 9.4.3 9.4.4
Direct Integration Approximation and Load Operators, 803 Stability Analysis, 806 Accuracy Analysis, 810 Some Practical Considerations, 813
9.4.5
Exercises, 822
Solution of Nonlinear Equations in Dynamic Analysis 9.5.1 9.5.2 9.5.3
Explicit Integration, 824 Implicit Integration, 826 Solution Using Mode Superposition, 828
9.5.4
Exercises, 829
824
Solution of Nonstructural Problems; Heat Transfer and Fluid Flows 9.6.I The aMethod of Time Integration, 830 9.6.2
830
Exercises, 836
cHAPTERTEN (— Preliminaries to the Solution of Eigenproblems 10.1
Introduction
10.2
Fundamental Facts Used in the Solution of Eigensystems 10.2.1 10.2.2
838
838
840
Properties of the Eigenvectors, 841 The Characteristic Polynomials of the Eigenproblem K&
= AMd
and of Its
Associated Constraint Problems, 845
10.2.3 10.2.4 10.2.5 10.2.6
10.3
10.4
Shifting, 851 Effect of Zero Mass, 852 Transformation of the Generalized Eigenproblem Kb Form, 854 Exercises, 860
Approximate Solution Techniques
861
10.3.1
Static Condensation, 861
10.3.2 10.3.3 10.3.4
RayleighRitz Analysis, 868 Component Mode Synthesis, 875 Exercises, 879
Solution Errors 10.4.1 10.4.2
880
Error Bounds, 880 Exercises, 886
.
= AMd
1o a Standard
xi
Contents
CHAPTER ELEVEN o —— Solution Methods for Eigenproblems 11.1
Introduction
11.2
Vector Iteration Methods 11.2.1 11.2.2 11.2.3 11.2.4 11.2.5 11.2.6 11.2.7
11.3
11.4
11.5
11.6
887
889
Inverse Iteration, 890 Forward Iteration, 897 Shifting in Vector Iteration, 899 Rayleigh Quotient Iteration, 904 Matrix Deflation and GramSchmidt Orthogonalization, 906 Some Practical Considerations Concerning Vector Iterations, 909 Exercises, 910
Transformation Methods 11.3.1 11.3.2 11.3.3 11.3.4
887
911
The Jacobi Method, 912 The Generalized Jacobi Method, 919 The HouseholderQRInverse Iteration Solution, 927 Exercises, 937
Polynomial Iterations and Sturm Sequence Techniques  938 11.4.1 11.4.2 11.4.3
Explicit Polynomial Iteration, 938 Implicit Polynomial Iteration, 939 Iteration Based on the Sturm Sequence Property, 943
11.4.4
Exercises, 945
The Lanczos Iteration Method
945
11.5.1
The Lanczos Transformation, 946
11.5.2 11.5.3
Iteration with Lanczos Transformations, 950 Ezxercises, 953
The Subspace Iteration Method 11.6.1 11.6.2 11.6.3 11.6.4 11.6.5 11.6.6
954
Preliminary Considerations, 955 Subspace lteration, 958 Starting Iteration Vectors, 960 Convergence, 963 Implementation of the Subspace Iteration Method, 964 Exercises, 978
CHAPTER TWELVE implementation of the Finite Element Method 12.1
12.2
Introduction
12.3
979
Computer Program Organization for Calculation of System Matrices 12.2.1 12.2.2 12.2.3
979
980
Nodal Point and Element Information Readin, 981 Calculation of Element Stiffness, Mass, and Equivalent Nodal Loads, 983 Assemblage of Matrices, 983
Calculation of Element Stresses
987
xii
124
Contents
Example Program STAP 12.4.1 12.4.2
12.5
988
Data Input to Computer Program STAP, 988 Listing of Program STAP, 995
Exercises and Projects
1009
References
1013
Index
1029
Preface
Finite element procedures are now an important and frequently indispensable part of engineering analysis and design. Finite element computer programs are now widely used in
practically all branches of engineering for the analysis of structures, solids, and fluids. My objective in writing this book was to provide a text for upperlevel undergraduate and graduate courses on finite element analysis and to provide a book for selfstudy by engineers and scientists. With this objective in mind, I have developed this book from my earlier publication Finite Element Procedures in Engineering Analysis (PrenticeHall, 1982). I have kept the same mode of presentation but have consolidated, updated, and strengthened the earlier writing to the current state of finite element developments. Also, I have added new sections, both to cover some important additional topics for completeness of the presentation and to facilitate (through exercises) the teaching of the material discussed in the book. This text does not present a survey of finite element methods. For such an endeavor, anumber of volumes would be needed. Instead, this book concentrates only on certain finite element procedures, namely, on techniques that I consider very useful in engineering practice and that will probably be employed for many years to come. Also, these methods are introduced in such a way that they can be taught effectively—and in an exciting manner—to students. An important aspect of a finite element procedure is its reliability, so that the method can be used in a confident manner in computeraided design. This book emphasizes this point throughout the presentations and concentrates on finite element procedures that are general and reliable for engineering analysis. Hence, this book is clearly biased in that it presents only certain finite element procedures and in that it presents these procedures in a certain manner. In this regard, the book reflects my philosophy toward the teaching and the use of finite element methods.
xiv
Preface
While the basic topics of this book focus on mathematical methods, an exciting and thorough understanding of finite element procedures for engineering applications is achieved only if sufficient attention is given to both the physical and mathematical characteristics of the procedures. The combined physical and mathematical understanding greatly enriches our confident use and further development of finite element methods and is therefore emphasized in this text. These thoughts also indicate that a collaborauon between engineers and mathematicians to deepen our understanding of finite element methods and to further advance in the fields of research can be of great benefit. Indeed, I am thankful to the mathematician Franco Brezzi for our research collaboration in this spirit, and for his valuable suggestions regarding this book. I consider it one of the greatest achievements for an educator to write a valuable book. In these times, all fields of engineering are rapidly changing, and new books for students are needed in practically all areas of engineering. I am therefore grateful that the Mechanical Engineering Department of M.L.T. has provided me with an excellent environment in which. to pursue my interests in teaching, research, and scholarly writing. While it required an immense effort on my part to write this book, I wanted to accomplish this task as a commitment to my past and future students, to any educators and researchers who might have an interest in the work, and, of course, to improve upon my teaching at M.I.T. I have been truly fortunate to work with many outstanding students at M.I.T., for which I am very thankful. It has been a great privilege to be their teacher and work with them. Of much value has also been that I have been intimately involved, at my company ADINA R & D, Inc., in the development of finite element methods for industry. This involvement has been very beneficial in my teaching and research, and in my writing of this book. A text of significant depth and breadth on a subject that came to life only a few decades ago and that has experienced tremendous advances, can be written only by an author who has had the benefit of interacting with many people in the field. I would like to thank all my students and friends who contributed—and will continue to contribute—to my knowledge and understanding of finite element methods. My interaction with them has given me great joy and satisfaction. I also would like to thank my secretary, Kristan Raymond, for her special efforts in typing the manuscript of this text. Finally, truly unbounded thanks are due to my wife, Zorka, and children, Ingrid and Mark, who, with their love and their understanding of my efforts, supported me in writing this book.
K. J. Bathe
Il CHAPTER ONE I
An Introduction to the Use of Finite Element Procedures
1.1 INTRODUCTION Finite element procedures are at present very widely used in engineering analysis, and we can expect this use to increase significantly in the years to come. The procedures are employed extensively in the analysis of solids and structures and of heat transfer and fluids, and indeed, finite element methods are useful in virtually every field of engineering analysis. The development of finite element methods for the solution of practical engineering problems began with the advent of the digital computer. That is, the essence of a finite element solution of an engineering problem is that a set of governing algebraic equations is established and solved, and it was only through the use of the digital computer that this process could be rendered effective and given general applicability. These two properties— effectiveness and general applicability in engineering analysis—are inherent in the theory used and have been developed to a high degree for practical computations, so that finite element methods have found wide appeal in engineering practice. ⋟ As is often the case with original developments, it is rather difficult to quote an exact “date of invention,” but the roots of the finite element method can be traced back to three
separate research groups: applied mathematicians—see R. Courant [A]; physicists—see
J. L. Synge [A]; and engineers—see J. H. Argyris and S. Kelsey [A]. Although in principle published already, the finite element method obtained its real impetus from the developments of engineers. The original contributions appeared in the papers by J. H. Argyris and S. Kelsey [A]; M. J. Turner, R. W. Clough, H. C. Martin, and L. J. Topp [A]; and R. W.
Clough [A]. The name “finite element” was coined in the paper by R. W. Clough [A]. Important early contributions were those of J. H. Argyris [A] and O. C. Zienkiewicz and Y. K. Cheung [A]. Since the early 1960s, a large amount of research has been devoted to the technique, and a very large number of publications on the finite element method is 1
2
An Introduction to the Use of Finite Element Procedures
Chap. 1
available (see, for example, the compilation of references by A. K. Noor [A] and the Finite Element Handbook edited by H. Kardestuncer and D. H. Norrie [A)). The finite element method in engineering was initially developed on a physical basis for the analysis of problems in structural mechanics. However, it was soon recognized that the technique could be applied equally well to the solution of many other classes of problems. The objective of this book is to present finite element procedures comprehensively and in a broad context for solids and structures, field problems (specifically heat transfer), and fluid flows. To introduce the topics of this book we consider three important items in the following sections of this chapter. First, we discuss the important point that in any analysis we always select a mathematical model of a physical problem, and then we solve that model. The finite element method is employed to solve very complex mathematical models, but it is important to realize that the finite element solution can never give more information than that
contained in the mathematical model. Then we discuss the importance of finite element analysis in the complete process of computeraided design (CAD). This is where finite element analysis procedures have their greatest utility and where an engineer is most likely to encounter the use of finite element methods. In the last section of this chapter we address the question of how to study finite element methods. Since a voluminous amount of information has been published on these techniques, it can be rather difficult for an engineer to identify and concentrate on the most important principles and procedures. Our aim in this section is to give the reader some guidance in studying finite element analysis procedures and of course also in studying the
various topics discussed in this book.
1.2 PHYSICAL PROBLEMS, MATHEMATICAL MODELS, AND THE FINITE ELEMENT SOLUTION The finite element method is used to solve physical problems in engineering analysis and design. Figure 1.1 summarizes the process of finite element analysis. The physical problem typically involves an actual structure or structural component subjected to certain loads. The idealization of the physical problem to a mathematical model requires certain assumptions that together lead to differential equations governing the mathematical model (see Chapter 3). The finite element analysis solves this mathematical model. Since the finite element solution technique is a numerical procedure, it is necessary to assess the solution accuracy. If the accuracy criteria are not met, the numerical (i.e., finite element) solution has to be repeated with refined solution parameters (such as finer meshes) until a sufficient accuracy is reached. It is clear that the finite element solution will solve only the selected mathematical model and that all assumptions in this model will be reflected in the predicted response. We cannot expect any more information in the prediction of physical phenomena than the information contained in the mathematical model. Hence the choice of an appropriate mathematical model is crucial and completely determines the insight into the actual physical problem that we can obtain by the analysis,
Sec. 1.2
Physical Problems, Mathematical Models, the Finite Element Solution
3
Change of physical problem
Physical problem
Mathematical model Governed by differential equations Assumptions on * Geometry
Improve
+ Kinematics Meterial law « Loading * Boundery conditions * Etc.
mathematical model
Finite element solution Choice of ¢ Finite elements
» Mesh density » Solution parameters Representation of
Finite element
* Loading
of
mathematicel model
e
solution
Boundery conditions
* Etc.
Refine mesh,
solution parameters, etc.
Assessment of eccurecy of finite element solution of mathematical model
Interpretation of results
analysis
Design improvements Structural optimlzation Figure 1.1
The process of finite element analysis
Let us emphasize that, by our analysis, we can of course only obtain insight into the physical problem considered: we cannot predict the response of the physical problem exactly because it is impossible to reproduce even in the most refined mathematical model all the information that is present in nature and therefore contained in the physical problem. Once a mathematical model has been solved accurately and the results have been interpreted; we may well decide to consider next a refined mathematical model in order to increase our insight into the response of the physical problem. Furthermore, a change in the physical problem may be necessary, and this in turn will also lead to additional mathematical models and finite element solutions (see Fig. 1.1).
The key step in engineering analysis is therefore choosing appropriate mathematical models. These models will clearly be selected depending on what phenomena are to be
4
An Introduction to the Use of Finite Element Procedures
Chap. 1
predicted, and it is most important to select mathematical models that are reliable and effective in predicting the quantities sought. To define the reliability and effectiveness of a chosen model we think of a verycomprehensive mathematical model of the physical problem and measure the response of our chosen model against the response of the comprehensive model. In general, the verycomprehensive mathematical model is a fully threedimensional description that also includes nonlinear effects. ⊾
Effectiveness of a mathematical model The most effective mathematical model for the analysis is surely that one which yields the required response to a sufficient accuracy and at least cost. Reliability of a mathematical model The chosen mathematical model is reliable if the required response is known to be predicted within a selected level of accuracy measured on the response of the verycomprehensive mathematical model. Hence to assess the results obtained by the solution of a chosen mathematical model, it may be necessary to also solve higherorder mathematical models, and we may well think of (but of course not necessarily solve) a sequence of mathematical models that include increasingly more complex effects. For example, a beam structure (using engineering terminology) may first be analyzed using Bernoulli beam theory, then Timoshenko beam theory, then twodimensional plane stress theory, and finally using a fully threedimensional continuum model, and in each case nonlinear effects may be included. Such a sequence of models is referred to as a hierarchy of models (see K. J. Bathe, N. S. Lee, and M. L. Bucalem
[A]). Clearly, with these hierarchical models the analysis will include ever more complex response effects but will also lead to increasingly more costly solutions. As is well known, a fully threedimensional analysis is about an order of magnitude more expensive (in computer resources and engineering time used) than a twodimensional solution. Let us consider a simple example to illustrate these ideas. ⊽ Figure 1.2(a) shows a bracket used to support a vertical load. For the analysis, we need to choose a mathematical model. This choice must clearly depend on what phenomena are to be predicted and on the geometry, material properties, loading, and support conditions of the bracket. We have indicated in Fig. 1.2(a) that the bracket is fastened to a very thick steel column. The description “very thick” is of course relative to the thickness ¢ and the height h of the bracket. We translate this statement into the assumption that the bracket is fastened to a (practically) rigid column. Hence we can focus our attention on the bracket by applying a “rigid column boundary condition” to it. (Of course, at a later time, an analysis of the column may be required, and then the loads carried by the two bolts, as a consequence of the load W, will need to be applied to the column.) We also assume that the load W is applied very slowly. The condition of time “very slowly” is relative to the largest natural period of the bracket; that is, the time span over
which the load W is increased from zero to its full value is much longer than the fundamental period of the bracket. We translate this statement into requiring a static analysis (and not a dynamic analysis).
Sec. 1.2
Physical Problems, Mathematicel Models, the Finite Element Solution
5
With these preliminary considerations we can now establish"an appropriate mathematical model for the analysis of the bracket—depending on what phenomena are to be predicted. Let us assume, in the first instance, that only the total bending moment at section AA in the bracket and the deflection at the load application are sought. To predict these quantities, we consider a beam mathematical model including shear deformations [see Fig. 1.2(b)] and obtain M=
WL 1.1
= 27,500 N cm
an
1WIL+rm)®  WL+ )
8 latioaaw = 3
EI
+ T 3AG 6
(1.2)
= (0.053 cm
Uniform Two bolts
thickness t
W= 1000 N L =275cm = 0.5cm
E =2x 107 N/cm? v
=03
h
=6.0cm
t
=04cm
Pin
(a) Physical problem of steel bracket :A ~
‘< I
ry=05cm
I
i I
! I
% I I I I
L+rN=28cm
(b) Beam model
Figure 1.2
Bracket to be analyzed and two mathematical models
An Introduction to the Use of Finite Element Procedures
6
Chap. 1
Areas with imposed zero displecements u, v
y,v
Load applied atpoint B
Zw
xu Equilibrium equations (see Example 4.2) OTex 
a
OTxy _
x
o
dx
ay
. . in domain of bracket
0T  8Ty _
7 = 0, T, = 0 on surfaces except at point B and at imposed zero displacements Stressstrain relation (see Table 4.3): Tax Ty  =
E 2
1 vy
v
0
€xx
1
0
€y
(1u2]{7v
00
Tay
E = Young’s modulus, » = Poisson’s ratio Straindisplacement relations (see Section 4.2): €xx
’
€yy
e
»
Yy
_w ay
ax
(c) Plene stress model
Figure 1.2
(continued)
where L and rv are given in Fig. 1.2(a), E is the Young’s modulus of the steel used, G is the shear modulus, 7 is the moment of inertia of the bracket arm (I = $5h3%1), A is the crosssectional area (A = ht), and the factor 2 is a shear correction factor (see Section 5.4.1). Of course, the relations in (1.1) and (1.2) assume linear elastic infinitesimal displacement conditions, and hence the load must not be so large as to cause yielding of the material and/or large displacements. Let us now ask whether the mathematical model used in Fig. 1.2(b) was reliable and effective. To answer this question, strictly, we should consider a verycomprehensive mathematical model, which in this case would be a fully threedimensional representation of the
Sec. 1.2
Physical Problems, Mathematical Models, the Finite Element Solution
7
full bracket. This model should include the two bolts fastening the bracket to the (assumed rigid) column as well as the pin through which the load W is applied. The threedimensional solution of this model using the appropriate geometry and material data would give the numbers against which we would compare the answers given in (1.1) and (1.2). Note that this threedimensional mathematical model contains contact conditions (the contact is between the bolts, the bracket, and the column, and between the pin carrying the load and the bracket) and stress concentrations in the fillets and at the holes. Also, if the stresses are high, nonlinear material conditions should be incluped in the model. Of course, an
analytical solution of this mathematical model is not available, and all we can obtain is a numerical solution. We describe in this book how such solutions can be calculated using finite element procedures, but we may note here already that the solution would be relatively expensive in terms of computer resources and engineering time used. Since the threedimensional comprehensive mathematical model is very likely too comprehensive a model (for the analysis questions we have posed), we instead may consider alinear elastic twodimensional plane stress model as shown in Fig. 1.2(c). This mathematical model represents the geometry of the bracket more accurately than the beam model and assumes a twodimensional stress situation in the bracket (see Section 4.2). The bending moment at section AA and deflection under the load calculated with this model can be expected to be quite close to those calculated with the verycomprehensive threedimensional model, and certainly this twodimensional model represents a higherorder model against which we can measure the adequacy of the results given in (1.1) and (1.2). Of course, an analytical solution of the model is not available, and a numerical solution must
be sought. Figures 1.3(a) to (¢) show the geometry and the finite element discretization used in
the solution of the plane stress mathematical model and some stress and displacement results obtained with this discretization. Let us note the various assumptions of this mathematical model when compared to the more comprehensive threedimensional model discussed earlier. Since a plane stress condition is assumed, the only nonzero stresses are 7xy, Tyy, and 7xy. Hence we assume that the stresses 7., 7y, and 7, are zero. Also, the actual bolt fastening and contact conditions between the steel column and the bracket are not included
(a) Geometry of bracket as obtained from a CAD program
Figure 1.3
Plane stress analysis of bracket in Fig. 1.2. AutoCAD was used to create the
geometry, and ADINA was used for the finite element analysis.
{b)
Mesh of ninenode elements cretization
used in finite element dis
(c)
Deflected shape. Deflections are drawn with a magnification factor of 100 together with the original configuration SIGMAP1
SMOOTHED
TIME 1.000
SIGMAPY TIME 1 000
—
17500. 12500.
7500. 2500.
5
17500,
E
125Q0. 7500.
2500.
{d)
Maximum principal stress near notch. Unsmoothed stress results are shown. The small breaks in the bands indicate that a reasonably accurate solution of the mathematical model has been obtained (see Section 4.3.6)
Figure 1.3
{e)
Maximum principal stress near notch. Smoothed stress results. (The averages of nodal point stresses are taken and interpolated over the elements.)
(continued)
Sec. 1.2
Physical Problems, Mathematicsl
Models, the Finite Element Solution
9
in the model, and the pin carrying the load into the bracket is not modeled. However, since our objective is only to predict the bending moment at section AA and the deflection at point B, these assumptions are deemed reasonable and of relatively little influence. Let us assume that the results obtained in the finite element solution of the mathematical model are sufficiently accurate that we can refer to the solution given in Fig. 1.3. as the solution of the plane stress mathematical model. Figure 1.3(c) shows the calculated deformed configuration. The deflection at the point of load application B as predicted in the plane stress solution is
8atiosaw = 0.064 cm
(1.3)
Also, the total bending moment predicted at section AA is
M=o = 27,500 N cm
(1.4)
Whereas the same magnitude of bending moment at section AA is predicted by the beam model and the plane stress model,' the deflection of the beam model is considerably less than that predicted by the plane stress model [because of the assumption that the beam in Fig. 1.2(b) is supported rigidly at its left end, which neglects any deformation between the beam end and the boits]. Considering these results, we can say that the beam mathematical model in Fig. 1.2(b) is reliable if the required bending moment is to be predicted within 1 percent and the deflection is to be predicted only within 20 precent accuracy. The beam model is of course also effective because the calculations are performed with very little effort. On the other hand, if we next ask for the maximum stress in the bracket, then the simple mathematical beam model in Fig. 1.2(b) will not yield a sufficiently accurate answer. Specifically, the beam model totally neglects the stress increase due to the fillets.> Hence a plane stress solution including the fillets is necessary. The important points to note here are the following.
1. The selection of the mathematical model must depend on the response to be predicted (i.e., on the questions asked of nature). 2. The most effective mathematical model is that one which delivers the answers to the
questions in a reliable manner (i.e., within an acceptable error) with the least amount of effort. 3. A finite element solution can solve accurately only the chosen mathematical model (e.g., the beam model or the plane stress model in Fig. 1.2) and cannot predict any more information than that contained in the model. 4. The notion of reliability of the mathematical model hinges upon an accuracy assessment of the results obtained with the chosen mathematical model (in response to the
questions asked) against the results obtained with the verycomprehensive mathematical model. However, in practice the verycomprehensive mathematical model is ! The bending moment at section AA in the plane stress model is calculated here froin the finite element nodal point forces, and for this statically determinate analysis problem the internal resisting moment must be equal to the externally applied moment (see Example 4.9).
20f course, the effect of the fillets could be estimated by the use of stress concentration factors that have been established from plane stress solutions.
10
An Introduction to the Use of Finite Element Procedures
Chap. 1
usually not solved, and instead engineering experience is used, or a more refined mathematical model is solved, to judge whether the mathematical model used was adequate (i.e., reliable) for the response to be predicted.
Finally, there is one further important general point. The chosen mathematical model may contain extremely high stresses because of sharp corners, concentrated loads, or other effects. These high stresses may be due solely to the simplifications used in the model when compared with the verycomprehensive mathematical model (or with nature). For example, the concentrated load in the plane stress model in Fig. 1.2(c) is an idealization of a pressure load over a small area. (This pressure would in nature be transmitted by the pin carrying the load into the bracket.) The exact solution of the mathematical model in Fig. 1.2(c) gives an infinite stress at the point of load application, and we must therefore expect a very large stress at point B as the finite element mesh is refined. Of course, this very large stress is an artifice of the chosen model, and the concentrated load should be replaced by a pressure load over a small area when a very fine discretization is used (see further discussion). Furthermore, if the model then still predicts a very high stress, a nonlinear mathematical model would be appropriate. Note that the concentrated load in the beam model in Fig. 1.2(b) does not introduce any solution difficulties. Also, the rightangled sharp corners at the support of the beam model, of course, do not introduce any solution difficulties, whereas such corners in a plane stress model would introduce infinite stresses. Hence, for the plane stress model, the corners
have to be rounded to more accurately represent the geometry of the actual physical bracket. We thus realize that the solution of a mathematical model may result in artificial difficulties that are easily removed by an appropriate change in the mathematical model to more closely represent the actual physical situation. Furthermore, the choice of a more encompassing mathematical model may result, in such cases, in a decrease in the required solution effort. While these observations are of a general nature, let us consider once again, specifically, the use of concentrated loads. This idealization of load application is extensively used in engineering practice. We now realize that in many mathematical models (and therefore also in the finite element solutions of these models), such loads create stresses of infinite value. Hence, we may ask under what conditions in engineering practice solution difficulties may arise. We find that in practice solution difficulties usually arise only when the finite element discretization is very fine, and for this reason the matter of infinite stresses under concentrated loads is frequently ignored. As an example, Fig. 1.4 gives finite element results obtained in the analysis of a cantilever, modeled as a plane stress problem. The cantilever is subjected to a concentrated tip load. In practice, the 6 X 1 mesh is usually considered sufficiently fine, and clearly, a much finer discretization would have to be used
to accurately show the effects of the stress singularities at the point of load application and at the support. As already pointed out, if such a solution is pursued, it is necessary to change the mathematical model to more accurately represent the actual physical situation of the structure. This change in the mathematical model may be important in selfadaptive finite element analyses because in such analyses new meshes are generated automatically and artificial stress singularities cause—artificially—extremely fine discretizations. We refer to these considerations in Section 4.3.4 when we state the general elasticity problem considered for finite element solution.
Finite Element Analysis as an Integral Part of ComputerAided Design
Sec. 1.3
1
W=0.1
4
7
?
E = 200,000 v=0.30 Thickness = 0.1
1.0
o
⊤ (a) Geometry, boundary conditions, and material data. Bernoulli beam theory results: § = 0.16, Tmax = 120
(b) Typical finite element discretization, 6 x 1 mesh of 9node elements; results are: 8 = 0.16, Tmax = 116
Figure 1.4
Analysis of a cantilever as a plane stress problem
In summary, we should keep firmly in mind that the crucial step in any finite element analysis is always choosing an appropriate mathematical model since a finite element solution solves only this model. Furthermore, the mathematical model must depend on the analysis questions asked and should be reliable and effective (as defined earlier). In the process of analysis, the engineer has to judge whether the chosen mathematical model has been solved to a sufficient accuracy and whether the chosen mathematical model was appropriate (i.e., reliable) for the questions asked. Choosing the mathematical model, solving the model by appropriate finite element procedures, and judging the results are the fundamental ingredients of an engineering analysis using finite element methods.
1.3 FINITE ELEMENT ANALYSIS AS AN INTEGRAL PART OF COMPUTERAIDED DESIGN Although a most exciting field of activity, engineering analysis is clearly only a support activity in the larger field of engineering design. The analysis process helps to identify good new designs and can be used to improve a design with respect to performance and cost. In the early use of finite element methods, only specific structures were analyzed, mainly in the aerospace and civil engineering industries. However, once the full potential of finite element methods was realized and the use of computers increased in engineering design environments, emphasis in research and development was placed upon making the use of finite element methods an integral part of the design process in mechanical, civil, and aeronautical engineering.
An Introduction to the Use of Finite Element Procedures
12
Chap. 1
Geomatry generation
Finite elament analysis
Kinematic
Interactive use of software
analysis
Automatic drafting
Figure 1.5
Management of process
The field of CAD/CAM viewed schematically
Figure 1.5 gives an overview of the steps in a typical computeraided design process. Finite element analysis is only a small part of the complete process, but it is an important
part.
We note that the first step in Figure 1.5 is the creation of a geometric representation of the design part. Many different computer programs can be employed (e.g., a typical and popular program is AutoCAD). In this step, the material properties, the applied loading and boundary conditions on the geometry also need to be defined. Given this information, a finite element analysis may proceed. Since the geometry and other data of the actual physical part may be quite complex, it is usually necessary to simplify the geometry and loading in order to reach a tractable mathematical model. Of course, the mathematical model should be reliable and effective for the analysis questions posed, as discussed in the preceding section. The finite element analysis solves the chosen mathematical model, which may be changed and evolve depending on the purpose of the analysis (see Fig. 1.1). Considering this process—which generally is and should be performed by engineering designers and not only specialists in analysis—we recognize that the finite element methods must be very reliable and robust. By reliability of the finite element methods we now’ mean that in the solution of a wellposed mathematical model, the finite element procedures should always for a reasonable finite element mesh give a reasonable solution, 3 Note that this meaning of “reliability of finite element methods” is different from that of “reliability of a mathematical model” defined in the previous section.
Sec. 1.3
Finite Element Analysis as an Integral Part of ComputerAided Design
13
and if the mesh is reasonably fine, an accurate solution should always be obtained. By robustness of the finite element methods we mean that the performance of the finite element procedures should not be unduly sensitive to the material data, the boundary conditions, and the loading conditions used. Therefore, finite element procedures that are not robust will also not be reliable. For example, assume that in the plane stress solution of the mathematical model in Fig. 1.2(c), any reasonable finite element discretization using a certain element type is employed. Then the solution obtained from any such analysis should not be hugely in error, that is, an order of magnitude larger (or smaller) than the exact solution. Using an unreliable finite element for the discretization would typically lead to good solutions for some mesh topologies, whereas with other mesh topologies it would lead to bad solutions. Elements based on reduced integration with spurious zero energy modes can show this unreliable behavior (see Section 5.5.6). Similarly, assume that a certain finite element discretization of a mathematical model gives accurate results for one set of material parameters and that a small change in the parameters corresponds to a small change in the exact solution of the mathematical model. Then the same finite element discretization should also give accurate results for the mathematical model with the small change in material parameters and not yield results that are very much in error. These considerations regarding effective finite element discretizations are very important and are discussed in the presentation of finite element discretizations and their stability and convergence properties (see Chapters 4 to 7). For use in engineering design, it is of utmost importance that the finite element methods be reliable, robust, and of course
efficient. Reliability and robustness are important because a designer has relatively little time for the process of analysis and must be able to obtain an accurate solution of the chosen mathematical model quickly and without “trial and error.” The use of unreliable finite element methods is simply unacceptable in engineering practice. An important ingredient of a finite element analysis is the calculation of error estimates, that is, estimates of how closely the finite element solution approximates the exact solution of the mathematical model (see Section 4.3.6). These estimates indicate whether a specific finite element discretization has indeed yielded an accurate response prediction, and a designer can then rationally decide whether the given results should be used. In the case that unacceptable results have been obtained using unreliable finite element methods, the difficulty is of course how to obtain accurate results. Finally, we venture to comment on the future of finite element methods in computeraided design. Surely, many engineering designers do not have time to study finite element methods in depth or finite element procedures in general. Their sole objective is to use these techniques to enhance the design product. Hence the integrated use of finite element methods in CAD in the future will probably involve to an increasingly smaller degree the scrutiny of finite element meshes during the analysis process. Instead, we expect that in linear elastic static analysis the finite element solution process will be automatized such that, given a mathematical model to be solved, the finite element meshes will be automatically created,
the solution will be calculated with and the desired solution accuracy, refined (without the analyst or the required solution accuracy has been
error estimates, and depending on the estimated errors the finite element discretization will be automatically designer ever seeing a finite element mesh) until the attained. In this automatic analysis process—in which
14
An Introduction to the Use of Finite Element Procedures
Chap. 1
of course the engineer must still choose the appropriate mathematical model for analysis— the engineer can concentrate on the design aspects while using analysis tools with great efficiency and benefit. While this design and analysis environment will be commonly available for linear elastic static analysis, the dynamic or nonlinear analysis of structures and fluids will probably still require from the engineer a greater degree of involvement and expertise in finite element methods for many years to come. With these remarks we do not wish to suggest overconfidence but to express a realistic outlook with respect to the valuable and exciting future use of finite element procedures. For some remarks on overconfidence in regard to finite element methods (which are still pertinent after almost two decades), see the article “A Commentary on Computational Mechanics” by J. T. Oden and K. J. Bathe [A].
1.4 A PROPOSAL ON HOW TO STUDY FINITE ELEMENT METHODS With a voluminous number of publications available on the use and development of finite element procedures, it may be rather difficult for a student or teacher to identify an effective plan of study. Of course, such a plan must depend on the objectives of the study and the time available. In very general terms, there are two different objectives: 1. To learn the proper use of finite element methods for the solution of complex problems, and 2. To understand finite element methods in depth so as to be able to pursue research on finite element methods. This book has been written to serve students with either objective, recognizing that the population of students with objective 1 is the much larger one and that frequently a student may first study finite element methods with objective 1 and then develop an increasing interest in the methods and also pursue objective 2. While a student with objective 2 will need to delve still much deeper into the subject matter than we do in this book, it is hoped that this book will provide a strong basis for such further study and that the “philosophy” of this text (see Section
1.3) with respect to the use of and research into finite element
methods will be valued. Since this book has not been written to provide a broad survey of finite element methods—indeed, for a survey, a much more comprehensive volume would be necessary— it is clearly biased toward a particular way of introducing finite element procedures and toward only certain finite element procedures. The finite element methods that we concentrate on in this text are deemed to be effective techniques that can be used (and indeed are abundantly used) in engineering practice. The methods are also considered basic and important and will probably be employed for a long time to come. The issue of what methods should be used, and can be best used in engineering practice, is important in that only reliable techniques should be employed (see Section 1.3), and this book concentrates on the discussion of only such methods.
Sec. 1.4
A Proposal on How to Study Finite Element Methods
15
In introducing finite element procedures we endeavor to explain the basic concepts and equations with emphasis on physical explanations. Also, short examples are given to demonstrate the basic concepts used, and exercises are provided for hand calculations and the use of a finite element computer program. The preceding thoughts lead to our proposal for the study of finite element procedures. If objective 1 is being pursued, the reader will be mostly interested in the basic formulations, the properties of the finite elements and solution algorithms, and issues of convergence and efficiency. An important point to keep in mind is that numerical finite element procedures are used to solve a mathematical model (with some reasonably small solution errors) (see Section 1.2). Hence, it is important for the user of a finite element computer program to always be able to judge the quality of the finite element results obtained in a solution. We demonstrate in this book how such judging is properly performed. However, if objective 2 is being pursued, much more depth in the formulations and numerical algorithms must be acquired. This text provides a broad basis for such study (but of course does not give all details of derivations and implementations). In either case, we believe that for a study of finite element methods, it is effective to use a finite element computer program while learning about the theory and numerical procedures. Then, at the same time as the theoretical concepts are being studied on paper, it will be possible to test and demonstrate these concepts on a computer. Based on this book, various courses can be taught depending on the aim of the instruction. A first course would be “Introduction to Finite Element Analysis,” which could be based on Sections 1.1t0 1.4,2.1t02.3,3.1t0 3.4,4.1t04.3,5.1t05.3, and 8.1 t0 8.2.2.
A more advanced course would be “Finite Element Procedures,” based on Sections 1.1. to 1.4, 3.1, 3.3, 4.1 to 4.4, 5.1 t0 5.6, and 8.1 to 8.3. A course on finite element methods for dynamic analysis would be “Computer Methods in Dynamics™ and could be based on Sections 1.1 to 1.4, 2.1 to 2.7, 3.1, 3.3, 4.1, 4.2,
51t05.3,8.110822,9.1t094,10.1, 10.2, 11.1, 11.2.1, 11.3.1, and 11.6. A course emphasizing the continuum mechanics principles for linear and nonlinear finite element analysis would be “Theory and Practice of Continuum Mechanics,” which could be based on Sections 1.1t0 1.4,3.1,3.3,4.1,4.2.1,4.2.2,5.1,5.2,5.3.1,5.3.3,5.3.5, 6.1,6.2,6.3.1,6.3.2,6.4.1,6.6,7.1,and 7.4, A course on the analysis of field problems and fluid flows would be “Finite Element Analysis of Heat Transfer, Field Problems, and Fluid Flows,” based on Sections 1.1 to 1.4, 3.1,3.3,7.1t0 7.4, 5.3, 5.5, and 4.5.1 t0 4.5.6. Note that the presentation in this course
would first provide the finite element formulations (Sections 7.1 to 7.4) and then numerical procedures and mathematical results. A course on nonlinear finite element analysis of solids and structures would be “Nonlinear Analysis of Solids and Structures” and could be based on Sections 1.1 to 1.4, 6.1 t0 6.8, 8.1, and 8.4.
A course
on the numerical
methods
used
in finite element
analysis
would
be
“Numerical Methods in Finite Element Analysis,” which could be based on Sections 1.1 to 14,2.1102.7,4.1,4.2.1,5.1,5.3,5.5,8.1t0 8.4,9.1 0 9.6, 10.1 t0 10.3,and 11.1 to 11.6. And, finally, a course on more theoretical aspects of finite element analysis could be offered, entitled “Theoretical Concepts of Finite Element Analysis,” based on Sections 1.1
to 1.4, 4.1 to 4.5, and 5.1 to 5.5.
16
An Introduction to the Use of Finite Element Procedures
Chap. 1
In addition, in most of these courses, the material in Chapter 12 would be useful. In general, such courses are best taught with homework that includes problem solutions, as suggested in the exercise sections of this book, as well as the completion of a project in which a finite element computer program is used. Various projects are proposed in the exercise sections of this book. One type of project requires the student to use the program STAP and to modify this program for certain new capabilities as suggested in Section 12.5. Such a project would be of value to students interested in learning about the actual detailed implementation of finite element methods. The other type of project is based on the use of a standard (commercial) finite element program—as it might be employed in engineering practice—to analyze an interesting physical problem. A valuable way to then proceed is to first solve a mathematical model for which an analytical exact solution is known and then introduce complicating features in the mathematical model that require a numerical solution. In this way, studies using different finite element discretizations and solution procedures in which the results are evaluated against a known exact solution can first be performed, which establishes confidence in the use of the methods. The required expertise and confidence then become valuable assets in the finite element solution of more complicated mathematical models.It is of particular importance during the study of finite element procedures that good judgment be developed and exercised concerning the quality of any finite element solution and of any mathematical model considered. Of course, exercising such judgments requires that the analyst be sufficiently familiar with possible mathematical models and available finite element procedures, and this aspect stimulates learning about mathematical models and methods of solution. The value of using a finite element computer code—the program STAP or a commercial finite element program as previously mentioned—is clearly in stimulating the process of learning and in reinforcing the understanding of what has been learned in class and from homework problem solutions. Using a powerful analysis program, in particular, also adds to the excitement of actually solving complicated analysis problems that heretofore could not be tackled.
Il CHAPTER TWO I
Vectors, Matrices, and Tensors
2.1 INTRODUCTION The use of vectors, matrices, and tensors is of fundamental importance in engineering analysis because it is only with the use of these quantities that the complete solution process can be expressed in a compact and elegant manner. The objective of this chapter is to present the fundamentals of matrices and tensors, with emphasis on those aspects that are important in finite element analysis. From a simplistic point of view, matrices can simply be taken as ordered arrays of numbers that are subjected to specific rules of addition, multiplication, and so on. It is of course important to be thoroughly familiar with these rules, and we review them in this chapter. However, by far more interesting aspects of matrices and matrix algebra are recognized when we study how the elements of matrices are derived in the analysis of a physical problem and why the rules of matrix algebra are actually applicable. In this context, the use of tensors and their matrix representations are important and provide a most interesting subject of study. Of course, only a rather limited discussion of matrices and tensors is given here, but we hope that the focused practical treatment will provide a strong basis for understanding the finite element formulations given later.
17
Vectors, Matrices, and Tensors
18
Chap. 2
2.2 INTRODUCTION TO MATRICES The effectiveness of using matrices in practical calculations is readily realized by considering the solution of a set of linear simultaneous equations such as S5x; — 4x; + —4x,
x3
= (
+ 6x2 —4dxs +
Xy —4x;
x4 =1
2.1
+ 6x3 —4x4 =0
X2 —4x3 + Sx4 =0
where the unknowns are x;, x2, x3, and x,. Using matrix notation, this set of equations is written as 5
—4
1
’_4
6
—4
1 0
4
6 1
4
01
[«
1
X2
0 1 0 0
_
—4f(x 5 [ xe
(2.2)
where it is noted that, rather logically, the coefficients of the unknowns (5, —4, 1, etc.) are grouped together in one array; the lefthandside unknowns (x:, x2, x3, and x4) and the
righthandside known quantities are each grouped together in additional arrays. Although written differently, the relation (2.2) still reads the same way as (2.1). However, using
matrix symbols to represent the arrays in (2.2), we can now write the set of simultaneous equations as
Ax=b
(2.3)
where A is the matrix of the coefficients in the set of linear equations, x is the matrix of unknowns, and b is the matrix of known quantities; i.e.,
5 A; and that p(v) approximates A; more closely than v approximates v;.
Sec. 2.6
The Rayleigh Quotient and the Minimax Characterization of Eigenvalues
63
Having introduced the Rayleigh quotient, we can now proceed to a very important principle, the minimax characterization of eigenvalues. We know from Rayleigh’s principle that p(v) = A (2.125) where v is any vector. In other words, if we consider the problem of varying v, we will always have p(v) = A,, and the minimum will be reached when v = v;, in which case
p(v1) = A,. Suppose that we now impose a restriction on v, namely that v be orthogonal to a specific vector w, and that we consider the problem of minimizing p(v) subject to this restriction. After calculating the minimum of p(v) with the condition v'w = 0, we could start varying w and for each new w evaluate a new minimum of p(v). We would then find that the maximum value of all the minimum values evaluated is A,. This result can be generalized to the following principle, called the minimax characterization of eigenvalues, A, = max
. VTAv {min — vy
r=1,....n
withv satisfying v'w; = Ofori = 1,...,r
(2.126)
— 1, r = 2. In (2.126) we choose vectors w;,
i=1,...,r — 1, and then evaluate the minimum of p(v) with v subject to the condition viw, = 0,i = 1, ..., r — 1. After calculating this minimum we vary the vectors w; and always evaluate a new minimum, The maximum value that the minima reach is A,. The proof of (2.126) is as follows. Let v=2
A
i=1
2.127)
and evaluate the righthand side of (2.126), which we call R, R
=
{min
max
P+ afA oA o+ kA, ﬂ__l__z__(_!__z_ﬁz_lj__—z(_x—} airt  +tartarmt ot
(2'128)
The coefficients o; must satisfy the conditions
wXav=0
j=1,...,r—1
(2.129)
i=1
Rewriting (2.128), we obtain a%(Ar _
R = max
.
{min  A,
+.0+

A1)
+
a72'+1(Ar
St

a;—l(Ar
Ar+1)
+
Faltal
But we can now see that for the condition a1
R =A
= a2

.+
A,_l) a%(Ar
ot

An)
(2.130)
=    = @, = 0, we have
(2.131)
and the condition in (2.129) can still be satisfied by a judicious choice for a,. On the other hand, suppose that we now choose w; = v;forj = 1, ..., r — 1. This would require that oy =0forj=1,...,r — 1, and consequently we would have R = A,, which completes
the proof. A most important property that can be established using the minimax characterization of eigenvalues is the eigenvalue separation property. Suppose that in addition to the
64
Vectors, Matrices, and Tensors
Chap. 2
problem Av = Av, we consider the problems Almylm) = ) (miylm)
(2.132)
.where A™ is obtained by omitting the last m rows and columns of A. Hence A™ is a squaresymmetric matrix of order (n — m). Using also the notation A® = A, A® = ), v® = v, the eigenvalue separation property states that the eigenvalues of the problem Al Dylm+) = \(m+Dym+D) geparate the eigenvalues of the problem in (2.132); i.e., we have A'(l m )
0
;
P
EA(L.1) = Ro
d)
u(x,0) =0
and the initial conditions are
()
du a_t(x’
0)

0
With the conditions in (d) and () the formulation of the problem is complete, and (c) can be
solved for the displacement response of the rod.
Although we considered in these examples specific problems that are governed by elliptic, parabolic, and hyperbolic differential equations, the problem formulations illustrate in a quite general way some basic characteristics. In elliptic problems (see Exam
110
Some Basic Concepts of Engineering Analysis
Chap. 3
ple 3.15) the values of the unknown state variables (or their normal derivatives) are given
on the boundary. These problems are for this reason also called boundary value problems, where we should note that the solution at a general interior point depends on the data at every point of the boundary. A change in only one boundary value affects the complete solution; for instance, in Example 3.15 the complete solution for ¢ depends on the actual value of hy. Elliptic differential equations generally govern the steadystate response of systems. Comparing the governing differential equations given in Examples 3.15 to 3.17 it is noted that in contrast to the elliptic equation, the parabolic and hyperbolic equations (Examples 3.16 and 3.17, respectively) include time as an independent variable and thus define propagation problems. These problems are also called initial value problems because the solution depends on the initial conditions. We may note that analogous to the derivation of the dynamic equilibrium equations of lumpedparameter models, the governing differential equations of propagation problems are obtained from the steadystate equations by including the “resistance to change” (inertia) of the differential elements. Conversely, the
parabolic and hyperbolic differential equations in Examples 3.16 and 3.17 would become elliptic equations if the timedependent terms were neglected. In this way the initial value problems would be converted to boundary value problems with steadystate solutions. We stated earlier that the solution of a boundary value problem depends on the data at all points of the boundary. Here lies a significant difference in the analysis of a propagation problem, namely, considering propagation problems the solution at an interior point may depend only on the boundary conditions of part of the boundary and the initial conditions over part of the interior domain. 3.3.2 Variational Formulations
The variational approach of establishing the governing equilibrium equations of a system was already introduced as an alternative to the direct approach when we discussed the analysis of discrete systems (see Section 3.2.1). As described, the essence of the approach is to calculate the total potential IT of the system and to invoke the stationarity of I, i.e., 611 = 0, with respect to the state variables. We pointed out that the variational technique can be effective in the analysis of discrete systems; however, we shall now observe that the variational approach provides a particularly powerful mechanism for the analysis of continuous systems. The main reason for this effectiveness lies in the way by which some boundary conditions (namely, the natural boundary conditions defined below) can be generated and taken into account when using the variational approach. To demonstrate the variational formulation in the following examples, we assume that the total potential IT is given and defer the description of how an appropriate I can be determined until after the presentation of the examples. The total potential IT is also called the functional of the problem. Assume that in the functional the highest derivative of a state variable (with respect to a space coordinate) is of order m; i.e., the operator contains at most mthorder derivatives. We call such a problem a C™! variational problem. Considering the boundary conditions of the problem, we identify two classes of boundary conditions, called essential and natural boundary conditions. The essential boundary conditions are also called geometric boundary conditions because in structural mechanics the essential boundary conditions correspond to prescribed
Sec. 3.3
Solution of ContinuousSystem Mathematical Models
1M
displacements and rotations. The order of the derivatives in the essential boundary conditions is, in a C™! problem, at most m — 1. The second class of boundary conditions, namely, the natural boundary conditions, are also called force boundary conditions because in structural mechanics the natural boundary conditions correspond to prescribed boundary forces and moments. The highest derivatives in these boundary conditions are of order m to 2m — 1. We will see later that this classification of variational problems and associated boundary conditions is most useful in the design of numerical solutions. In the variational formulations we will use the variational symbol 8, already briefly employed in (3.1). Let us recall some important operational properties of this symbol; for more details, see, for example, R. Courant and D. Hilbert [A]. Assume that a function F for
a given value of x depends on v (the state variable), dv/dx, . . ., dPv/dx?, where p = 1,2, .... Then the first variation of F is defined as oF
=
oF
oF
—
dv dv +
dv —
d(dv/dx) 6
e
+
—
oF
dfo —_
.
(dPv/dxP) 6
(3.72)
This expression is explained as follows. We associate with v(x) a function € n(x) where € is a constant (independent of x) and 7(x) is an arbitrary but sufficiently smooth function that is zero at and corresponding to the essential boundary conditions. We call n(x) a variation in v, that is n(x) = 8v(x) [and of course € 17(x) is then also a variation in v] and also have for the required derivatives
that is, the variation of a derivative of v is equal to the derivative of the variation in v. The expression (3.7a) then follows from evaluating d F[o
P
P
+en,—%€—@,...,Lod+P—€l)i
8F = lim
 F(o,%,...,%)
 x
(3.7b)i
Considering (3.7a) we note that the expression for 6F looks like the expression for the total differential dF; that is, the variational operator & acts like the differential operator with respect to the variables v, dv/dx, ..., d°v/dx". These equations can be extended to multiple functions and state variables, and we find that the laws of variations of sums, products, and so on, are completely analogous to the corresponding laws of differentiation. For example, let F and Q be two functions possibly dependent on different state variables; then ⊳
3(F+Q)=3F+5Q;
S(FQ)=(3F3Q+F(3Q);
3(ﬂ"=n(F)""5F
In our applications the functions usually appear within an integral sign; and so, for example, we also use 6 ∫ F(x) dx = ∫ 8F(x)dx We shall employ these rules extensively in the variational derivations and will use one important condition (which corresponds to the properties of 7 stated earlier), namely, that the variations of the state variables and of their (m — 1)st derivatives must be zero at and
112
Some Basic Concepts of Engineering Analysis
Chap. 3
corresponding to the essential boundary conditions, but otherwise the variations can be arbitrary. Consider the following examples. EXAMPLE 3.18: in
Example
The functional governing the temperature distribution in the slab considered
3.16
iS
0
2
6x
o
0
a
and the essential boundary condition is 0. =6,
where
6 = 60, 1)
(b)
and
= O(L, t)
q” is the heat generated per unit volume, and otherwise the same notation as in Example 3.16 is used. Invoke the stationarity condition on II to derive the governing heat conduction equation and the natural boundary condition. This is a C? variational problem; i.e., the highest derivative in the functional in (a) is of order 1, or m =
1. An essential boundary condition, here given in (b), can therefore correspond
only to a prescribed temperature, and a natural boundary condition must correspond to a prescribed temperature gradient or boundary heat flow input. To invoke the stationarity condition 811 = 0, we can directly use the fact that variations and differentiations are performed with the same rules. That is, using (3.7a) we obtain
Ly
oa0\[.
a0
J; (kEC)(Sax> dx — ∫ 60qu—800qo—0
©
where also 8(86/dx) = 986/ dx . The same result is also obtained when using (3.7b), which gives here
con) v~ [ (2 o [0 L
8I1 = lim e—0
1
(a0
d
_"_____.____________—__ €
(36Y
"1
5
dx
L
—
By
€ 1
a0an
L
,
377)2]
A/
≔∫
o
—
—
B dx
eml.
€
€—0
L
dx
J~L
—
0 o k_a__ndx_f X 0x
L andx_nqu
o
where 1o = 7.0 and we would now substitute 88 for 7. Now using integration by parts, we obtain from (c) the following equation: 2
−∫ (ka—o )60dx+ka—0ox o
dx?
@
x=L
®
2The divergence theorem is used (see Examples 4.2 and 7.1).
—_—
®
Sec. 3.3
Solution of ContinuousSystem Mathematical Models
13
To extract from (d) the governing differential equation and natural boundary condition, we use the argument that the variations on @ are completely arbitrary, except that there can be no variations on the prescribed essential boundary conditions. Hence, because 6, is prescribed, we
have 86, = 0 and term ) in (d) vanishes. Considering next terms (D) and (3), assume that 86, = 0 but that 88 is otherwise nonzero (except at x = 0, where we have a sudden jump to a zero value). If (d) is to hold for any nonzero 80, we need to have?
©
k——+4q°=0
Conversely, assume that 66 is zero everywhere except at x = 0; i.e., we have 86, # O; then (d) is valid only if a0 k—ox o
+go=0 9o
f)
The expression in (f) represents the natural boundary condition. The governing differential equation of the propagation problem is obtained from (e), specifying here that
q° = —pc
(8
3%0
Hi ence (e) reduces d to
36
k—; %2 = pc pc— o
We may note that until the heat capacity effect was introduced in the formulation in (g), the equations were derived as if a steadystate problem (and with g® timedependent a pseudo steadystate problem) was being considered. Hence, as noted previously, the formulation of the propagation problem can be obtained from the equation governing the steadystate response by simply taking into account the timedependent “inertia term.”
EXAMPLE 3.19: The functional and essential boundary condition governing the wave propagation in the rod considered in Example 3.17 are 1'[
=
L1
ou\? _—
∫∅ 2EA
and
dx

£
J;’ uf®
B
_
dx — u,R
(a)
u =0
®)
where the same notation as in Example 3.17 is used, 4o = u(0, #), u, = u(L, t), and /@ is the body
force per unit length of the rod. Show that by invoking the stationarity condition on II the governing differential equation of the propagation problem and the natural boundary condition can be derived. We proceed as in Example 3.18. The stationarity condition 8II = 0 gives L
u
6u>
J;’ (EAEC)(SEC
dx
L
5
J;’ ouf®dx
— uR=20
Writing d8u/dx for 8 (9u/dx), recalling that EA is constant, and using integration by parts yields L
3 — R] Su, — EAZ
2
−∫ o (EAa—"2+ f") Bu dx + [EAa—u ax ax
x=L
a0x
3We in effect imply here that the limits of integration are not 0 to L but 0* to L~
Suo x=0
=
()
114
Some Basic Concepts of Engineering Analysis
Chap. 3
To obtain the governing differential equation and natural boundary condition we use, in essence, the same argument as in Example 3.18; i.e., since 6u, is zero but Su is arbitrary at all other points, we must have
o*u
EA ‘("E + f5=0
(©)
ou EA —ox , =R
d)d
and
In this problem we have f2 = —Ap 9%u/0t* and hence (c) reduces to the problemgoverning differential equation
The natural boundary condition was stated in (d). Finally, it may be noted that the problem in (a) and (b) is a C° variational problem; i.e.,
m = 1 in this case. EXAMPLE 3.20:
The functional governing static buckling of the column in Fig. E3.20 is IT=
1 ("
(d*wY
EJ;, EI(d—x2>
P
dx
2J;
(" (dwY
1,
(E)
dx + EkWL
(a)
where w, = wx=, and the essential boundary conditions are dw Wx=o
=u,
E
o
=
0
(b)
Invoke the stationarity condition 811 = 0 to derive the problemgoverning differential equation and the natural boundary conditions. Flexural stiffness
El
Spring stiffness
Figure E3.20
Column subjected to a compressive load
This problem is a C! variational problem, i.e., m = 2, because the highest derivative in the functional is of order 2. The stationarity condition 811 = 0 yields L
∫
L
EIw”Sw”dxo
PJ
w’6w’dx+kwL8wL=0 o
Sec. 3.3
Solution of ContinuousSystem
Mathematical Models
where we use the notation w' = dw/dx, and so on. But 8w”
115
= d(8w')/dx, and EI is constant;
hence, using integration by parts, we obtain L
L
∫ EIw”Sw”dx=EIw”6w’(’;— EIJ o
w”’6w’dx o
If we continue to integrate by parts fy w
m
8w’ dx and also integrate by parts fy w' w' dx,
we obtain L
∫
(EIin+Pw”)6wdx+(EIw”Sw’)L—(Elw”Sw’)o
↼∅−−−−⋎−−−−∣↘−−∙⋎−−⇃
@
@
®
— [(Etw"™ + Pw') éw]. + [(EIw™ + Pw') 8w]o + kw, 6w = 0
@
®
©
®
Since the variations on w and w' must be zero at the essential boundary conditions, we have
dwo = 0 and 8w{ = 0. It follows that terms (3) and (5) are zero. The variations on w and w’ are arbitrary at all other points, hence to satisfy (c) we conclude, using the earlier arguments (see Example 3.18), that the following equations must be satisfied: term 1:
EIw"™ + Pw" =
(d)
term 2:
EW", =0
©
(EIw" + Pw' — kw)e—r =
(f)
terms 4 and 6:
The problemgoverning differential equation is given in (d), and the natural boundary conditions are the relations in (e) and (f). We should note that these boundary conditions correspond to the physical conditions of moment and shear equilibrium at x = L.
We have illustrated in the preceding examples how the problemgoverning differential equation and the natural boundary conditions can be derived by invoking the stationarity of the functional of the problem. At this point a number of observations should be made. First, considering a C™~! variational problem, the order of the highest derivative present in the problemgoverning differential equation is 2m. The reason for obtaining a derivative of order 2m in the problemgoverning differential equation is that integration by parts is employed m times. A second observation is that the effect of the natural boundary conditions was always included as a potential in the expression for I1. Hence the natural boundary conditions are implicitly contained in II, whereas the essential boundary conditions have been stated separately. Our objective in Examples 3.18 to 3.20 was to derive the governing differential equations and natural boundary conditions by invoking the stationarity of a functional, and for this purpose the appropriate functional was given in each case. However, an important question then arises: How can we establish an appropriate functional corresponding to a given problem? The two previous observations and the mathematical manipulations in Examples 3.18 to 3.20 suggest that to derive a functional for a given problem we could start with the governing differential equation, establish an integral equation, and then proceed backward in the mathematical manipulations. In this derivation it is necessary to use integration by parts, i.e., the divergence theorem, and the final check would be that the
116
Some
Basic Concepts of Engineering Analysis
Chap. 3
stationarity condition on the IT derived does indeed yield the governing differential equations. This procedure is employed to derive appropriate functionals in many cases (see Section 3.3.4 and Chapters 4 and 7, and for further treatment see, for example, R. Courant and D, Hilbert [A], S. G. Mikhlin [A], K. Washizu [B], and J. T. Oden and J. N. Reddy [A]). In this context, it should also be noted that in considering a specific problem, there does not generally exist a unique appropriate functional, but a number of functionals are applicable. For instance, in the solution of structural mechanics problems, we can employ the principle of minimum potential energy, other displacementbased variational formulations, the HuWashizu or HellingerReissner principles, and so on (see Section 4.4.2). Another important observation is that once a functional has been established for a certain class of problems, the functional can be employed to generate the governing equations for all problems in that class and therefore provides a general analysis tool. For example, the principle of minimum potential energy is general and is applicable to all problems in linear elasticity theory. Based simply on a utilitarian point of view, the following observations can be made in regard to variational formulations.
1. The variational method may provide a relatively easy way to construct the systemgoverning equations. This ease of use of a variational principle depends largely on the fact that in the variational formulation scalar quantities (energies, potentials, and so
on) are considered rather than vector quantities (forces, displacements, and so on).
2. A variational approach may lead more directly to the systemgoverning equations and boundary conditions. For example, if a complex system is being considered, it is of advantage that some variables that need to be included in a direct formulation are not considered in a variational formulation (such as internal forces that do no net work).
3. The variational approach provides some additional insight into a problem and gives an independent check on the formulation of the problem. 4. For approximate solutions, a larger class of trial functions can be employed in many cases if the analyst operates on the variational formulation rather than on the differential formulation of the problem; for example, the trial functions need not satisfy the natural boundary conditions because these boundary conditions are implicitly contained in the functional (see Section 3.3.4).
This last consideration has most important consequences, and much of the success of the finite element method hinges on the fact that by employing a variational formulation, a larger class of functions can be used. We examine this point in more detail in the next section and in Section 3.3.4. 3.3.3 Weighted Residual Methods; Ritz Method
In previous sections we have discussed differential and variational formulations of the governing equilibrium equations of continuous systems. In dealing with relatively simple systems, these equations can be solved in closed form using techniques of integration, separation of variables, and so on. For more complex systems, approximate procedures of solution must be employed. The objective in this section is to survey some classical techniques in which a family of trial functions is used to obtain an approximate solution. We
Sec. 3.3
Mathematical
Solution of ContinuousSystem
Models
117
shall see later that these techniques are very closely related to the finite element method of analysis and that indeed the finite element method can be regarded as an extension of these classical procedures. Consider the analysis of a steadystate problem using its differential formulation Ladd) =1
(3.8)
in which L,,, is a linear differential operator, ¢ is the state variable to be calculated, and » is the forcing function. The solution to the problem must also satisfy the boundary conditions =
B’[¢]
1, 2,
i =
qilat boundary 5;5
(39)
. o
We shall be concerned, in particular, with symmetric and positive definite operators that satisfy the symmetry condition ∫ (Lz,,,[u])odD= ∫ (Lam[v])u dD D
(3.10)
D
and the condition of positive definiteness (3.11)
∫ (Lz,,,[u])udD>0
where D is the domain of the operator and « and v are any functions that satisfy homogeneous essential and natural boundary conditions. To clarify the meaning of relations (3.8) to (3.11), we consider the following example. EXAMPLE 3.21: The steadystate response of the bar shown in Fig. E3.17 is calculated by solving the differential equation
o*u —
—EA A
=
(a)
0
ax2
subject to the boundary conditions
U=0 = 0;
ou
EA
i
R
(b)
Identify the operators and functions of (3.8) and (3.9) and check whether the operator L., is
symmetric and positive definite. Comparing (3.8) with (a), we see that in this problem 62
L2m=—EA'(;x—2‘,
b=u
r=0
Similarly, comparing (3.9) with (b), we obtain
B,
g
=1
B
a
2
= EA—; ax
7
=0
=R
To identify whether the operator L,,, is symmetric and positive definite, we consider the case R = 0. This means physically that we are concerned only with the structure itself and not
118
Some
Basic Concepts of Engineering Analysis
Chap. 3
with the loading applied to it. For (3.10) we have
=
ou
*
axv
0
—FA—
avl*
FAu
—
ox
“
x Ox
o
x
x?2
0
—
J'L
%0
EA—ud “ax
dx?
0
0
©
Since the boundary conditions are ¥ = v = Oatx = 0 and du/dx = dv/dx = Qatx = L, we have L
J;
—
0%u EA —— oz’ dx
L =
J;
—
E. A
820 ——— axt
dx
and the operator is symmetric. We can also directly conclude that the operator is positive definite because from (c) we obtain L
L
9%u
a
E
This Ritz analysis therefore yields the approximate solution
_ 129
u=—7x
0341 7
o = 129 — 0.682x;
X
,
(h)
0= m
HOT[fEm — s H@U vim)
— k™H™T] v
(4.26)
166
Formulation of the Finite Element Method
Chap. 4
In this case the vectors f*™ no longer include inertia and velocitydependent damping forces, U is a vector of the nodal point velocities (i.e., the first time derivative of U), and k'™ is the damping property parameter of element m. The equilibrium equations are, in this case,
MU + CU + KU = R
@.27)
where C is the damping matrix of the structure; i.e., formally,
c=2
,, KPHOTHE dye m
(4.28)
m)
=
C(m)
In practice it is difficult, if not impossible, to determine for general finite element assemblages the element damping parameters, in particular because the damping properties are frequency dependent. For this reason, the matrix C is in general not assembled from element damping matrices but is constructed using the mass matrix and stiffness matrix of the complete element assemblage together with experimental results on the amount of damping. Some formulations used to construct physically significant damping matrices are described in Section 9.3.3. A complete analysis, therefore, consists of calculating the matrix K (and the matrices M and C in a dynamic analysis) and the load vector R, solving for the response U from (4.17) [or U, U, U from (4.24) or (4.27)], and then evaluatmg the stresses using (4.12). We should emphasize that the stresses are simply obtained using (4.12)—hence only from the initial stresses and element displacements—and that these values are not corrected for externally applied element pressures or body forces, as is common practice in the analysis of frame structures with beam elements (see Example 4.5 and, for example, S. H. Crandall,
N. C. Dahl, and T. J. Lardner [A]). In the analysis of beam structures, each element represents a onedimensional stress situation, and the stress correction due to distributed loading is performed by simple equilibrium considerations. In static analysis, relatively long beam elements can therefore be employed, resulting in the use of only a few elements (and degrees of freedom) to represent a frame structure. However, a similar scheme would require, in general two and threedimensional finite element analysis, the solution of boundary value problems for the (large) element domains used, and the use of fine meshes for an accurate prediction of the displacements and strains is more effective. With such fine discretizations, the benefits of even correcting approximately the stress predictions for the effects of distributed element loadings are in general small, although for specific situations of course the use of a rational scheme can result in notable improvements. To illustrate the above derivation of the finite element equilibrium equations, we consider the following examples. EXAMPLE 4.5: Establish the finite element equilibrium equations of the bar structure shown in Fig. E4.5. The mathematical model to be used is discussed in Examples 3.17 and 3.22. Use the twonode bar element idealization given and consider the following two cases: 1. Assume that the loads are applied very slowly when measured on the time constants
(natural periods) of the structure, 2. Assume that the loads are applied rapidly. The structure is initially at rest.
Crosssectional
A={(1+n/40)
area A= 1cm?
————
1001,
100cm —— 80cm
E = Young's modulus
−−−⊣
p = mass density
100 cm ——— 80 cm −−⊣
(b)Elementassemblageinglobalsystem
(c)Element1,fx"=f2(l*)N/cm3 9 cm?
fi(t), R(t) 1.0
1
2
3
4
5
Time
(e} Time variation of loads
Figure E4.5
Twoclement bar assemblage
167
168
Formulation of the Finite Element Method In
Chap. 4
the formulation of the finite element equilibrium equations we employ the general
equations (4.8) to (4.24) but use that the only nonzero stress is the longitudinal stress in the bar.
Furthermore, considering the complete bar as an assemblage of 2 twonode bar elements corresponds to assuming a linear displacement variation between the nodal points of each element. The first step is to construct the matrices H™ and B™ for m = 1, 2. We recall that although the displacement at the left end of the structure is zero, we first include the displacement at that surface in the construction of the finite element equilibrium equations.
Corresponding to the displacement vector UT = [U;
)
= 
B
1 ———
[
100
100
2)
1]
80)
( 1 —_——
=
B
0]
§  X [0
Us), we have
1
® H
U,
[0
80
80
1
80]
The material property matrices are
C® = E;
C?=E
where E is Young’s modulus for the material. For the volume crosssectional areas of the elements. We have
AV = lcm’
integrations we need the
AP = (1 + 4i0>2 cm’?
When the loads are applied very slowly, a static analysis is required in which the stiffness matrix K and load vector R must be calculated. The body forces and loads are given in Fig. E4.5. We therefore have
K
=
≺⋯⋮∫∘
1%
11 ] — —  —=
+
100 [ 100 100 0] dx
1  + =} =
%
80 [0
40)
⋮∫∘ (1
0
80 so] d
% [
or
11 = 
K=WO
e[ =%
11 1
0 1
0
135] © +?4_0
0
00
24 24
24 154
0 13
0
13
13
0
0
0
1
1
o0 1
1
(a
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
169
and also,
x 1  100 _
Rp = {(1) J;
100
0 )
x
100
(1) dx + J;
%
(1
x
+ 4—())
V1
0 ]
=3

. 80
1
(—(—)>dx}fz(t)
30
150
186  f2(¢)
(b
68 0 Re =
0
A@)
©
100
To obtain the solution at a specific time ¢*, the vectors R and Rc must be evaluated corresponding to ¢*, and the equation KU
t=t* =
RBt=t"
+
RCz=z*
(d)
yields the displacements at *. We should note that in this static analysis the displacements at time t* depend only on the magnitude of the loads at that time and are independent of the loading history. Considering now the dynamic analysis, we also need to calculate the mass matrix. Using the displacement interpolations and (4.25), we have
100 M=(1)pj
0
]  X 100 x
5 [(1——)
100
5 _—
100
0] dx
100
0 0 +pJ;
80
x
\?
x
x
(1+4—0>1—%[0
(1_8_0>
x
%]dx
x 80
Hence
200
100
0
M=‘63 100
584
336
0
336
1024
Damping was not specified; thus, the equilibrium equations now to be solved are
MU(7) + KU(?) = Ra(r) + Rc(r)
©
where the stiffness matrix K and load vectors R and Rc have already been given in (a) to (c).
Using the initial conditions
Ulwo=0;
Ulwo=0
()
these dynamic equilibrium equations must be integrated from time O to time ¢* in order to obtain the solution at time z* (see Chapter 9).
170
Formulation of the Finite Element Method
Chap. 4
To actually solve for the response of the structure in Fig. E4.5(a), we need to impose U, = 0 for all time . Hence, the equations (d) and (e) must be amended by this condition (see Section 4.2.2). The solution of (d) and (¢) then yields Ux(z), Us(#), and the stresses are obtained
using
@
m=1,2
7w = C"B"U();,
These stresses will be discontinuous between the elements because constant element strains are assumed. Of course, in this example, since the exact solution to the mathematical model can be computed, stresses more accurate than those given by (g) could be evaluated within each element. In static analysis, this increase in accuracy could simply be achieved, as in beam theory, by adding a stress correction for the distributed element loading to the values given in (g). However, such a stress correction is not straightforward in general dynamic analysis (and in any two and threedimensional practical analysis), and if a large number of elements is used to represent the structure, the stresses using (g) are sufficiently accurate (see Section 4.3.6). EXAMPLE 4.6: Consider the analysis of the cantilever plate shown in Fig. E4.6. To illustrate the analysis technique, use the coarse finite element idealization given in the figure (in a practical analysis more finite elements must be employed (see Section 4.3). Establish the matrices H?,
B, and C?, The cantilever plate is acting in plane stress conditions. For an isotropic linear elastic material the stressstrain matrix is defined using Young’s modulus E and Poisson’s ratio v (see Table 4.3),
Cco = 1 E
.
 v
1
v
0
v
1
1 0
00
The displacement transformation matrix H® displacements to the nodal point displacements, [u(xv
)’)]
@
 v 2
=
of element 2 relates the element internal
H(Z)U
(a)
o(x, y) where U is a vector listing all nodal point displacements of the structure, Ur=[U,
U.
Us
Us
...
Uy
Us)
(b)
(As mentioned previously, in this phase of analysis we are considering the structural model without displacement boundary conditions.) In considering element 2, we recognize that only the
displacements at nodes 6, 3, 2, and 5 affect the displacements in the element. For computational purposes it is convenient to use a convention to number the element nodal points and corresponding element degrees of freedom as shown in Fig E4.6(c). In the same figure the global structure degrees of freedom of the vector U in (b) are also given. To derive the matrix H® in (a) we recognize that there are four nodal point displacements each for expressing u(x, y) and o(x, y). Hence, we can assume that the local element displacements # and v are given in the following form of polynomials in the local coordinate variables x and y: u(x,y) = ai + aax + azy + auxy
o(x,y) = B1 + Bax + Bsy + Baxy
©
Sec. 4.2
171
Formulation of the DisplacementBased Finite Element Method
o  @  O Element
Thickness t=0.1
2cm
X U
4
7
U
⊢−−≀ ∁⊓↿↠⊢−≀ ∁⊓↿−−⊣
4cm
(b) Finite element idealization (plane stress condition)
(a) Cantilever plate
.
2= Us
vi=
I
U12
upym U
Uy m Uny
El
Element nodal point 4 =
@
ement
structure nodal point 5
5 4
3
uz=Us
— iy = Uy
≀∁⊓↿⇁∙↑ V35U4
Vi =
U1o
(c) Typical twodimensional fournode element defined in local coordinate system
Figure E4.6
Finite element plane stress analysis
The unknown coefficients ai, . . . , Bs, which are also called the generalized coordinates, will be expressed in terms of the unknown element nodal point displacements u;, . .., us and ¥1, . . . , 4. Defining 0
=[w
u
us
us
 vr
vz
03
04
d)
we can write (c) in matrix form:
[ ] = o=
0
Formulation of the Finite Element Method
172
where
P
and
o
_ db=[1
0of]
_[®[0
=[e;
v
a5
x
i P
y
B
Chap. 4
xy] B
B
Equation (¢) must hold for all nodal points of the element; therfore, using (d), we have
i = Aa
)
A
21]
A= [0‘
in which
ATl
and
1 1
1 1
o1
o
1
1
1
1 1
1 1 1
Solving from (f) for a and substituting into (€), we obtain
H = ®A™!
(®
where the fact that no superscript is used on H indicates that the displacement interpolation matrix is defined corresponding to the element nodal point displacements in (d),
H = 1[(1 +x00+y 4
Ax010+y
0
T090y
0
0
1+x01y
0
0
0
0
0
F+90+y) AH1+y) AH1y
1+0y)
]
b
The displacement functions in H could also have been established by inspection. Let H; be the (i, j)th element of H; then Hy; corresponds to a function that varies linearly in x and y [as required in ()], isunity atx = 1,y = 1, and is zero at the other three element nodes. We discuss the construction of the displacement functions in H based on these thoughts in Section 5.2. With H given in (h) we have Us
U H?
= [0 0
Us
Uy { 0 O
U3
U
V2
U}
Us
Usi
V4
U
Up
1} ⋮
U
Ui
0 O
0
⋮
H14H13
0

Hu
Hy
v; ‘
Figure E4.8
Chap. 4
Element 2 of bar analyzed
in Example 4.5
These are the exact element internal displacements. The element end forces required to subject the bar to these displacements are du
ki 12 = —EA — x,y du
©
k2 22 = EA — P Substituting from (b) into (c) we have
3E ks
=
3E
%;
kiz
=
_%
Hence we have, using the symmetry of the element matrix and equilibrium to establish &,; and ki,
3
K= %E[_l
1
1
1]
@
The same result is of course obtained using the principle of virtual displacements with the displacement (b). We
note that the stiffness coefficient in (d) is smaller than the corresponding
value
obtained in Example 4.5 (3 E/80 instead of 13 E/240). The finite element solution in Example 4.5 overestimates the stiffness of the structure because the assumed displacements artificially constrain the motion of the material particles (see Section 4.3.4). To check that the internal equilibrium is indeed not satisfied, we substitute the finite element solution (given by the displacement
assumption in Example 4.5) into (a) and obtain
£l d =
(1+35) ) *
V1 + =) —
The solution of truss and beam structures, using the exact displacements corresponding to unit nodal point displacements and rotations to evaluate the stiffness matrices, gives analysis results that for the selected mathematical model satisfy all three requirements of mechanics exactly: differential equilibrium for every point of the structure (including nodal point equilibrium), compatibility, and the stressstrain relationships. Hence, the exact (unique) solution for the selected mathematical model is obtained.
We may note that such an exact solution is usually pursued in static analysis, in which
the exact stiffness relationships are obtained as described in Example 4.8, but an exact
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
177
solution is much more difficult to reach in dynamic analysis because in this case the distributed mass and damping effects must be included (see, for example, R. W. Clough and J. Penzien [A]). However, although in a general (static or dynamic) finite element analysis, differential equilibrium is not exactly satisfied at all points of the continuum considered, two important properties are always satisfied by the finite element solution using a coarse or a fine mesh. These properties are (see Fig. 4.2) 1. Nodal point equilibrium 2. Element equilibrium.
Element m
Sum of forces F™ equilibrate ~————»!
\\
⊮∶∘↾∁⊜⊱∣⋍≺⋒⊞⊜
↟ \/
~—
——
externally applied loads
!
Figure 4.2
Nodal point and element equilibrium in a finite element analysis
in equilibrium
178
Formulation of the Finite Element Method
Chap. 4
Namely, consider that a finite element analysis has been performed and that we calculate for each finite element m the element nodal point force vectors F
= J'
B Tgm) gyim Vim)
(4.29)
where 7™ = C™e€™, Then we observe that according to property 1, At any node, the sum of the element nodal point forces is in equilibrium with the externally applied nodal loads (which include all effects due to body forces, surface tractions, initial stresses, concentrated loads, inertia and damping forces, and reactions). And according to property 2, Each element m is in equilibrium under its forces F™,
Property and we have
1 follows simply because
(4.27) expresses the nodal point equilibrium
Y F™ = KU
(4.30)
The element equilibrium stated in property 2 is satisfied provided the finite element displacement interpolations in H™ satisfy the basic convergence requirements, which include the condition that the element must be able to represent the rigid body motions (see Section 4.3). Namely, let us consider element m subjected to the nodal point forces F™ and
impose virtual nodal point displacements corresponding to the rigid body motions. Then for each virtual element rigid body motion with nodal point displacements @, we have ﬁ' TF‘(’")
=
J’
⋁≼⋯≱
B2
⋜≺⋯≻⊺⊤≺⋯≻∠↿⋁≺⋯⋟∶∘
qym = ∫ V(m)
because here € ™ = 0. Using all applicable rigid body motions we therefore find that the forces F™ are in equilibrium. Hence, a finite element analysis can be interpreted as a process in which
1. The structure or continuum is idealized as an assemblage of discrete elements connected at nodes pertaining to the elements. 2. The externally applied forces (body forces, surface tractions, initial stresses, concentrated loads, inertia and damping forces, and reactions) are lumped to these nodes using the virtual work principle to obtain equivalent externally applied nodal point forces. 3. The equivalent externally applied nodal point forces (calculated in 2) are equilibrated by the element nodal point forces that are equivalent (in the virtual work sense) to the element internal stresses; i.e., we have
2 Fm =R m
4. Compatibility and the stressstrain material relationship are satisfied exactly, but instead of equilibrium on the differential level, only global equilibrium for the com
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
179
plete structure, at the nodes, and of each element m under its nodal point forces F™ is satisfied. Consider the following example. EXAMPLE
4.9:
The finite element solution to the problem in Fig. E4.6, with P =
100, E =
2.7 X 105, v = 0.30, r = 0.1, is given in Fig. E4.9. Clearly, the stresses are not continuous between elements, and equilibrium on the differential level is not satisfied. However,
1. Show that £ F™ = R and calculate the reactions. 2. Show that the element forces F for element 4 are in equilibrium. The fact that 5 F™
= R follows from the solution of (4.17), and R consists of the sum
of all nodal point forces. Hence, this relation can also be used to evaluate the reactions. Referring to the nodal point numbering in Fig. E4.6(b), we find for node 1: reactions R, = 100.15
R, = 41.36 for node 2:
reactions R, = 2.58 — 2.88 = —0.30 R, = 16.79 + 5.96 = 22.74 (because of rounding) for node 3: reactions R, =
—99.85
R, = 35.90 for node 4:
horizontal force equilibrium: —42.01
+ 42.01 = 0
vertical force equilibrium: —22.90 + 22.90 = 0 for node 5: horizontal force equilibrium: —60.72 — 12.04 + 44.73 + 28.03 = 0 vertical force equilibrium: —35.24 — 35.04 + 19.10 + 51.18 = 0 for node 6:
horizontal force equilibrium: 57.99 — 57.99 = 0 vertical force equilibrium: —6.81
+ 6.81 = 0
And for nodes 7, 8, and 9, force equilibrium is obviously also satisfied, where at node 9 the
element nodal force balances the applied load P = 100.
Finally, let us check the overall force equilibrium of the model: horizontal equilibrium:
100.15 — 030 — 99.85 = 0 vertical equilibrium:
41.36 + 22.74 + 3590 — 100 =0
Formulation of the Finite Element Method
180 1129
Chap. 4
1123
o) 95.92
—972.1
931.4
{a) Exploded view of elements showing stresses 7\7. Note the stress discontinuities between elements and the nonzero stresses along the free edges
zim)
y
(b) Exploded view of elements showing stresses t(y"y” Figure E4.9
Solution results for problem considered in Example 4.6 (rounded to digits shown)
426.9

' y
420.3
.
x
5
‘
467.8
35.90 99.85
6.81 57.99
——
—
—288
44'73
↑ 5.96
2.58
100.0
57.99
0
e——
——
2803
35.04
29.97
↑ 51.18
—35 24
↑ 42.01
19.10
ok 
16.79
6.81
~42.01
60 72 12.04
↑ 41.36
↑ ~22.90
↑
22.90
I
E f
→−
↓
→
~29.97
↑
(d) Exploded view of elements showing element nodal point forces equivalent {in the virtual work sense) to the element stresses. The nodal point forces are at each node in equilibrium with the applied forces (including the reactions) Figure E4.9
(continued)
181
182
Formulation of the Finite Element Method
Chap. 4
moment equilibrium (about node 2):
—100 X 4 + 100.15 X2
+ 9985
X2
=0
It is important to realize that this force equilibrium will hold for any finite element mesh, however
coarse the mesh may be, provided properly formulated elements are used (see Section 4.3). Now consider element 4: horizontal equilibrium:
0 — 57.99 + 28.03 + 29.97 = 0 (because of rounding) vertical equilibrium:
~100 + 6.81 + 51.18 + 42,01 = 0 moment equilibrium (about its local node 3):
—100 X 2 + 5799
X2 + 4201
x2=0
Hence the element nodal forces are in equilibrium.
Element Local and Structure Global Degrees of Freedom
The derivations of the element matrices in Example 4.6 and 4.7 show that it is expedient to first establish the matrices corresponding to the local element degrees of freedom. The construction of the finite element matrices, which correspond to the global assemblage degrees of freedom [used in (4.19) to (4.25)] can then be directly achieved by identifying the global degrees of freedom that correspond to the local element degrees of freedom. However, considering the matrices H™, B™, K™, and so on, corresponding to the global assemblage degrees of freedom, only those rows and columns that correspond to element degrees of freedom have nonzero entries, and the main objective in defining these specific matrices was to be able to express the assemblage process of the element matrices in a theoretically elegant manner. In the practical implementation of the finite element method, this elegance is also present, but all element matrices are calculated corresponding only to the element degrees of freedom and are then directly assembled using the correspondence between the local element and global assemblage degrees of freedom. Thus, with only the element local nodal point degrees of freedom listed in @i, we now write (as in Example 4.6)
where the entries in the vector u are the element displacements measured in any convenient local coordinate system. We then also have
Considering the relations in (4.31) and (4.32), the fact that no superscript is used on
the interpolation matrices indicates that the matrices are defined with respect to the local element degrees of freedom. Using the relations for the element stiffness matrix, mass
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
183
matrix, and load vector calculations as before, we obtain
K=
J' B’CB dV
(4.33)
v
M = ∫⋁ pHTHdV
(4.34)
Rs = ∫⋁ HTdev
(4.35)
Rs = J; HfS 4s
(4.36)
R = ∫⋁ BT‘r’dV
(4.37)
where all variables are defined as in (4.19) to (4.25), but corresponding to the local element degrees of freedom. In the derivations and discussions to follow, we shall refer extensively
to the relations in (4.33) to (4.37). Once the matrices given in (4.33) to (4.37) have been
calculated, they can be assembled directly using the procedures described in Example 4.11 and Chapter 12. In this assemblage process it is assumed that the directions of the element nodal point displacements 11 in (4.31) are the same as the directions of the global nodal point displacements U, However, in some analyses it is convenient to start the derivation with element nodal point degrees of freedom fi that are not aligned with the global assemblage degrees of freedom. In this case we have ~
u = Hi
(4.38)
i=Ta
(4.39)
and
where the matrix T transforms the degrees of freedom @ to the degrees of freedom i and (4.39) corresponds to a firstorder tensor transformation (see Section 2.4); the entries in column j of the matrix T are the direction cosines of a unit vector corresponding to the jth degree of freedom in @ when measured in the directions of the @ degrees of freedom. Substituting (4.39) into (4.38), we obtain
H = HT
(4.40)
Chap. 4
Formulation of the Finite Element Method
184
Thus, identifying all finite element matrices corresponding to the degrees of freedom @ with a curl placed over them, we obtain from (4.40) and (4.33) to (4.37),
K =TKT;,
M =TMT
Rs = T'Ry;
Rs = T'R;
@.41) R, = T'R,
We note that such transformations are also used when boundary displacements must be imposed that do not correspond to the global assemblage degrees of freedom (see Section 4.2.2). Table 4.1 summarizes some of the notation that we have employed. We demonstrate the presented concepts in the following examples. TABLE 4.1
@
Summary of some notation used
u™=H"Uoru™ =HU where u™ = displacements within element m as a function of the element coordinates U = nodal point displacements of the total element assemblage [from equation (4.17)
onward we simply use U]. (®)
u=Hi where u = u™ and it is implied that a specific element is considered
@ = nodal point displacements of the element under consideration; the entries of # are those displacements in U that belong to the element. u=Hi
where @i = nodal point displacements of an element in a coordinate system other than the global system (in which U is defined).
[u(x)]=%
cos po H¥”
p=1 be
w(x,y, 0) = > p=1
sin po HW?
®)
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
211
where for the triangular elements, referring to Example 4.17, H=[1
x
ylA{
©
and the @, %, and W” are the element unknown generalized nodal point displacements corre
sponding to mode p. We should note that we superimpose in (b) the response measured in individual harmonic displacement distributions. Using (b), we can now establish the straindisplacement matrix of the element. Since we are dealing with a threedimensional stress distribution, we use the expression for threedimensional strain distributions in cylindrical coordinates: du
or do
By €=
u
Low
r
r a0
du ——
do —
+
@
dy
Oor
ow
l1dou
w
o
rae
r
Y
Ve
wdy 1r a0 where
€=[e
€
€
Vel
©
Substituting from (b) into (d) we obtain a straindisplacement matrix B, for each value of
p, and the total strains can be thought of as the superposition of the strain distributions contained in each harmonic. The unknown nodal point displacements can now be evaluated using the usual procedures. The equilibrivm equations corresponding to the generalized nodal point displacements U?, V¥, Wi i=1,...,N(Nisequal to the total number of nodes) andp = 1, . . ., p. are evaluated as given in (4.17) to (4.22), where we now have
U=
[U7
U7
...
U]
(f)
and
U=
Vi WU
Wi
@)
In the calculations of K and Rs we note that because of the orthogonality properties 27
∫
sinnBsinm0d0=0
n#m
0 (h)
27
cosnBcosm0d0=0
∫
n+m
0
the stiffness matrices corresponding to the different harmonics are decoupled from each other. Hence, we have the following equilibrium equations for the structure:
KU
= R%
p=1...,p
(@)
where K? and R are the stiffness matrix and load vector corresponding to the pth harmonic.
212
Formulation of the Finite Element Method
Chap. 4
Solution of the equations in (i) gives the generalized nodal point displacements of each element, and (b) then yields all element internal displacements. In the above displacement solution we considered only the symmetric load contributions. But an analogous analysis can be performed for the antisymmetric load harmonics of (a) by simply replacing in (b) to (i) all sine and cosine terms by cosine and sine terms, respectively. The complete structural response is then obtained by superimposing the displacements corresponding to all harmonics. Although we have considered only surface loading in the discussion, the analysis can be extended using the same approach to include body force loading and initial stresses. Finally, we note that the computational effort required in the analysis is directly proportional to the number of load harmonics used. Hence, the solution procedure is very efficient if the loading can be represented using only a few harmonics (e.g., wind loading) but may be inefficient when many harmonics must be used to represent the loading (e.g., a concentrated force).
4.2.4 Lumping of Structure Properties and Loads
A physical interpretation of the finite element procedure of analysis as presented in the previous sections is that the structure properties—stiffness and mass—and the loads, internal and external, are lumped to the discrete nodes of the element assemblage using the virtual work principle. Because the same interpolation functions are employed in the
fe
elong
edge 2—1
3.0 20 1.0
"
3 fB g
plane stress element,
18 elong
thickness = 0.5
edges 41
10
sgmg34
10
30
4
end 32
e X
R=10.5 0.5 0.5 0.5 0.0 0.0 0.0 0.0 Figure 4.6
Body force distribution and corresponding lumped body force vector R; of a
rectangular element
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
213
calculation of the load vectors and the mass matrix as in the evaluation of the stiffness matrix, we say that “consistent” load vectors and a consistent mass matrix are evaluated. In this case, provided certain conditions are fulfilled (see Section 4.3.3), the finite element solution is a Ritz analysis. It may now be recognized that instead of performing the integrations leading to the consistent load vector, we may evaluate an approximate load vector by simply adding to the actually applied concentrated nodal forces Rc additional forces that are in some sense equivalent to the distributed loads on the elements. A somewhat obvious way of constructing approximate load vectors is to calculate the total body and surface forces corresponding to an element and to assign equal parts to the appropriate element nodal degrees of freedom. Consider as an example the rectangular plane stress element in Fig. 4.6 with the variation of the body force shown. The total body force is equal to 2.0, and hence we obtain the lumped body force vector given in the figure. In considering the derivation of an element mass matrix, we recall that the inertia forces have been considered part of the body forces. Hence we may also establish an approximate mass matrix by lumping equal parts of the total element mass to the nodal points. Realizing that each nodal mass essentially corresponds to the mass of an element contributing volume around the node, we note that using this procedure of lumping mass, we assume in essence that the accelerations of the contributing volume to a node are constant and equal to the nodal values. An important advantage of using a lumped mass matrix is that the matrix is diagonal, and, as will be seen later, the numerical operations for the solution of the dynamic equations of equilibrium are in some cases reduced very significantly. EXAMPLE 4.21: Evaluate the lumped body force vector and the lumped mass matrix of the element assemblage in Fig. E4.5. The lumped mass matrix is 100
M=pf
80
%00
⋯∘≣∘≀⊅∊⊹∣⊃∫
or
M=§
(1+i>
150 0 0
0 670 0
0
& 0]dx
00 }
40
0
00
0
0
2000
0 0 520
Similarly, the lumped body force vector is
R; = ⋃⋮⊖⊖≺↕≻ i
=3
g
1) dx + J:w (1 ⊹⊐⋮−⊖⋗≳
150
g (11—0> dx)fz(t) 2
25022 £
It may be noted that, as required, the sums of the elements in M and Rz in both this
example and in Example 4.5 are the same.
214
Formulation of the Finite Element Method
Chap. 4
When using the load lumping procedure it should be recognized that the nodal point loads are, in general, calculated only approximately, and if a coarse finite element mesh is employed, the resulting solution may be very inaccurate. Indeed, in some cases when higherorder finite elements are used, surprising results are obtained. Figure 4.7 demonstrates such a case (see also Example 5.12).
Thickness = 1 cm
y
p = 300 N/em? E =3 x 107 Njem? v=03
Integration
Nio
B
Nlo
{b) Finite element mode! with consistent loading
o
ini
Figure 4.7
odel
e ot
' 2 Nfem*) of units have (All stresses
(3 x 3 Gauss points are used, see Table 5.7)
Some sample analysis results with and without consistent loading
Considering dynamic analysis, the inertia effects can be thought of as body forces. Therefore, if a lumped mass matrix is employed, little might be gained by using a consistent load vector, whereas consistent nodal point loads should be used if a consistent mass matrix is employed in the analysis. 4.2.5
Exercises
4.1. Use the procedure in Example 4.2 to formally derive the principle of virtual work for the onedimensional bar shown.
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method Y,
215
A(x)
E = Young's modulus
The differential equations of equilibrium are Ei(A%>+ ax ax
B =0
a
EAZ ax
=R x=L
4.2. Consider the structure shown. (a) Write down the principle of virtual displacements by specializing the general equation (4.7) to this case. (b) Use the principle of virtual work to check whether the exact solution is
7% Use the following (iii) u(x) = apx>. (¢)
three
virtual
72 ={=4
(73
F 24x\ —]—
73L> e
displacements:
(i) u(x) = aox,
(ii) w(x) = aox?,
Solve the governing differential equations of equilibrium, a du E—[A—) = 8x( E)x) 0
EA ⋮↾⊻ ax
=F x=L
(d) Use the three different virtual displacement patterns given in part (b), substitute into the
principle of virtual work using the exact solution for the stress [from part (c)], and explicitly show that the principle holds.
Z
ol
] F = total force exerted on right end E = Young's modulus
Al(x) = Apl1  x/4L)
Formulation of the Finite Element Method
216
Chap. 4
4.3. Consider the bar shown. (a) Solve for the exact displacement response of the structure. (b) Show explicitly that the principle of virtual work is satisfied with the displacement patterns (i) # = ax and (i) ¥ = ax> (c) Identify a stress T, for which the principle of virtual work is satisfied with pattern (ii) but not with pattern (i).
A= Ag(43x/L)
f? = constant force per unit length Young's modulus £
4.4. For the twodimensional body shown, use the principle of virtual work to show that the body forces are in equilibrium with the applied concentrated nodal loads.
2 =101 + 2x) N/m? 2 =201 + y) N/m® Ry =60 Rz
N
=45N
R; =
I5N
Unit thickness 2m
B fy 8
1m
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
217
4.5, Idealize the bar structure shown as an assemblage of 2 twonode bar elements.
(a) Calculate the equilibrium equations KU = R, (b) Calculate the mass matrix of the element assemblage.
e 60
80
w—
20 Aol
f8(x) = 0.1f, force/unit volume E = Young's modulus
f'lx
4.6. Consider the disk with a centerline hole of radius 20 shown spinning at a rotational velocity of o radians/second.
20
60

80
E =Young's modulus p = mass density v = Poisson’s ratio
Idealize the structure as an assemblage of 2 twonode elements and calculate the steadystate
(pseudostatic) equilibrium equations. (Note that the strains are now du/dx and u/x, where u/x is the hoop strain.) 4.7. Consider Example 4.5 and the state at time ¢ = 2.0 with U,(f) = O at all times. (a) Use the finite element formulation given in the example to calculate the static nodal point
displacements and the element stresses.
(b) Calculate the reaction at the support.
Formulation of the Finite Element Method
218 (c)
Chap. 4
Let the calculated finite element solution be «FE. Calculate and plot the error r measured in
satisfying the differential equation of equilibrium, i.e.,
o[
r
= E[—A
ou F")]
[ax (
ox
+
!
f2A
(d) Calculate the strain energy of the structure as evaluated in the finite element solution and compare this strain energy with the exact strain energy of the mathematical model.
4.8. The twonode truss element shown, originally at a uniform temperature, 20°C, is subjected to a temperature variation 0 = (10x + 20)°C Calculate the resulting stress and nodal point displacement. Also obtain the analytical solution, assuming a continuum, and briefly discuss your results.
'
E = 200,000 .
A=1
a=1x10" (per °C)
4.9. Consider the finite element analysis illustrated.
Young's modulus E Poisson's rstio v = 0.30
Plane stress condition (thickness t). All elements sre 4node elements
(a) Begin [141
by establishing Uy
U
U2
Uz
the typical U3
Us
matrix
B
of an element
for the vector 47 =
04]
(b) Calculate the elements of the K matrix, Kuy,u,, Kusu;, Ky, v,
and Kygy,, of the structural
assemblage. (c) Calculate the nodal load R, due to the linearly varying surface pressure distribution,
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
219
U 2
1
3in v3
u3 3 4in ——»‘
4.10. Consider the simply supported beam shown. (a) Assume that usual beam theory is employed and use the principle of virtual work to evaluate the reactions R, and R;.
(b) Now assume that the beam is modeled by a fournode finite element. Show that to be able
to evaluate R, and R, as in part (a) it is necessary that the finite element displacement functions can represent the rigid body mode displacements.

5
 R1
Ry
4.11. The fournode plane stress element shown carries the initial stresses
7. = 0 MPa
7, = 10 MPa 75, = 20 MPa
220
Formulation of the Finite Element Method
Chap. 4
(a) Calculate the corresponding nodal point forces R;. (b) Evaluate the nodal point forces Rs equivalent to the surface tractions that correspond to the element stresses. Check your results using elementary statics and show that Ry is equal to R, (c)
evaluated in part (a). Explain why this result makes sense. Derive a general result: Assume that any stresses are given, and R; and Rg are calculated.
What conditions must the given stresses satisfy in order that R, = Ry, where the surface tractions in Ry are obtained from equation (b) in Example 4.2?
e
T
⋁
1
30 mm
Young's modulus £ Poisson's rstio v Thickness = 0.5 mm
4.12. The fournode plane strain element shown is subjected to the constant stresses Tex = 20 psi
Ty = 10 psi Ty = 10 psi Calculate the nodal point displacements of the element. 

↿⊺ 2in
y
3
∆≀⊥
x
Young's modulus E = 30 x 10° psi Poisson's ratio v = 0.30
221
Formulation of the DisplacementBased Finite Element Method
Sec. 4.2
4.13. Consider element 2 in Fig. E4.9. (a) Show explicitly that
F® = ∫
v
B(z)r.,.(z)dvtz)
(b) Show that the element nodal point forces F®@ are in equilibrium. 4.14. Assume that the element stiffness matrices K, and Kj corresponding to the element displacements shown have been calculated. Assemble these element matrices directly into the global structure stiffness matrix with the displacement boundary conditions shown.
Individual elements
usz
Structural assemblage and degrees of freedom
KA
=
an
4n
Gz
Gu
ais
Qe  W
bu
bz
bz
bu
bis
bie
w
az
Gn
Gp
Gu
G5
G
U1
by
bxn
by
by
bis
b
vr
a3
a3
A3
G
G
Q3
Uz
by,
bn
by
by
bis
b 
6
Ay
Q4
Q43
Qa4
Qg5
Qas 
U2
by
by
bss
bas
biss
bas 
w2
QGss
Gse  Us
bsy
bsy
bss
bsa
bss
bss  vz
Qg5
Qe  U3
bey
bex
bes
bss
bes
bes  &
Gsy
Gs;
Gs3
Qs4
As1
Q2
Qe3
QGss
KB
_=
4.15. Assume that the element stiffness matrices K, and K corresponding to the element displacements shown have been calculated. Assemble these element matrices directly into the global structure stiffness matrix with the displacement boundary conditions shown.
Ka=
U
222
Formulation of the Finite Element Method 7]
Chap. 4
w1
6, U
vi
un
u
A
vs us
B
va
6, Uy
v:
* u
2
4,16. Consider Example 4.11. Assume that at the support A, the roller allows a displacement only along a slope of 30 degrees to the horizontal direction. Determine the modifications necessary in the solution in Example 4.11 to obtain the structure matrix K for this situation.
(a) Consider imposing the zero displacement condition exactly. (b) Consider imposing the zero displacement condition using the penalty method. Quadrilateral plane stress element
4.17. Consider the beam element shown. Evaluate the stiffness coefficients ky; and k2. (a) Obtain the exact coefficients from the solution of the differential equation of equilibrium
(using the mathematical model of Bernoulli beam theory). (b) Obtain the coefficients using the principle of virtual work with the Hermitian beam functions (see Example 4.16). hix) = hg(1+ x/L)
Young's modulus £ Unit thickness
Sec. 4.2
Formulation of the DisplacementBased Finite Element Method
223
4.18. Consider the twoelement assemblage shown. (a) Bvaluate the stiffness coefficients X,;, K4 for the finite element idealization. (b) Calculate the load vector of the element assemblage.
Element 1
Element 2
E = 200,000 v=03
Plane stress, thickness = 0.1
4.19. Consider the twoelement assemblage in Exercise 4.18 but now assume axisymmetric conditions. The yaxis is the axis of revolution. (a) Evaluate the stiffness coefficients X,;, K4 for the finite element idealization. (b) Evaluate the corresponding load vector.
4.20. Consider Example 4.20 and let the loading on the structure be R, = f,(¢) cos 6. (a) Establish the stiffness matrix, mass matrix, and load vector of the threenode element
E = 200,000 v=03
p = mass density
Formulation of the Finite Element Method
224
Chap. 4
shown. Establish explicitly all matrices you need but do not perform any multiplications and integrations. (b) Explain (by physical reasoning) that your assumptions on «, v, w make sense. 4.21. An inviscid fluid element (for acoustic motions) can be obtained by considering only volumetric strain energy (since inviscid fluids provide no resistance to shear). Formulate the finite element
fluid stiffness matrix for the fournode plane element shown and write out all matrices required. Do not actually perform any integrations or matrix multiplications. Hint: Remember that
p=—BAV/Vandt" =[r,.
7,
Ty
T
=[—p
0
—p
b‘ l
—plandAV/V = €. + ¢,
Thickness t Bulk modulus8
]
w o
4.22, Consider the element assemblages in Exercises 4.18 and 4.19. For each case, evaluate a lumped mass matrix (using a uniform mass density p) and a lumped load vector. 4.23. Use a finite element program to solve the model shown of the problem in Example 4.6.
(a) Print out the element stresses and element nodal point forces and draw the “exploded element views” for the stresses and nodal point forces as in Example 4.9, (b) Show that the element nodal point forces of element 5 are in equilibrium and that the element nodal point forces of elements 5 and 6 equilibrate the applied load. (c) Print out the reactions and show that the element nodal point forces equilibrate these reactions. (d) Calculate the strain energy of the finite element model. P=100
Eight constantstrsin trisngles
Sec. 4.3
Convergence of Analysis Results
225
4.24. Use a finite element program to solve the model shown of the problem in Example 4.6. Print out the element stresses and reactions and calculate the strain energy of the model. Draw the “exploded element views” for the stresses and nodal point forces. Compare your results with those for Exercise 4.23 and discuss why we should not be surprised to have obtained different results (although the same kind and same number of elements are used in both idealizations). P=100
Eight constantstrain triangles
43 CONVERGENCE
OF ANALYSIS RESULTS
Since the finite element method is a numerical procedure for solving complex engineering problems, important considerations pertain to the accuracy of the analysis results and the convergence of the numerical solution. The objective in this section is to address these issues. We start by defining in Section 4.3.1 what we mean by convergence. Then we consider in a rather physical manner the criteria for monotonic convergence and relate these criteria to the conditions in a Ritz analysis (introduced in Section 3.3.3). Next, some important properties of the finite element solution are summarized (and proven) and the rate of convergence is discussed. Finally, we consider the calculation of stresses and the evaluation of error measures that indicate the magnitude of the error in stresses at the completion of an analysis. We consider in this section displacementbased finite elements leading to monotonically convergent solutions. Formulations that lead to a nonmonotonic convergence are considered in Sections 4.4 and 4.5. 4.3.1 The Model Problem and a Definition of Convergence
Based on the preceding discussions, we can now say that, in general, a finite element analysis requires the idealization of an actual physical problem into a mathematical model and then the finite element solution of that model (see Section 1.2). Figure 4.8 summarizes these concepts. The distinction given in the figure is frequently not recognized in practical analysis because the differential equations of motion of the mathematical model are not dealt with, and indeed the equations may be unknown in the analysis of a complex problem,
226
Formulation of the Finite Element Method
Chap. 4
Actual Physical Problem Geometric domain Material Loading Boundary conditions

Mathematical Model (which corresponds to a mechanical idealization) Kinematics, e.g., truss plane stress threedimensional irch
Yields:
Kirchhoff plate
_
Material, e.g.,
_
Governing differential
__
equation(s) of motion
isotropic linear
e.g.,
elastic MooneyRivlin rubber
etc.
ai (EA ‘7_“) =  plx)
Loading, e.g.,
concentrated centrifugal etc.
x 9x and principle of virtual work equation
Boundary
prescribed
(see Example 4.2)
Conditions, e.g., displecements ↓
etc.
Yields:
Finite Element Solution Choice of elements and solution procedures
Figure 4.8
_
Approximate solution of the mathematical model (that is, approximate response of mechanical idealization)
Finite element solution process
such as the response prediction of a threedimensional shell. Instead, in a practical analysis, a finite element idealization of the physical problem is established directly. However, to study the convergence of the finite element solution as the number of elements increases, itis valuable to recognize that a mathematical model is actually implied in the finite element representation of the physical problem. That is, a proper finite element solution should converge (as the number of elements is increased) to the analytical (exact) solution of the differential equations that govern the response of the mathematical model. Furthermore, the convergence behavior displays all the characteristics of the finite element scheme because the differential equations of motion of the mathematical model express in a very precise and compact manner all basic conditions that the solution variables (stress, displacement, strain, and so on) must satisfy. If the differential equations of motion are not known, as in a complex shell analysis, and/or analytical solutions cannot be obtained, the convergence of the finite element solutions can be measured only on the fact that all basic kinematic, static, and constitutive conditions contained in the mathematical model must ultimately (at convergence) be satisfied. Therefore, in all discussions of the convergence of
finite element solutions we imply that the convergence to the exact solution of a mathematical model is meant. Here it is important to recognize that in linear elastic analysis there is a unique exact solution to the mathematical model. Hence if we have a solution that satisfies the governing
Sec. 4.3
Convergence of Analysis Results
227
mathematical equations exactly, then this is the exact solution to the problem (see Section 4.3.4). In considering the approximate finite element solution to the exact response of the mathematical model, we need to recognize that different sources of errors affect the finite
element solution results. Table 4.4 summarizes various general sources of errors. Round off errors are a result of the finite precision arithmetic of the computer used; solution errors in the constitutive modeling are due to the linearization and integration of the constitutive relations; solution errors in the calculation of the dynamic response arise in the numerical integration of the equations of motion or because only a few modes are used in a mode superposition analysis; and solution errors arise when an iterative solution is obtained because convergence is measured on increments in the solution variables that are small but not zero. In this section, we will discuss only the finite element discretization errors, which are due to interpolation of the solution variables. Thus, in essence, we consider in this section a model problem in which the other solution errors referred to above do not arise: a linear elastic static problem with the geometry represented exactly with the exact calculation of the element matrices and solution of equations, i.e., also negligible roundoff errors. For ease of presentation, we assume that the prescribed displacements are zero. Nonzero displacement boundary conditions would be imposed as discussed in Section 4.2.2, and such boundary conditions do not change the properties of the finite element solution. For this model problem, let us restate for purposes of our discussion the basic equation of the principle of virtual work governing the exact solution of the mathematical model
∫ ∊⊺⇅↿⋁∶ ∫ ﬁsf’fsfds+ ∫ ﬁTf”dV v
TABLE
4.4
Sf
(4.62)
v
Finite element solution errors
Error
Error occurrence in
See section
Discretization
Use of finite element interpolations for geometry and solution variables Evaluation of finite element matrices using numerical integration
4.2.1 423,53
Use of nonlinear material
6.6.3
Numerical integration in space Evaluation of constitutive,
models
5.5 6.8.4
6.6.4
relations” Solution of dynamic equilibrium equations Solution of finite element equations by iteration Roundoff
Direct time integration, mode superposition
9.29.4
GaussSeidel, conjugate gradient, NewtonRaphson, quasiNewton methods, eigensolutions Setting up equations and
8.3,84 9.5 10.4
their solution
8.2.6
228
Formulation of the Finite Element Method
Chap. 4
We recall that for T to be the exact solution of the mathematical model, (4.62) must
hold for arbitrary virtual displacements U (and corresponding virtual strains € ), with Wizero at and corresponding to the prescribed displacements. A short notation for (4.62) 1s Find the displacements u (and corresponding stresses 7) such that a(u, v) = (f, v)
for all admissible v
(4.63)
Here a(,") is a bilinear form and (f,) is a linear form®>—these forms depend on the mathematical model considered—u is the exact displacement solution, v is any admissible virtual displacement [“admissible” because the functions v must be continuous and zero at and corresponding to actually prescribed displacements (see (4.7)], and f represents the forcing functions (loads % and £®). Note that the notation in (4.63) implies an integration process. The bilinear forms a(,) that we consider in this section are symmetric in the sense that a(u, v) = a(v, u).
From (4.63) we have that the strain energy corresponding to the exact solution u is 1/2 a(u, u). We assume that the material properties and boundary conditions of our model problem are such that this strain energy is finite. This is not a serious restriction in practice but requires the proper choice of a mathematical model. In particular, the material properties must be physically realistic and the load distributions (externally applied or due to displacement constraints) must be sufficiently smooth. We have discussed the need of modeling the applied loads properly already in Section 1.2 and will comment further on it in Section 4.3.4. Assume that the finite element solution is u,: this solution lies of course in the finite
element space given by the displacement interpolation functions (k denoting here the size of the generic element and hence denoting a specific mesh). Then we define “convergence” to mean that au—w,uu)—0
ash—0
(4.64)
or, equivalently [see (4.90)], that a(uy,, u,) — a(u, u)
ash—0
Physically, this statement means that the strain energy calculated by the finite element solution converges to the exact strain energy of the mathematical model as the finite element mesh is refined. Let us consider a simple example to show what we mean by the bilinear form a(,). ® The bilinearity of af:,) refers to the fact that for any constants y, and 7,
a(yiu; + v, v) = ya(uy, v) + yauy, v) a(u, v,
+ v.v2) = ya(u, vi) + pau, vo)
and the linearity of (f,.) refers to the fact rthat for any constants y, and y,,
& v +7%v) = n€ v,) + 7
v).
Sec. 4.3
229
Convergence of Analysis Results
EXAMPLE 4.22:
Assume
that a simply
supported
prestressed
membrane,
with
(constant)
prestress tension 7, subjected to transverse loading p is to be analyzed (see Fig. E4.22). Establish for this problem the form (4.63) of the principle of virtual work.
Hinged on all edges
Tension
T,
Figure E4.22
Prestressed membrane
The principle of virtual work gives for this problem
W)™ [ ax
Lff
ay
T
ax
dxdy = prxdy w d. f_wxy
ay
where w(x, y) is the transverse displacement. The lefthand side of this equation gives the bilinear form a(v, ), with v = W, ¥ = w, and the integration on the righthand side gives (f, v).
Depending on the specific (properly formulated) displacementbased finite elements used in the analysis of the model problem defined above, we may converge monotonically or nonmonotonically to the exact solution as the number of finite elements is increased. In the following discussion we consider the criteria for the monotonic convergence of solutions. Finite element analysis conditions that lead to nonmonotonic convergence are discussed in Section 4.4. 4.3.2 Criteria for Monotonic Convergence
For monotonic convergence, the elements must be complete and the elements and mesh must be compatible. If these conditions are fulfilled, the accuracy of the solution results will increase continuously as we continue to refine the finite element mesh. This mesh re
finement should be performed by subdividing a previously used element into two or more elements; thus, the old mesh will be “embedded” in the new mesh. This means mathemat
ically that the new space of finite element interpolation functions will contain the previously used space, and as the mesh is refined, the dimension of the finite element solution space will be continuously increased to contain ultimately the exact solution. The requirement of completeness of an element means that the displacement functions of the element must be able to represent the rigid body displacements and the constant strain states.
Formulation of the Finite Element Method
230
Chap. 4
The rigid body displacements are those displacement modes that the element must be able to undergo as a rigid body without stresses being developed in it. As an example, a twodimensional plane stress element must be able to translate uniformly in either direction
3

of its plane and to rotate without straining. The reason that the element must be able to undergo these displacements without developing stresses is illustrated in the analysis of the cantilever shown in Fig. 4.9: the element at the tip of the beam—for any element size— must translate and rotate stressfree because by simple statics the cantilever is not subjected to stresses beyond the point of load application. The number of rigid body modes that an element must be able to undergo can usually be identified without difficulty by inspection, but it is instructive to note that the number of element rigid body modes is equal to the number of element degrees of freedom minus the number of element straining modes (or natural modes). As an example, a twonoded truss has one straining mode (constant strain state), and thus one, three, and five rigid body modes in one, two, and threedimensional conditions, respectively. For more complex finite
(a) Rigid body modes of a plane stress element
Distributed load p
Rigid body translation and rotation; element must be
Y,
stressfree for any element size
(b) Analysis to itlustrate the rigid body mode conditon
Figure 4.9
Use of plane stress element in analysis of cantilever
Sec. 4.3
Convergence of Analysis Results
231
elements the individual straining modes and rigid body modes are displayed effectively by representing the stiffness matrix in the basis of eigenvectors. Thus, solving the eigenproblem
Kb = rd
(4.65)
K® = ®A
(4.66)
we have (see Section 2.5)
where ® is a matrix storing the eigenvectors &, . . ., ¢, and A is a diagonal matrix storing the corresponding eigenvalues, A = diag(A;). Using the eigenvector orthonormality property we thus have
KD
= A
(4.67)
We may look at A as being the stiffness matrix of the element corresponding to the eigenvector displacement modes. The stiffness coefficients A,. . ., A, display directly how stiff the element is in the corresponding displacement mode. Thus, the transformation in (4.67) shows clearly whether the rigid body modes and what additional straining modes are present.'® As an example, the eight eigenvectors and corresponding eigenvalues of a fournode element are shown in Fig. 4.10. The necessity for the constant strain states can be physically understood if we imagine that more and more elements are used in the assemblage to represent the structure. Then in the limit as each element approaches a very small size, the strain in each element approaches a constant value, and any complex variation of strain within the structure can be approximated. As an example, the plane stress element used in Fig. 4.9 must be able to represent two constant normal stress conditions and one constant shearing stress condition. Figure 4.10 shows that the element can represent these constant stress conditions and, in
addition, contains two flexural straining modes. The rigid body modes and constant strain states that an element can represent can also be directly identified by studying the element straindisplacement matrix (see Example 4.23). The requirement of compatibility means that the displacements within the elements and across the element boundaries must be continuous. Physically, compatibility ensures that no gaps occur between elements when the assemblage is loaded. When only translational degrees of freedom are defined at the element nodes, only continuity in the displacements u, v, or w, whichever are applicable, must be preserved. However, when rotational degrees of freedom are also defined that are obtained by differentiation of the transverse displacement (such as in the formulation of the plate bending element in Example 4.18), it is also necessary to satisfy element continuity in the corresponding first displacement derivatives. This is a consequence of the kinematic assumption on the displacements over the depth of the plate bending element; that is, the continuity in the displacement w and the
derivatives aw/0x and/or ow/dy along the respective element edges ensures continuity of displacements over the thickness of adjoining elements. Compatibility is automatically ensured between truss and beam elements because they join only at the nodal points, and compatibility is relatively easy to maintain in '0Note also that since the finite element analysis overestimates the stiffness, as discussed in Section 4.3.4, the “smaller” the eigenvalues, the more effective the element.
Formulation of the Finite Element Method
232
Thickness = 1.0
1.0
l'*—"l −↑↽
Young's −−−−∣
Chap. 4
———————————
1
!
i
i
modulus = 1.0 1
⋮ :
1.0 i
−−−−⇃
Poisson's ratio = 0.30
I :
Rigid body mode 41=0
Rigid body mode ;=0
Stretching mode A; = 0.769
Uniform extension mode Ag = 1.43
Figure 410
Eigenvalues and eigenvectors of fournode plane stress element
twodimensional plane strain, plane stress, and axisymmetric analysis and in threedimensional analysis, when only u, v, and w degrees of freedom are used as nodal point variables. However, the requirements of compatibility are difficult to satisfy in plate bending analysis, and particularly in thin shell analysis if the rotations are derived from the transverse displacements. For this reason, much emphasis has been directed toward the
development of plate and shell elements, in which the displacements and rotations are
Sec. 4.3
Convergence of Analysis Results
233
variables (see Section 5.4). With such elements the compatibility requirements are just as easy to fulfill as in the case of dealing only with translational degrees of freedom.
Whether a specific element is complete and compatible depends on the formulation used, and each formulation need be analyzed individually. Consider the following simple example.
EXAMPLE 4.23: complete.
Investigate if the plane stress element used in Example 4.6 is compatible and
We have for the displacements of the element,
u(x,y) = ay + aox + azy + axy
o(x,y) = B, + Bax + Bsy + Bixy
Observing that the displacements within an element are continuous, in order to show that the element is compatible, we need only investigate if interelement continuity is also preserved when an element assemblage is loaded. Consider two elements interconnected at two node points (Fig. E4.23) on which we impose two arbitrary displacements. It follows from the displacement assumptions that the points (i.e., the material particles) on the adjoining element edges displace linearly, and therefore continuity between the elements is preserved. Hence the element is compatible.
3
Particles on element edges remelin together
ode
x u Figure E4.23
Compatibility of plane stress element
Considering completeness, the displacement functions show that a rigid body translation in the x direction is achieved if only a, is nonzero. Similarly, a rigid body displacement in the y direction is imposed by having only 3, nonzero, and for a rigid body rotation a; and 8, must be nonzero only with 8, = —as. The same conclusion can also be arrived at using the matrix E that relates the strains to the generalized coordinates (see Example 4.6). This matrix also
shows that the constant strain states are possible. Therefore the element is complete.
234
Formulation of the Finite Element Method
Chap. 4
4.3.3 The Monotonically Convergent Finite Element Solution: A Ritz Solution
We observed earlier that the application of the principle of virtual work is identical to using the stationarity condition of the total potential of the system (see Example 4.4). Considering also the discussion of the Ritz method in Section 3.3.3, we can conclude that monotonically convergent displacementbased finite element solutions are really only applications of this method. In the finite element analysis the Ritz functions are contained in the element displacement interpolation matrices H™, m = 1, 2,. . ., and the Ritz parameters are the unknown nodal point displacements stored in U. As we discuss further below, the mathe
matical conditions on the displacement interpolation functions in the matrices H™, in order that the finite element solution be a Ritz analysis, are exactly those that we identified earlier using physical reasoning. The correspondence between the analysis methods is illustrated in Examples 3.22 and 4.5. Considering the Ritz method of analysis with the finite element interpolations, we
have
I1 = LU'KU — U'R
(4.68)
where II is the total potential of the system. Invoking the stationarity of IT with respect to the Ritz parameters U; stored in U and recognizing that the matrix K is symmetric, we obtain (4.69)
KU=R
The solution of (4.69) yields the Ritz parameters, and then the displacement solution in the domain considered is
u™ = HmU:
m=12,...
(4.70)
The relations in (4.68) to (4.70) represent a Ritz analysis provided the functions used satisfy certain conditions. We defined in Section 3.3.2 a C™! variational problem as one in which the variational indicator of the problem contains derivatives of order m and lower. We then noted that for convergence the Ritz functions must satisfy the essential (or geometric) boundary conditions of the problem 1nv?olv1ng derivatives up to order m — 1), but that the functions do not need to satisfy the natural (or force) boundary conditions involving derivatives of order m to (2m — 1) because these conditions are implicitly contained in the variational indicator I1. Therefore, in order for a finite element solution to be a Ritz analysis, the essential boundary conditions must be completely satisfied by the finite element nodal point displacements and the displacement interpolations between the nodal points. However, in selecting the finite element displacement functions, no special attention need be given to the natural boundary conditions because these conditions are imposed with the load vector and are satisfied approximately in the Ritz solution. The accuracy with which the natural or force boundary conditions are satisfied depends on the specific Ritz functions employed, but this accuracy can always be increased by using a larger number of functions, i.e., a larger number of finite elements to model the problem. In the classical Ritz analysis the Ritz functions extend over the complete domain considered, whereas in the finite element analysis the individual Ritz functions extend only over subdomains (finite elements) of the complete region. Hence, there must be a question
as to what conditions must be fulfilled by the finite element interpolations with regard to
Sec. 4.3
Convergence of Analysis Results
235
continuity requirements between adjacent subdomains. To answer this question we consider the integrations that must be performed to evaluate the coefficient matrix K. We recognize that in considering a C™~! problem we need continuity in at least the (m — 1)st derivatives of the Ritz trial functions in order that we can perform the integrations across the element boundaries. However, this continuity requirement corresponds entirely to the element compatibility conditions that we discussed in Section 4.3.2. For example, in the analysis of fully threedimensional problems only the displacements between elements must be continuous, whereas in the analysis of plate problems formulated using the Kirchhoff plate theory we also need continuity in the first derivatives of the displacement functions.
In summary, therefore, for a C™~! problem [C™™! = continuity on trial functions and their derivatives up to order (m — 1)], in the classical Ritz analysis the trial functions are .selected to satisfy exactly all boundary conditions that involve derivatives up to order
(m — 1). The same holds in finite element analysis, but in addition, continuity in the trial functions and their derivatives up to order (m — 1) must be satisfied between elements in
order for the finite element solution to correspond to a Ritz analysis. Although the classical Ritz analysis procedure and the displacementbased finite element method are theoretically identical, in practice, the finite element method has important advantages over a conventional Ritz analysis. One disadvantage of the conventional Ritz analysis is that the Ritz functions are defined over the whole region considered. For example, in the analysis of the cantilever in Example 3.24, the Ritz functions spanned from x = 0 to x = L. Therefore, in the conventional Ritz analysis, the matrix K is a full
matrix, and as pointed out in Section 8.2.3, the numerical operations required for solution of the resulting algebraic equations are considerable if many functions are used. A particular difficulty in a conventional Ritz analysis is the selection of appropriate Ritz functions since the solution is a linear combination of these functions. In order to solve accurately for large displacement or stress gradients, many functions may be needed. However, these functions also unnecessarily extend over the regions in which the displacements and stresses vary rather slowly and where not many functions are needed. Another difficulty arises in the conventional Ritz analysis when the total region of interest is made up of subregions with different kinds of strain distributions. As an example, consider a plate that is supported by edge beams and columns. In such a case, the Ritz functions used for one region (e.g., the plate) are not appropriate for the other regions (i.e., the edge beams and columns), and special displacement continuity conditions or boundary relations must be introduced. ~ The few reasons given already show that the conventional Ritz analysis is, in general, not particularly computeroriented, except in some cases for the development of specialpurpose programs. On the other hand, the finite element method has to a large extent removed the practical difficulties while retaining the advantageous properties of the conventional Ritz method. With regard to the difficulties mentioned above, the selection of Ritz functions is handled by using an adequate element library in the computer program. The use of relatively many functions in regions of high stress and displacement gradients is possible simply by using many elements, and the combination of domains with different kinds of strain distributions is possible by using different kinds of elements to idealize the domains. It is this generality of the finite element method, and the good mathematical foundation, that have made the finite element method the very widely used analysis tool in today’s engineering environments.
Formulation of the Finite Element Method
236
Chap. 4
4.3.4 Properties of the Finite Element Solution
Let us consider the general linear elasticity problem and its finite element solution and identify certain properties that are useful for an understanding of the finite element method. We shall use the notation summarized in Table 4.5. The elasticity problem can be written as follows (see, for example, G. Strang and G. F. Fix [A], P. G. Ciarlet [A], or F. Brezzi and M. Fortin [A]). Find u € V such that
awv)=@®v
Vvev
“71)
where the space V is defined as V=
{v v € L3(Vol); g—'; € L Vo), i,j = 1,2, 3,015, = 0,i = 1,2, 3}
(4.72)
o
Here L*(Vol) is the space of square integrable functions in the volume, “Vol”, of the body being considered, L*(Vol) = {wl w is defined in Vol and ∫ (i (w,)2> avol =  wiiZzey < +°°} Vol
TABLE 4.5 solutions
(4.73)
\i=1
Notation used in discussion of the properties and convergence of finite element
Symbol
Meaning
a(.,.)
Bilinear form corresponding to model problem being considered (see Example 4.22)
f
Load vector Exact displacement solution to mathematical model; an element of the space V
v Vi v e
V, Vu Vol L?
Displacements; an element of the space V
Finite element solution, an element of the space Vj, Finite element displacements; an element of the space V, For all An element of
Spaces of functions [see (4.72) and (4.84)] Volume of body considered Space of a square integrable functions [see (4.73)]
€,
Error between exact and finite element solution, e, = u — u,
3 C
There exists Contained in
⊊
Contained in but not equal to
Il e
Energy norm [see (4.74)]
inf
We take the infimum.
sup
We take the supremum.
Sec. 4.3
Convergence of Analysis Results
237
Hence, (4.72) defines a space of functions corresponding to a general threedimensional analysis. The functions in the space vanish on the boundary ., and the squares of the functions and of their first derivatives are integrable. Corresponding to V, we use the energy norm  vl = a(v, +)
(4.74)
which actually corresponds to twice the strain energy stored in the body when the body is subjected tothe displacement field v. We assume in our discussion that the structure considered in (4.71) is properly
supported, corrresponding to the zero displacement conditions on S,,, so that  v% is greater than zero for any v different from zero. In addition, we shall also use the Sobolev norms of orderm = Oand m = 1 defined as m
0:
Ul vlo? = ∫ ⋅≼∑≺∅↙≻≄⋟≴↿⋁∘∣
(4.75)
3
i=1
vl = (vl + ∫
(4.76) Vol
For our elasticity problem the norm of order 1 is used,'' and we have the following two important properties for our bilinear form a. Continuity: IM>0suchthat
Vv,
v, €V,
a(v, Vz)l
=M
∥ Vi ∥∎ ∥ ∇≖∥∣
(477)
Ellipticity:
Ja>0suchthat Vv EV,
a(v,v) = a v}
(4.78)
where the constants o and M depend on the actual elasticity problem being considered, including the material constants used, but are independent of v. "'In our discussion, we shall also use the PoincaréFriedrichs inequality, namely, that for the analysis problems we consider, for any v we have
) 5 ( 3 ( e = o j [ (B er
where c is a constant (see, for example, P. G. Ciarlet [A]).
Formulation of the Finite Element Method
238
Chap. 4
The continuity property is satisfied because reasonable norms are used in (4.77), and the ellipticity property is satisfied because a properly supported (i.e., stable) structure is
being considered (see P. G. Ciarlet [A] for a mathematical proof). Based on these properties we have
all vl = (alv, v)? =  vl
(4.79)
where ¢ and ¢, are constants independent of v, and we therefore have that the energy norm is equivalent to the 1norm (see Section 2.7). In mathematical analysis the Sobolev norms are commonly used to measure rates of convergence (see Section 4.3.5), but in practice the energy norm is frequently more easily evaluated [see (4.97)]. Because of (4.79), we can say that convergence can also be defined, instead of using (4.64), as
akh—0
la—wfi—=0
and the energy norm in problem solutions will converge with the same order as the 1norm. We examine the continuity and ellipticity of a bilinear form a in the following example. EXAMPLE 4.24: Consider the problem in Example 4.22. Show satisfies the continuity and ellipticity conditions.
that the bilinear form a
Continuity follows because'?
alw, wa) = ∫ ⊺≼≣≞≝ ⊹ ∂⊻⊻⋗ dx dy ad
A
ox
Ox
dy
9oy
=[G GG = G =G+ ] eer 
o
(A + (5 )erw ) < e s [ ]+ (2] o ol [ G+ (3o et s [+ (3o
Ellipticity requires that
However, the PoincaréFriedrichs inequality,
where ¢ is a constant, ensures that (a) is satisfied.
> Here we use the Schwarz inequality, which says that for vectorsaand b,  a  b =  a, bJ, where  
is defined in (2.148).
Sec. 4.3
Convergence of Analysis Results
239
The above statements on the elasticity problem encompass one important point already mentioned earlier: the exact solution to the problem must correspond to a finite strain energy, see (4.64) and (4.79). Therefore, for example,—strictly—we do not endeavor to solve general two or threedimensional elasticity problems with the mathematical idealization of point loads (the solution for a point load on a half space corresponds to infinite strain energy, see for instance S. Timoshenko and J. N. Goodier [A]). Instead, we represent the loads in the elasticity problem closer to how they actually act in nature, namely as smoothly distributed loads, which however can have high magnitudes and act over very small areas. Then the solution of the variational formulation in (4.71) is the same as the
solution of the differential formulation. Of course, in our finite element analysis so long as the finite elements are much larger than the area of load application, we can replace the distributed load over the area with an equivalent point load, merely for efficiency of solution; see Section 1.2 and the example in Fig. 1.4. An important observation is that the exact solution to our elasticity problem is unique. Namely, assume that u, and u, are two different solutions; then we would have
and
a(u, v) = (f, v)
Yveyv
(4.80)
a(uz, v) = (£, v)
Vvev
(4.81)
⋮
Subtracting, we obtain ⇩≼∐≖−∐∅⋅⋁⋟⇌∘
∀⋁∈∨
and takingv = u, — w,, wehave a(u, — w,,u; — u;) = 0. Using (4.79) withv =u;
(4.82)
— u,,
we obtain u, — u, , = 0, which means u, = u,, and hence we have proven that our assumption of two different solutions is untenable. Now let V; be the space of finite element displacement functions (which correspond to the displacement interpolations contained in all element displacement interpolation
matrices H™) and let v, be any element in that space (i.e., any displacement pattern that can be obtained by the displacement interpolations). Let u, be the finite element solution; hence u, is also an element in V, and the specific element that we seek. Then the finite element solution of the problem in (4.71) can be written as
Find u,
€
V, such that
(4.83)
and for the elements of this space we use the energy norm (4.74) and the Sobolev norm (4.76). Of course, V,, C V. The relation in (4.83) is the principle of virtual work for the finite element discretization corresponding to V,.. With this solution space, the continuity and ellipticity conditions (4.77) and (4.78) are satisfied, using v, € Vi, and a positive definite stiffness matrix is obtained for any V.
Formulation of the Finite Element Method
240
Chap. 4
We should note that Vi, corresponds to a given mesh, where 4 denotes the generic element size, and that in the discussion of convergence we of course consider a sequence of spaces V, (a sequence of meshes with decreasing k). We illustrate in Figure 4.11 the elements of V, for the discretization dealt with in Example 4.6.
Nodal point
Element number
S T A L Aerial view of basis functions for space V, used in analysis of cantilever plate Figure 4.11 of Example 4.6. The displacement functions are plotted upwards for ease of display, but each function shown is applicable to the u and v displacements. An element of V, is any linear combination of the 12 displacement functions. Note that the functions correspond to the element displacement interpolation matrices H™, discussed in Example 4.6, and that the displacements at nodes 1, 2 and 3 are zero.
Sec. 4.3
Convergence of Analysis Results
241
Considering the finite element solution u, and the exact solution u to the problem, we have the following important properties. Property 1. solution u, be e,
Let the error between the exact solution u and the finite element €=U
— U,
(4.85)
Then the first property is
a(e,., V;,)
=0
v
v.
eV,
(486)
The proof is obtained by realizing that the principle of virtual work gives a(u, V;,) and
=
a(u;., V},)
(f, V;,)
=
Vv,
(f, Vh)
Y
EV,
v,
(487)
eV,
(488)
so that by subtraction we obtain (4.86). We may say that the error is “orthogonal inac. , .)” to all v, in V,. Clearly, as the space V, increases, with the larger space always containing the smaller space, the solution accuracy will increase continuously. The next two properties are based on Property 1.
Property 2.
The second property is a(u,, uy) < a(u, u)
(4.89)
We prove this property by considering a(u, u) = a(u, + e, Uy + €)
= a(us, wy) + 2a(u,, e,) + ales, e)
(4.90)
= a(uy, ws) + aley, ey)
where we have used (4.86) with v, = u,. The relation (4.89) follows because a(e;, €:) > 0 for any e, # 0 (since for the properly supported structure vz > 0 for any nonzero v). Hence, the strain energy corresponding to the finite element solution is always smaller than or equal to the strain energy corresponding to the exact solution. Property 3.
The third property is
a(e;.,
e;,)
=
a(u

Vp
U
—
V;,)
Y
v,
EV,
(491)
For the proof we use that for any w, in V,, we have a(e;.
+
Wy,
€,
+
W;,)
=
a(e,.,
e;,)
+
a(w;.,
W;,)
(492)
242
Formulation of the Finite Element Method
HCHCC,
a(e,., e;.)
=
a(e,.
+
wy,
e,
+
W;,)
Chap. 4 (493)
Choosing w, = u, — v, gives (4.91). This third property says that the finite element solution u, is chosen from all the possible displacement patterns v, in V; such that the strain energy corresponding to u — 1, is the minimum. Hence, in that sense, the “energy distance” between u and the elements in Vi is minimized by the solution u, in Vj. Using (4.91) and the ellipticity and continuity of the bilinear form, we further obtain
allu  w } =a =
 v, u—u
inf a(u — vy, U — vp)
(4.94)
VAEV,
= Mv:ggh la = wiulli flu — vl
where “inf” denotes the infimum (see Table 4.5). If we let d(u, V;) = lim u — v,;, we WY recognize that we have the property la — wl = c d(u, Vi)
(4.95)
where ¢ is a constant, ¢ = V M/a, independent of h but dependent on the material properties."® This result is referred to as Cea’s lemma (see, for example, P. G. Ciarlet [A]). The above three properties give valuable insight into how the finite element solution is chosen from the displacement patterns possible within a given finite element mesh and what we can expect as the mesh is refined. We note, in particular, that (4.95), which is based on Property 3, states that a sufficient condition for convergence with our sequence of finite element spaces is that for
any u € V we have lims—¢infjju — vii = 0. Also, (4.95) can be used to measure the rate of convergence as the mesh is refined by introducing an upper bound on how d(u, V) changes with the mesh refinement (see Section 4.3.5). Further, Properties 2 and 3 say that at the finite element solution the error in strain
energy is minimized within the possible displacement patterns of a given mesh and that the strain energy corresponding to the finite element solution will approach the exact strain energy (from below) as increasingly finer meshes are used (with the displacement patterns of the finer mesh containing the displacement patterns of the previous coarser mesh). We can also relate these statements to earlier observations that in a finite element solution the stationarity of the total potential is established (see Section 4.3.2). That is, for a given mesh and any nodal displacements Uy, we have
In lU any
=
%U:nyKUany

UZnyR
(496)
13 There is a subtle point in considering the property (4.95) and the condition (4.156) discussed later; namely, while (4.95) is always valid for any values of bulk and shear moduli, the constant ¢ becomes very large as the bulk modulus increases, and the property (4.95) is no longer useful. For this reason, when the bulk modulus « is very large, we need the new property (4.156) in which the constant is independent of «, and this leads to the infsup condition.
Sec. 4.3
243
Convergence of Analysis Results
The finite element solution U is obtained by invoking the stationarity of II to obtain KU=R
At the finite element displacement solution U we have the total potential II and strain
energy U (4.97)
Therefore, to evaluate the strain energy corresponding to the finite element solution, we only need to perform a vector multiplication. To show with this notation that within the given possible finite element displacements (i.e., within the space V) Il is minimized at the finite element solution U, let us calculate
IT at U + €, where € is any arbitrary vector, IMy+e = 3(U + €)K(U
=TIy
+ €) ~ (U + €)'R
+ €"(KU — R) + 1e'Ke
(4.98)
=1I IU + %ETKE
where we used that KU
= R and the fact that K is a symmetric matrix. However, since K
is positive definite, IT y is the minimum of I1 for the given finite element mesh. As the mesh is refined, I1 will decrease and according to (4.97) U will correspondingly increase. Considering (4.89), (4.91), and (4.97), we observe that in the finite element solution the displacements are (on the “whole”) underestimated and hence the stiffness of the mathematical model is (on the “whole”) overestimated. This overestimation of the stiffness is (physically) a result of the “internal displacement constraints” that are implicitly imposed on the solution as a result of the displacement assumptions. As the finite element discretization is refined, these “internal displacement constraints” are reduced, and convergence to the exact solution (and stiffness) of the mathematical model is obtained. To exemplify the preceding discussion, Figure 4.12 shows the results obtained in the
analysis of an ad hoc test problem for twodimensional finite element discretizations. The problem is constructed so as to have no singularities. As we discuss in the next section, in this case the full (maximum) order of convergence is obtained with a given finite element in a sequence of uniform finite element meshes (in each mesh all elements are of equal square size). Figure 4.12 shows the convergence in strain energy when a sequence of uniform meshes of ninenode elements is employed for the solutions. The meshes are constructed by starting with a 2 X 2 mesh of square elements of unit side length (for which A = 1), then subdividing each element into four equal square elements (for which 2 = 3,) to obtain the second mesh, and continuing this process. We clearly see that the error in the strain energy decreases with decreasing element size h, as we would expect according to (4.91). We compare the order of convergence seen in the finite element computations with a theoretically established value in the next section.
244
Formulation of the Finite Element Method
Chap. 4
VeV 1, +1)
(+1, +1)
E = 200,000 v=0.30
(+1,1)
=1,1) N elements per side, N=2, 4, 8,...
(a) Square domain considered
u=c1(1x3 (1y? e cos kx v=2c1(1x2)(1y?) e* sin kx ¢y =constant; k=5
(b) Exact inplane displacements Obtain the finite element solution for the body loads ff and f§, where
£BX = _ (éfﬁ‘ ax
atxy)
—
oy
8 (_v_v . _Vi‘) T
ay
ot
ax
and 7y, 7y,, T, are the stresses corresponding to the exact inplane displacements given in (b).
{c) Test problem Adhoc test problem for plane stress (or plane strain, axisymmetric) elements. Figure 4.12 We use, for h small, E — E, = ¢ h* and hence log (E — E,) = log ¢ + a log h (see also (4.101)). The numerical solutions give « = 3.91.
4.3.5 Rate of Convergence
In the previous sections we considered the conditions required for monotonic convergence of the finite element analysis results and discussed how in general convergence is reached, but we did not mention the rate at which convergence occurs. As must be expected, the rate of convergence depends on the order of the polynomials used in the displacement assumptions. In this context the notion of complete polynomials is useful. Figure 4.13 shows the polynomial terms that should be included to have complete polynomials in x and y for twodimensional analysis. It is seen that all possible terms of the form x“y® are present, where @ + B = k and k is the degree to which the polynomial is complete. For example, we may note that the element investigated in Example 4.6 uses a
Sec. 4.3
Convergence of Analysis Results
245
logqo (E~ Ep)
1
2
0
logyo h
(d) Solution for plane stress problem Figure 4.12
(continued)
polynomial displacement that is complete to degree 1 only. Figure 4.13 also shows important notation for polynomial spaces. The spaces Pi correspond to the complete polynomials up to degree k. They can also be thought of as the basis functions of triangular elements: the functions in P; correspond to the functions of the linear displacement (constant strain)
triangle (see Example 4.17); the functions in P, correspond to the functions of the parabolic displacement (linear strain) triangle (see Section 5.3.2); and so on. In addition, Fig. 4.13 shows the polynomial spaces Q«, k = 1, 2, 3, which correspond
to the 4node, 9node, and 16node elements, referred to as Lagrangian elements because the displacement functions of these elements are Lagrangian functions (see also Section 5.5.1).
In considering threedimensional analysis of course a figure analogous to Fig. 4.13 could be drawn in which the variable z would be included. Let us think about a sequence of uniform meshes idealizing the complete volume of the body being considered. A mesh of a sequence of uniform meshes consists of elements of equal size—square elements when the polynomial spaces Qi are used. Hence, the
parameter A can be taken to be a typical length of an element side. The sequence is obtained by taking a starting mesh of elements and subdividing each element with a natural pattern to obtain the next mesh, and then repeating this process. We did this in solving the ad hoc test problem in Fig. 4.12. However, considering an additional analysis problem, for example, the problem in Example 4.6, we would in Fig. 4.11 subdivide each fournode element into four equal new fournode elements to obtain the first refined mesh; then we would
Chap. 4
Formulation of the Finite Element Method
246
Degree of complete polynomial
Corresponding space
............................ 1P1
........................ 2 i
P2
3........................... Ps e
P 5
 Figure 4.13
]
Polynomial terms in twodimensional analysis, Pascal triangle
subdivide each element of the first refined mesh into four equal new fournode elements to obtain the second refined mesh; and so on. The continuation of this subdivision process would give the complete sequence of meshes. To obtain an expression for the rate of convergence, we would ideally use a formula giving d(u, V;) in (4.95) as a function of h. However, such a formula is difficult to obtain, and it is more convenient to use interpolation theory and work with an upper bound on d(ll,
V;,)
Let us assume that we employ elements with complete polynomials of degree k and that the exact solution u to our elasticity problem is “smooth” in the sense that the solution
(S 23]
satisfies the relation'
3
L33 n=2
i=1
r+s+t=n
(———a"“" >2] dVol}m < oxi
dx3
(4.99)
ax'3
where of course k = 1. Therefore, we assume that all derivatives of the exact solution up to order (k + 1) in (4.99) can be calculated.
A basic result of interpolation theory is that there exists an interpolation function u; €
V, such that
lu = wl = ¢ hr*ulls
(4.100)
where h is the mesh size parameter indicating the “size” of the elements and ¢ is a constant independent of . Typically, & is taken to be the length of the side of a generic element or the diameter of a circle encompassing that element. Note that w, is not the finite element solution in V, but merely an element in V, that geometrically corresponds to a function 4We then have u is an element of the Hilbert space H**!.
Sec. 4.3
Convergence of Analysis Results
247
close to u. Frequently, as nere, we iew w,, at the hnite element nodes, take the value of the exact solution u.
Using (4.100) and Property 3 discussed in Section 4.3.4 [see (4.91)], we can now show that the rate of convergence of the finite element solution u, to the exact theory of elasticity solution u is given by the error estimate lu = wf = c r*[lufln
4.101)
where c is a constant independent of 4 but dependent on the material properties. Namely, using (4.95) and (4.100), we have ha — w.h = cd(u, Vi)
= ¢ é kM ull which gives (4.101) with a new constant c. For (4.101), we say that the rate of convergence is given by the complete righthandside expression, and we say that the order of convergence is k or, equivalently, that we have o(h*) convergence. Another way to look at the derivation of (4.101)—which is of course closely related to the previous derivation—is to use (4.79) and (4.91). Then we have
1 le = wiflh = —[a(w ~ w,, u ~ w)] 1
= cl [a(u — w;, u — u))”?
!
(4.101b)
A
c
∁−≣∥∥ −∥∣∥≖
≤∁∣⇂∣≴∥∥∥∣≺≁≖ Hence, we see directly that to obtain the rate of convergence, we really only expressed the distance d(u,V,) in terms of an upper bound given by (4.100). In practice, we frequently simply write (4.101) as
flu = wufl = ch*
(4.102)
and we now recognize that the constant ¢ used here is independent of 4 but depends on the solution and the material properties [because c in (4.101a) and c», ¢, in (4.101b) depend on the material properties]. This dependence on the material properties is detrimental when (almost) incompressible material conditions are considered because the constant then be
comes very large and the order of convergence k results in good accuracy only at very small (impractical) values of . For this reason we need in that case the property (4.95) with the constant independent of the material properties, and this requirement leads to the condition (4.156) (see Section 4.5).
248
Formulation of the Finite Element Method
Chap. 4
The constant ¢ also depends on the kind of element used. While we have assumed that the element is based on a complete polynomial of order &, different kinds of elements within that class in general display a different constant ¢ for the same analysis problem (e.g., triangular and quadrilateral elements). Hence, the actual magnitude of the error may be considerably different for a given A, while the order with which the error decreases as the mesh is refined is the same. Clearly, the magnitude of the constant ¢ can be crucial in practical analysis because it largely determines how small 4 actually has to be in order to reach an acceptable error. These derivations of course represent theoretical results, and we may question in how far these results are applicable in practice. Experience shows that the theoretical results indeed closely represent the actual convergence behavior of the finite element discretizations being considered. Indeed, to measure the order of convergence, we may simply consider the equal sign in (4.102) to obtain
log (Jlu—wlj) =logec + klogh
(4.103)
Then, if we plot from our computed results the graph of log ( u — w,;) versus log &, we find that the resulting curve indeed has the approximate slope £ when # is sufficiently small.
Evaluating the Sobolev norm may require considerable effort, and in practice, we may use the equivalence of the energy norm with the 1norm. Namely, because of (4.79), we see that (4.101) also holds for the energy norm on the left side, and this norm can frequently be evaluated more easily [see (4.97)]. Figure 4.12 shows an application. Note that the error in strain energy can be evaluated simply by subtracting the current strain energy from the strain energy of the limit solution (or, if known, the exact solution) [see (4.90)]. In the solution in Fig. 4.12 we obtained an order of convergence (of the numerical results) of 3.91, which compares very well with the theoretical value of 4 (here k = 2 and the strain energy is the energy norm squared). Further results of convergence for this ad hoc problem are given in Fig. 5.39 (where distorted elements with numerically integrated stiffness matrices are considered).
The relation in (4.101) gives, in essence, an error estimate for the displacement gradient, hence for the strains and stresses, because the primary contribution in the 1norm will be due to the error in the derivatives of the displacements. We will primarily use (4.101) and (4.102) but also note that the error in the displacements is given by
e = wifo < ¢ A*** ffuffers
(4.104)
Hence, the order of convergence for the displacements is one order higher than for the strains. These results are intuitively reasonable. Namely, let us think in terms of a Taylor series analysis. Then, since a finite element of “dimension #” with a complete displacement expansion of order k can represent displacement variations up to that order exactly, the local error in representing arbitrary displacements with a uniform mesh should be o(h**!). Also, for a C™~! problem the stresses are calculated by differentiating the displacements m times,
and therefore the error in the stresses is o(h**'™™). For the theory of elasticity problem considered above, m = 1, and hence the relations in (4.101) and (4.104) are what we might
expect.
Sec. 4.3
249
Convergence of Analysis Results
EXAMPLE 4.25: Consider the problem shown in Figure E4.25. Estimate the error of the finite element solution if linear twonode finite elements are used.
Constant crosssectional area A Young's modulus £
{s) Bsr subjected to losd per unit length F3(x) = ax
{b) Solutions (for finite element solution three elements sre used)
Figure E4.25
Analysis of bar
The finite element problem in this case is to calculate u, € Vj, such that
Yor EV
(EA uh, 0h) = (f?, va)
with
Vi = {u,.  ow € LA(Vol), %" € LA(Vol), Onlemo = 0}
To estimate the error we use (4.91) and can directly say for this simple problem L
∫
L
(u’—u;’.)zdxs ∫ o
W' — ur)?dx
(a)
o
where u is the exact solution, u, is the finite element solution, and u, is the interpolant, meaning that «, is considered to be equal to u at the nodal points. Hence, our aim is now to obtain an upper
bound on [y (u' — uf)? dx. Consider an arbitrary element with end points x; and x;., in the mesh. Then we can say that for the exact solution #(x) and x; < x = xi11,
w(x) =u'ly, + (x = x)u" e==
Formulation of the Finite Element Method
250
Chap. 4
where x = x. denotes a chosen point in the element and X is also a point in the element. Let us
choose an x. where u’., = u;, which can always be done because
wi(x) = u(x), wr(xie1) = w(xi41) Then we have for the element
lu'(x) — ut] = h( max lu"l) O=x= O=x= I
Po
x
P,=1i(3x21) Py = 1(5x% — 3x)
P, = 1(35x* — 30x + 3) (n+
)P,y
= @n + xP,
— nP,,
Calculate the stiffness matrix (K), and corresponding load vector of the element for p = 1. Let us first note that these interpolation functions fulfill the requirements of monotonic convergence: the displacement continuity between elements is enforced, and the functions are complete (they can represent the rigid body mode and the constant strain state). This follows because the functions in (a) fulfill these requirements and the functions in (b) merely add
higherorder displacement variations within the element with 2, = Qatx = = 1,i = 3. The stiffness matrix and load vector of the element are obtained using (4.19) and (4.20).
Hence, typical elements of the stiffness matrix and load vector are
o
dhdhy
K, i = ∫∙↕ AE——d. dx dx’x
R?=f
(dd
f(xX)h: dx
The evaluation of (d) gives 1
Ky ==
AE
1
1
Zero 1
2.,
entries
©
"2
"2(p+1)> 5 MPa.
(c)
“eyeball norm” for the accuracy of the stress prediction 7}; achieved with a given finite element mesh. In linear analysis, the finite element stress values can be calculated using the relation 7" = CBii at any point in the element; however, this evaluation is relatively expensive and hardly possible in general nonlinear analysis (including material nonlinear effects). An adequate approach is to use the integration point values to bilinearly interpolate over the corresponding domain of the element. Figure 4.16 illustrates an example in twodimensional analysis. An alternative procedure for obtaining an approximation to the error in the calculated stresses 7/ is to first find some improved values (7})impr. and then evaluate and display ∆⊤∣∎∫ =
Tl}'lj 
(Tg')impn
(4108)
The display can again be achieved effectively using the isoband procedure discussed above. Improved values might be found by simply averaging the stress values obtained at the nodes using the procedure indicated in Fig. 4.16 or by using a least squares fit over the integration point values of the elements (see E. Hinton and J. S. Campbell [A]). The least squares procedure might be applied over patches of adjacent elements or even globally over
Sec. 4.3
257
Convergence of Analysis Results Domain over which stresses are interpolated bilinearly using the four Gauss point values (3 x 3 Gauss integration is used)
b=‘/§a {see Section 5.5.3)
Gauss point
2a
Figure 4.16
Interpolation of stresses from Gauss point stresses
a whole mesh. However, if the domain over which the least squares fit is applied involves many stress points, the solution will be expensive and, in addition, a large error in one part of the domain may affect rather strongly the least squares prediction in the other parts. Another consideration is that when using the direct stress evaluation in (4.12), the stresses are frequently more accurate at the numerical integration points used to evaluate the element matrices (see Section 5.5) than at the nodal points. Hence, for a least squares fit, it can be of value to use functions of order higher than that of the stress variations obtained from the assumed displacement functions because in this way improved values can be expected. We demonstrate the nodal point and least squares stress averaging in the following example. EXAMPLE 4.27: Consider the mesh of ninenode elements shown in Fig. E4.27. Propose reasonable schemes for improving the stress results by nodal point averaging and least squares fitting. Let 7be a typical stress component. One simple and frequently effective way of improving the stress results is to bilinearly extrapolate the calculated stress components from the integration points of each element to node i. In this way, for the situation and node { in Fig. E4.27, four
values for each stress component are obtained. The mean value, say (7%) iean, of these four values is then taken as the value at nodal point i. After performing similar calculations for each nodal point, the improved value of the stress component over a typical element is (Th)impr.
=
∑
∣↿∎⋅≕≺⊺⇂∣⋟⋮↿↿⊧⊟∥
(a)
where the 4; are the displacement interpolation functions because the averaged nodal values are deemed to be more accurate than the values obtained simply from the derivatives of the displacements (which would imply that an interpolation of one order lower is more appropriate).
The key step in this scheme is the calculation of (7*)iean. Such an improved value can also be extracted by using a procedure based on least squares. Consider the eight nodes closest to node i, plus node i, and the values of the stress component of interest at the 16 integration points closest to node i (shown in Fig. E4.27). Let (T")lutege. be the known values of the stress component at the integration points,j = 1, ..., 16, and let (7*)%.4 be the unknown
values at the nine nodes (of the domain corresponding to the
258
Formulation of the Finite Element Method
Chap. 4
Integration point
Figure E4.27
Mesh of ninenode elements. Integration points near node i are also shown.
integration points). We can use the least squares procedure (see Section 3.3.3) to calculate the
values (7")%.q.. by minimizing the errors between the given integration point values and the values
calculated at the same points by interpolation from the nodal point values (T")hoges.
—J——[E (Yinege. — (frh)l.'m,g,y] =0 AT, oes =1
(b
k=1,...,9 9
where
(i"h)lmtegr.
=
∑ k=1
hk
(Th ﬁodes
(C)
at integr. point j
Note that in (c) we evaluate the interpolation functions at the 16 integration stations shown in Fig. E4.27. The relations in (b) and (c) give nine equations for the values (T")foses, k = 1, ..., 9. We solve for these values but accept only the value at node i as the improved stress value, which is now our value for (%) {,an in (a). The same basic procedure is used for all nodes to arrive at nodal “mean” values, so that (a) can be used for all elements.
A least squares procedure clearly involves more computations, and in many cases the simpler scheme of merely extrapolating the Gauss values and averaging at the nodes as described above is adequate. Of course, we presented in Fig. E4.27 a situation of four equal square elements. In practice, the elements are generaily distorted and fewer or more elements may couple into the node i. Also, element noncorner nodes and special mesh topologies at boundaries need to be considered. We emphasize that the calculation of an error measure and its display is a most important aspect of a finite element solution. The quality of the finite element stress solution 74 should be known. Once the error is acceptably small, values of stresses that have been smoothed, for example, by nodal point or least squares averaging, can be used with confidence.
Sec. 4.3
Convergence of Analysis Results
259
These error measures are based on the discontinuities of stresses between elements.
However, for highorder elements (of order 4 and higher), such discontinuities can be small and yet the solution is not accurate because the differential equations of equilibrium of stresses within the elements are not satisfied to sufficient accuracy. In this case the error measure should also include the element internal stress imbalance (4.106). Once an error measure for the stresses has been calculated in a finite element solution
and the errors are deemed to be too large, a procedure needs to be used to establish a new mesh
(with a refined discretization in certain areas, derefinement in other areas, and
possibly new element interpolation orders). This process of new mesh selection can be automatized to a large degree and is important for the widespread use of finite element analysis in CAD 4.3.7
(see Section 1.3).
Exercises
4.25, Calculate the eight smallest eigenvalues of the fournode shell element stiffness matrix available in a finite element program and interpret each eigenvalue and corresponding eigenvector. (Hint: The eigenvalues of the element stiffness matrix can be obtained by carrying out a frequency solution with a mass matrix corresponding to unit masses for each degree of freedom.) 4.26. Show that the strain energy corresponding to the displacement error e,, where e, = u — u,, is equal to the difference in the strain energies, corresponding to the exact displacement solution u and the finite element solution u,. 4.27. Consider the analysis problem in Example 4.6. Use a finite element program to perform the convergence study shown in Fig. 4.12 with the ninenode and fournode (Lagrangian) elements. That is, measure the rate of convergence in the energy norm and compare this rate with the theoretical results given in Section 4.3.5. Use N = 2, 4, 8, 16, 32; consider N = 32 to be the
limit solution, and use uniform and graded meshes. 4.28. Perform an analysis of the cantilever problem shown using a finite element program. Use a twodimensional plane stress element idealization to solve for the static response. (a) Use meshes of fournode elements. (b) Use meshes of ninenode elements.
In each case construct a sequence of meshes and identify the rate of convergence of strain energy. P
1
6 — E = 200,000 v=03
1
260
Formulation of the Finite Element Method
Chap. 4
Also, compare your finite element solutions with the solutions using BernoulliEuler and
Timoshenko beam theories (see S. H. Crandall, N. C. Dahl, and T. J. Lardner [A] and Section 5.4.1).
4.29. Consider the threenode bar element shown. Construct and plot the displacement functions of the element for the following two cases: for case 1:
h.=1latnodei,i=1,2,3
=(Qatnodej # i for case 2:
hi=1latnodei,i=1,2
=Q0atnodej#i,j=1,2 hs =
1 at node3
hs = Qatnode 1, 2
We note that the functions for case 1 and case 2 contain the same displacement variations, and hence correspond to the same displacement space. Also, the sets of functions are hierarchical because the threenode element contains the functions of the twonode element. 1
1
3
2
’_,;
+1
4.30. Consider the eightnode element shown. Idenﬁfy the terms of the Pascal triangle present in the element interpolations.
B =31+ 00 +y,h=4{1200+y) hy =51 — 01— y) he = 3(1 + (1  y) hs = (1 + y)a(x), hs = $(1 — x)a(y) By =51
— y)(x), hs = 3(1 + x)a(y)
where ¢, is defined in Example 4.26.
Sec. 4.4
Incompatible and Mixed Finite Element Models
4.31. A pelement
of order p = 4 is obtained by using
261
the following
displacement
functions.
h,, i =1, 2,3, 4, as for the basic fournode element (with corner nodes only; see Example 4.6). hi,i = 5,..., 16 to represent side modes.
side 1:
AV =11
+ y)(x);
i=5913;=23,4
side 2:
A® =301
— x)¢i(y);
i=6,10,14;=2,3,4
side 3:
A
=51
— y)di(x);
i=17,11,15j=2,3,4
side 4:
A® =11
+ X)(y);
i=8,12,16;j=2,3,4
where ¢., ¢, and ¢4 have been defined in Example 4.26. A;; to represent an internal mode ki = (1 — x)(1
= y?)
Identify the terms of the Pascal triangle present in the element interpolations.
Side 2
3
Side 3
4
4.32. Consider the analysis problem in Example 4.6. Use a finite element program to solve the problem with the meshes of ninenode elements in Exercise 4.27 and plot isobands of the von Mises stress and the pressure (without using stress smoothing). Hence, the isobands will dispiay stress discontinuities between elements. Show how the bands converge to continuous stress bands over the cantilever plate.
4.4 INCOMPATIBLE AND MIXED FINITE ELEMENT MODELS In the previous sections we considered the displacementbased finite element method, and the conditions imposed so far on the assumed displacement (or field) functions were com
pleteness verges in pleteness condition
and compatibility. If these conditions are satisfied, the calculated solution conthe strain energy monotonically (i.e., onesided) to the exact solution. The comcondition can, in general, be satisfied with relative ease. The compatibility can also be satisfied without major difficulties in C° problems, for example, in
262
Formulation of the Finite Element Method
Chap. 4
plane stress and plane strain problems or in the analysis of threedimensional solids such as dams. Yet, in the analysis of shell problems, and in complex analyses in which completely different finite elements must be used to idealize different regions of the structure, compatibility may be quite impossible to maintain. However, although the compatibility requirements are violated, experience shows that good results are frequently obtained. Also, in the search for finite elements it was realized that for shell analysis and the analysis of incompressible media, the pure displacementbased method is not efficient. The difficulties in developing compatible displacementbased finite elements for these problems that are computationally effective, and the realization that by using variational approaches many more finite element discretizations can be developed, led to large research efforts. In these activities various classes of new types of elements have been proposed, and the amount of information available on these elements is voluminous. We shall not present the various formulations in detail but only briefly outline some of the major ideas that have been used and then concentrate upon a formulation for a large class of problems—the analysis of almost incompressible media. The analysis of plate and shell structures using many of the concepts outlined below is then further addressed in Chapter 5. 4.4.1
Incompatible DisplacementBased Models
In practice, a frequently made observation is that satisfactory finite element analysis results have been obtained although some continuity requirements between displacementbased elements in the mesh employed were violated. In some instances the nodal point layout was such that interelement continuity was not preserved, and in other cases elements were used that contained interelement incompatibilities (seg Example 4.28). The final result was the same in either case, namely, that the displacemenhts or their derivatives between elements
were not continuous to the degree necessary to satisfy all compatibility conditions discussed in Section 4.3.2. Since in finite element analysis using incompatible (nonconforming) elements the requirements presented in Section 4.3.2 are .not satisfied, the calculated total potential energy is not necessarily an upper bound to the exact total potential energy of the system, and consequently, monotonic convergence is not ensured. However, having relaxed the objective of monotonic convergence in the analysis, we still need to establish conditions that will ensure at least a nonmonotonic convergence. Referring to Section 4.3, the element completeness condition must always be satisfied, and it may be noted that this condition is not affected by the size of the finite element. We recall that an element is complete if it can represent the physical rigid body modes (but the element matrix has no spurious zero eigenvalues) and the constant strain states. However, the compatibility condition can be relaxed somewhat at the expense of not obtaining a monotonically convergent solution, provided that when relaxing this requirement, the essential ingredients of the completeness condition are not lost. We recall that as the finite element mesh is refined (i.e., the size of the elements gets smaller), each element
should approach a constant strain condition. Therefore, the second condition on convergence of an assemblage of incompatible finite elements, where the elements may again be of any size, is that the elements together can represent constant strain conditions. We should
Sec. 4.4
Incompatible and Mixed Finite Element Models
263
note that this is not a condition on a single individual element but on an assemblage of
elements. That is, although an individual element is able to represent all constant strain states, when the element is used in an assemblage, the incompatibilities between elements may prohibit constant strain states from being represented. We may call this condition the completeness condition on an element assemblage. As a test to investigate whether an assemblage of nonconforming elements is com
plete, the patch test has been proposed (see B. M. Irons and A. Razzaque [Al]). In this test a specific element is considered and a patch of elements is subjected to the minimum displacement boundary conditions to eliminate all rigid body modes and to the boundary nodal point forces that by an analysis should result in constant stress conditions. If for any patch of elements the element stresses actually represent the constant stress conditions and all nodal point displacements are correctly predicted, we say that the element passes the patch test. Since a patch may also consist of only a single element, this test ensures that the element itself is complete and that the completeness condition is also satisfied by any element assemblage, The number of constant stress states in a patch test depends of course on the actual number of constant stress states that pertain to the mathematical model; for example, in plane stress analysis three constant stress states must be considered in the patch test, whereas in a fully threedimensional analysis six constant stress states should be possible. Fig. 4.17 shows a typical patch of elements used in numerical investigations for various problems. Here of course only one mesh with distorted elements is considered, whereas in fact any patch of distorted elements should be analyzed. This, however, requires an analytical solution. If in practice the element is complete and the specific analyses shown here produce the correct results, then it is quite likely that the element passes the pitcit test. When considering displacementbased elements with incompatibilities, if #f¢ patch test is passed, convergence is ensured (although convergence may not be monatonic and convergence may be slow). The patch test is used to assess incompatible finite element meshes, and we may note that when properly formulated displacementbased elements are used in compatible meshes, the patch test is automatically passed. Figure 4.18(a) shows a patch of eightnode elements (which are discussed in detail in Section 5.2). The tractions corresponding to the plane stress patch test are also shown. The elements form a compatible mesh, and hence the patch test is passed. However, if we next assign to nodes 1 to 8 individual degrees of freedom for the adjacent elements [e.g., at node 2 we assign two u and v degrees of freedom each for elements 1 and 2] such that the displacements are not tied together at these nodes (and therefore displacement incompatibilities exist along the edges), the patch test is not passed. Figure 4.18(b) gives some results of the solution. The example in Fig. 4.18(b) uses, in essence, an element that was proposed by E. L. Wilson, R. L. Taylor, W. P. Doherty, and J. Ghaboussi [A]. Since the degrees of freedom of the midside nodes of an element are not connected to the adjacent elements, they can be statlcally condensed out at the element level (see Section 8.2.4) and a fournode element is
obtained. However, as indicated in Fig. 4.18(b), this element does not pass the patch test. In the following example, we consider the element in more detail, first as a square element
{10, 10)
Z
(10, 0)
(a) Patch of elements, twodimensional elements, plate bending elements, or side view of threedimensionel elements. Each quedrilaterel domain
A
aT
T
represents an element; for triangular and tetrahedral elements, each quadrilateral domain is further subdivided
Plane stress and plene strain: T, Ty, Txy
Axisymmetric; here perform the test with R —
constent; In threedimensional analysis the
additionel three stress conditions 75, Tz, 7yz constant are tested
↡ My
YI
@'
T
{This test also produces bending) Plate bending (see Section 5.4.2) (b) Stress conditions to be tested Figure 4.17
264
Patch tests
Constant normal traction t, Constant shear traction t,
Constant normal traction t,
y.v
(a) Patch test of compatible mesh of 8noda elements (discussed in Section 5.3.1). The patch test is paased; that ia, all calculated element stresses are equal to tha applied tractions
va{" = v displacement at node 2 associated with element 1, etc.
T =2.74 Tyy = 0.19 Ty = 0.08 Ty = 0.44
Tyy
E=2.1x108
=~0.11
Tyy=0.10
v=03 Thickness = 1.0
Ty = 0.08 Tyy =0.58 Ty = 0.02
10.0 Txx =—0.52 Tyy = =0.40 Tyy = 0.85
y
Analytical solution: %
Txx =
3.44
Tyy=004
Ty=018
Txx =
2.0
1,.€
corresponding to term 2:
— (fv E’E dV)‘r]
([map
o  ([
corresponding to term 3:
61"1’
’
epm
6 Then we can write
u = Hg;
€&x = B,i
©
Sec. 4.4
Incompatible and Mixed Finite Element Models
275
Substituting into (b) and (c), we obtain Kuu
Kue
ﬁ
R
[K,L K][e]  [ 03] where
K. = ∫ B[EBde;
@
K.. = ∫ ≖⋝≖⊺∊∌≨∖⊱↙↿∨
⋁
⋁
Rs= ∫ HTdeV
∥
↧⊂∊∊− j (B )'GBdV,
V
V
We can now use static condensation on & to obtain the final element stiffness matrix: K=
Kuu

KueK;El
Kze
In our case, we have
0
zfL ——
zfL ~—
B = (1]
 @z 2@ 2 Gh L
th so that
2 K= —Gh L
—
Gh 2
—Gh L
2L
~Gh 2
Gh 2
2 Gh L
~Gh 2
c @8 @y 2
12L
2
12L
©
12L
It is interesting to note that a pure displacement formulation would give a very similar stiffness matrix. The only difference is that the circled terms would be GhL/3 on the diagonal and GhL/6 in the offdiagonal locations. However, the element predictive capability of the pure displacementbased formulation is drastically different, displaying a behavior that is much too stiff when the element is thin (we discuss this phenomenon in Sections 4.5.7 and 5.4.1). Note that if we assume a displacement vector corresponding to section rotations only,
i=00
a
0
a
then using () the element displays bending stiffness only, whereas the pure displacementbased element shows an erroneous shear contribution. Let us finally note that the stiffness matrix in (€) corresponds to the matrix obtained in the mixed interpolation approach discussed in detail in Section 5.4.1. Namely, if we use the last equation in (d), which corresponds to the equation (c), we obtain AS
Yxz
=
=
Wz—W_91+02
L
2
276
Formulation of the Finite Element Method
Chap. 4
which shows that the assumed shear strain value is equal to the shear strain value at the midpoint of the beam calculated from the nodal point displacements.
As pointed out above, the HuWashizu principle provides the basis for the derivation of various variational principles, and many different mixed finite element discretizations can be designed. However, whether a specific finite element discretization is effective for practical analysis depends on a number of factors, particularly on whether the method is general for a certain class of applications, whether the method is stable with a sufficiently high rate of convergence, how efficient the method is computationally, and how the method compares to alternative schemes. While mixed finite element discretizations can offer some advantages in certain analyses, compared to the standard displacementbased discretization, there are two large areas in which the use of mixed elements is much more efficient
than the use of pure displacementbased elements. These two areas are the analysis of almost incompressible media and the analysis of plate and shell structures (see the following sections and Section 5.4). 4.4.3 Mixed Interpolation—Displacement/Pressure Formulations for Incompressible Analysis
The displacementbased finite element procedure described in Section 4.2 is very widely used because of its simplicity and general effectiveness. However, there are two problem areas in which the pure displacementbased finite elements are not sufficiently effective, namely, the analysis of incompressible (or almost incompressible) media and the analysis of plates and shells. In each of these cases, a mixed interpolation approach—which can be thought of as a special use of the HuWashizu variational principle (see Example 4.30)—is far more efficient. We discuss the mixed interpolation for beam, plate, and shell analyses in Section 5.4, and we address here the analysis of incompressible media. Although we are dealing with the solution of incompressible solid media, the same basic observations are also directly applicable to the analysis of incompressible fluids (see Section 7.4). For example, the elements summarized in Tables 4.6 and 4.7 (later in this section) are also used effectively in fluid flow solutions. The Basic Differential Equations for Incompressible Analysis
In the analysis of solids, it is frequently necessary to consider that the material is almost incompressible. For example, some rubberlike materials, and materials in inelastic conditions, may exhibit an almost incompressible response. Indeed, the compressibility effects may be so small that they could be neglected, in which case the material would be idealized as totally incompressible. A basic observation in the analysis of almost incompressible media is that the pressure is difficult to predict accurately. Depending on how close the material is to being incompressible, the displacementbased finite element method may still provide accurate solutions, but the number of elements required to obtain a given solution accuracy is usually far greater than the number of elements required in a comparable analysis involving a compressible material.
Sec. 4.4
Incompatible and Mixed Finite Element Models
277
To identify the basic difficulty in more detail, let us again consider the threedimensional body in Fig. 4.1. The material of the body is isotropic and is described by Young’s modulus E and Poisson’s ratio v. Using indicial notation, the governing differential equations for this body are (see Example 4.2) wy+
f2=0
throughout the volume V of the body
(4.123)
on S¢
f¥
= w; = u
(4.122)
on S,
(4.124)
If the body is made of an almost incompressible material, we anticipate that the volumetric strains will be small in comparison to the deviatoric strains, and therefore we use the constitutive relations in the form (see Exercise 4.39)
where « is the bulk modulus, K=
£
4.126
3 = 29)
(4.126)
€év is the volumetric strain, €y =
Av
€
.
.
.
= —V( = €, + €, + €, in Cartesian coordinates)
(4.127)
8, is the Kronecker delta,
’
128
€;; are the deviatoric strain components,
(4.129)
€ = € — ⊲⋮−∋∀ 8y and G is the shear modulus, =
E
m
(4.130)
We also have for the pressure in the body, P = —K€ where
p=
XX
+
+
T—;" ( = — I_?__T_
2z
(4.131) .
.
.
in Cartesian coordmates)
(4.132)
Formulation of the Finite Element Method
278
Chap. 4
Now let us gradually increase k (by increasing the Poisson ratio v to approach 0.5).
Then, as k increases, the volumetric strain ey decreases and becomes very small. In fact, in total incompressibility » is exactly equal to 0.5, the bulk modulus is infinite, the volumetric strain is zero, and the pressure is of course finite (and of the order of the applied boundary tractions). The stress components are then expressed as [see (4.125) and
(4.131)]
and the solution of the governing differential equations (4.122) to (4.124) now involves using the displacements and the pressure as unknown variables.
In addition, special attention need also be given to the boundary conditions in (4.123) and (4.124) when material incompressibility is being considered and the displacements are prescribed on the complete surface of the body, i.e., when we have the special case S, = §, S¢ = 0. If the material is totally incompressible, a first condition is that the prescribed displacements 7 must be compatible with the zero volumetric strain throughout the body. This physical observation is expressed as € =0
hence,
throughout V
(4.134)
∫ e,,dV= ∫ uwendS=0
(4.135)
v
S
where we used the divergence theorem and n is the unit normal vector on the surface of the body. Hence, the displacements prescribed normal to the body surface must be such that the
volume of the body is preserved. This condition will of course be automatically satisfied if the prescribed surface displacements are zero (the particles on the surface of the body are not displaced). Assuming that the volumetric strain/boundary displacement compatibility is satisfied, for the case S. = S, the second condition is that the pressure must be prescribed at some point in the body. Otherwise, the pressure is not unique because an arbitrary constant pressure does not cause any deformations. Only when both these conditions are fulfilled is the problem well posed for solution. Of course, the condition of prescribed displacements on the complete surface of the body is a somewhat special case in the analysis of solids, but we encounter an analogous situation frequently in fluid mechanics. Here the velocities may be prescribed on the complete boundary of the fluid domain (see Chapter 7). Although we considered here a totally incompressible medium, it is clear that these considerations are also important when the material is only almost incompressible—a
violation of the conditions discussed will lead to an illposed problem statement. Of course, these observations also pertain to the use of the principle of virtual work. Let us consider the simple example shown in Fig. 4.19. Since only volumetric strain energy is present, the principle of virtual work gives for this case, ∫ EVKeVdV= v
−∫
8¢
Bsp*dS
(4.136)
101010101010101010101010]
() I L3CCCICICICICICICI{
Bulk modulus x Pressure p
279
I/
≣
Y, v
Incompatible and Mixed Finite Element Models
(D O0O0O0O0O0O000000)
Sec. 4.4
Figure 4.19
Block of material in plane
strain condition, subjected to uniform
surface pressure p*
If the bulk modulus « is finite, we obtain directly from (4.136), *
and
oS = —’lK—
4.137)
p = p*
(4.138)
However, if ks infinite, we need to use instead of (4.136) the following form of the principle of virtual work, with the pressure p unknown, ∫ Ev(—p)dV= V
−∫
Esp*dS
(4.139)
Sf
and we again obtainp = p*. Of course, the solution of (4.139) does not use the constitutive relation but only the equilibrium condition. The Finite Element Solution of Almost Incompressible Conditions
The preceding discussion indicates that when pursuing a pure displacementbased finite element analysis of an almost incompressible medium, significant difficulties must be expected. The very small volumetric strain, approaching zero in the limit of total incompressibility, is determined from derivatives of displacements, which are not as accurately predicted as the displacements themselves. Any error in the predicted volumetric strain will appear as a large error in the stresses, and this error will in turn also affect the displacement prediction since the external loads are balanced (using the principle of virtual work) by the stresses. In practice, therefore, a very fine finite element discretization may be required to obtain good solution accuracy.
280
Formulation of the Finite Element Method
p = 0.006 MPa
↜∢∎∎∎∎∎∎∎∎∎∎∎∎∎∎∎∎∎ 10 mm
v=0.3, 0.499
S ymmetry about
centerline of bracket
bt
ot
i
E =55 MPa
(a)
Geometry, material data, applied loading, and the coarse sixteen element mesh
Figure 4.20 Analysis of cantilever bracket in plane strain conditions. Ninenode displacementbased elements are used. The 16 X 64 = 1024eleinent inesh is obtained by dividing each element of the 16element mesh into 64 elements. Maximuin principal stress results are shown using the band representation of Fig. 4.15. Also, (01)max is the predicted naxiinuin
value of the maxiinum principal stress, and & is defined in (a).
Chap. 4
Sec. 4.4
Incompatible and Mixed
Finite Element Models
281
1EMA
1
Tomn
g
{o1)max = 0.5050 8 = 15682
t
g
1
f
+
T
e
{o1)max = 0.6056 = 1.669
M
[T
L (b}
Displacementbased element solution resuits for the case Poisson’s ratio v = 0.30. Sixteen element and 16 X 64 element mesh resulis Figure 4.20
(continued)
Figure 4.20 shows some results obtained in the analysis of a cantilever bracket subjected to pressure loading. We consider plane strain conditions and the cases of Poisson’s ratio v = 0.30 and v = 0.499. In all solutions, ninenode displacementbased elements have been used (with 3 X 3 Gauss integration; see Section 5.5.5). A coarse mesh and a very
fine mesh are used, and Fig. 4.20(a) shows the coarse idealization using only 16 elements. The solution results for the maximum principal stress oy are shown using the isoband representation discussed in Section 4.3.6. Here we have selected the bandwidth so as to be able to see the rather poor performance of the displacementbased element when the Poisson ratio is close to 0.5. Figure 4.20(b) shows that when
» = 0.30, the element stresses are
reasonably smooth across boundaries for the coarse mesh and very smooth for the fine mesh. Indeed, the coarse idealization gives a quite reasonable stress prediction. However,
Formulation
282
of the Finite Element
Method
Chap. 4
SIGMAT TTMF .G00 03500
0.14%0
(0‘1)max
&=
=
1.955
1.044
== {o1)max
&=
(c)
1.343
1.363
Displacementbased element solution results for the case Poisson’s ratio v = 0.499. Sixteen
element and
16 X 64 element
Figure 4.20
when
=
mesh
results
(continued)
» = 0.499, the same meshes of ninenode displacementbased elements result into
poor stress predictions [see Fig. 4.20(c)]. Large stress fluctuations are seen in the individual
elements of the coarse mesh and the fine mesh.'® Hence, in summary, we see here that the
displacementbased element used in the analysis is effective when » = 0.3, but as v approaches 0.5, the stress prediction becomes very inaccurate. This discussion indicates what is very desirable, namely, a finite element formulation
which gives essentially the same accuracy in results for a given mesh irrespective of what Poisson’s ratio is used, even when v is close to 0.5. Such behavior is observed if for the finite 5We discuss briefly in Section 5.5.6 the use of “reduced integration.” If in this analysis the reduced integration of 2 X 2 Gauss integration is attempted, the solution cannot be obtained because the resulting stiffness matrix is singular.
Sec. 4.4
Incompatibla and Mixed Finite Element Models
283
element formulation the predictive capability of displacements and stresses is independent of the bulk modulus used. We refer to finite element formulations with this desirable behavior as nonlocking, whereas otherwise the finite elements are locking. The term “locking” is based upon experiences in the analysis of beams, plates, and shells (see Section 5.4.1), where an inappropriate formulation—one that locks—results in displacements very much smaller than those intuitively expected for a given mesh (and calculated with an appropriate formulation; see, for example, Fig. 5.20). In the analysis of almost incompressible behavior, using a formulation that locks, the displacements are not necessarily that much in error but the stresses (the pressure) are very inaccurate. We note that the pure displacement formulation generally locks in almost incompressible analysis. These statements are discussed more precisely in Section 4.5. Effective finite element formulations for the analysis of almost incompressible behavior that do not lock are obtained by interpolating displacements and pressure. Figure 4.21 shows the results obtained in the analysis of the cantilever bracket in Fig. 4.20 with a displacement/pressure formulation referred to as u/p formulation using the 9/3 element (see below
for the explanation
of the formulation
and the element).
We
see that the
isobands of the maximum principal stress have in all cases the desirable degree of smoothness and that the stress prediction does not deteriorate when Poisson’s ratio v approaches 0.5. To introduce the displacement/pressure formulations, we recall that in a pure displacement formulation, the evaluation of the pressure from the volumetric strain is difficult when « is large (in comparison to G) and that when a totally incompressible condition is considered, the pressure must be used as a solution variable [see (4.133)]. It therefore appears reasonable to work with the unknown displacements and pressure as solution variables when almost incompressible conditions are analyzed. Such analysis procedures,
if properly formulated, should then also be directly applicable to the limit of incompressible conditions. The basic approach of displacement/pressure finite element formulations is therefore to interpolate the displacements and the pressure. This requires that we express the principle of virtual work in terms of the independent variables u and p, which gives
∫ E'Tde— ∫ Gpdv=a v
(4.140)
v
where, as usual, the overbar indicates virtual quantities, R corresponds to the usual external virtual work [ is equal to the righthand side of (4.7)], and S and €' are the deviatoric stress and strain vectors, S=1+ € =€~
pd 1
∃∊⋁∂
(4.141) (4.142)
where 8 is a vector of the Kronecker delta symbol [see (4.128)]. Note that using the definition of p in (4.131), a uniform compressive stress gives a positive pressure and that in the simple example in Fig. 4.19, only the volumetric part of the
internal virtual work contributed. In (4.140) we have separated and then summed the deviatoric strain energy and the bulk strain energy. Since the displacements and pressure are considered independent vari
Formulation of the Finite Element Method
284
Chap. 4
SIGMAP TIMF 1.000 C.anee
{o)max = 0.4940 5 = 1.625
e
tt
{o1)max
=
0.6000
= 1.680
(a)
Bands of maximum principal stress. Case of Poisson's ratio v = 0.30. Sixteen and 16 X 64 element mesh results
in wn sho is ket Brac ns. itio cond in stra e plan in ket brac er ilev cant of ysis Anal 4.21 ure Fig Fig. 4.20(a). Same meshes as in Fig. 4.20 are used but with the ninenode mixed interpolated element (the 9/3 element). Compare the results shown with those given in Fig. 4.20.
ables, we need another equation to connect these two solution variables. This equation is provided by (4.131) written in integral form (see Example 4.31), (4.143)
(see s ple nci pri al ion iat var m fro ly mal for e mor d ive der be also can These basic equations
g owin foll the in ons ati equ c basi the ive der We . [A]) Key W. S. and [A] nn rma Her L. R. example from the HuWashizu functional.
Sec. 4.4
Incompatible and Mixed Finite Element Models
285 BlGMa
1] )L aen
I
).
(01)max = 0.4983 = 1.349
G1GVA T S1ME 1L e
1]
UL 3ho0
(01)max = 0.5998 8= 1.393
(b)
Bands of maximum principal stress. Case of Poisson's ratio v = 0.499. Sixteen and 16 X 64 element mesh resulis
Figure 4.21
(continued)
EXAMPLE 4.31: Derive the u/p formulation from the HuWashizu variational principle. The derivation is quite analogous to the presentation in Example 4.30 where we considered a mixed interpolation for a beam element. We start by letting T = Ce in (4.114) to obtain the HellingerReissner functional,
1
,’!‘,R(u,e)= −∫
3 € Ce dv + ∫ ∊⊺∁∂⋮⇂∥↿∨⇀ v
v
.
J uw'ft v 4
— ∫
∐⊾⋝⋅↗⊺↥⋅⋅⋝∎↗⋅∠↿⊧⋝
(a)
Sj
where we assume that the displacement boundary conditions are satisfied exactly (hence, also the displacement variations will be zero on the surface of prescribed displacements).
286
Formulation of the Finite Element Method
Chap. 4
Next we establish the deviatoric and volumetric contributions and postulate that the deviatoric contribution will be evaluated from the displacements. Hence, we can specialize (a) into 2
#(u, p) = ∫ ↨∈∣⊤∁∣∈∣ dv — ∫ ↨∣↗− dv — ∫ pevdV — ∫ u’f”dV— v2
v2
K
v
v
∫
s
w7t
ds (b)
where the prime denotes deviatoric quantities, €y is the volumetric strain evaluated from the displacements, p is the pressure, and k is the bulk modulus. Note that whereas in (a) the independent variables are u and €, in (b) the independent variables are u and p. Invoking the stationarity of II%; with respect to the displacements and the pressure, we obtain
∫ SG’TC'e'dV—
∫ péey dV = R
v
v
p
and
_ + ev)Sp dav =20
f(— v\K
where R corresponds to the virtual work of the externally applied loading [see (4.7)]. It is interesting to note that we may also think of (b) as the total potential in terms of the displacements and the pressure plus a Lagrange multiplier term that enforces the constraint between the volumetric strains and the pressure, =
I
1
= −∫
—e’TC’e’dV+ v
−∫ s
2
1
p?
∫ —p—dV— v
2 K
∫ u’f? dv v
(C)
uszfsde— ∫ A(ev + I—’) av v
K
In (c) the last integral represents the Lagrange multiplier constraint, and we find A = p.
To arrive at the governing finite element equations, we can now use (4.140) and (4.143) as in Section 4.2.1, but in addition to interpolating the displacements we also interpolate the pressure p. The discussion in Section 4.2.1 showed that we need to consider the formulation of only a single element; the matrices of an assemblage of elements are then formed in a standard manner. Using, as in Section 4.2.1,
u = Hi
(4.144)
we can calculate
€ = Boii;
& = Byii
(4.145)
The additional interpolation assumption is
p=Hp
(4.146)
where the vector P lists the pressure variables [(see the discussion following (4.148)].
= I3I 0
Substituting from (4.144) to (4.146) into (4.140) and (4.143), we obtain
where
K. = ∫ B{,C'BDdV V
Sec. 4.4
287
Incompatible and Mixed Finite Element Models
K.,
K} =~ ∫
€H,,dv
(4.148)
KPP
Il
v

1 T
∫ VH,,KH,,dV
andC’isthestressstrainmatrixforthedeviatoricstressandstraincomponents. The relations in (4.144) to (4.148) give the basic equations for formulating elements with displacements and pressure as variables. The key question for the formulation is now, What pressure and displacement interpolations should be used to arrive at effective elements? For example, if the pressure interpolation is of too high a degree compared to the displacement interpolation, the element may again behave as a displacementbased element and not be effective. Considering for the moment only the pressure interpolation, the following two main possibilities exist and we label them differently.
The u/p formulation. In this formulation, the pressure variables pertain only to the specific element being considered. In the analysis of almost incompressible media (as so far discussed), the element pressure variables can be statically condensed out prior to the element assemblage. Continuity of pressure is not enforced between elements but will be a result of the finite element solution if the mesh used is fine enough.
The u/pc formulation. The letter “’c” denotes continuity in pressure. In this formulation, the element pressure is defined by nodal pressure variables that pertain to adjacent elements in the assemblage. The pressure variables therefore cannot be statically condensed out prior to the element assemblage. Continuity of pressure between elements is directly enforced and will therefore always be a result of the solution irrespective of whether the mesh used is fine or coarse. Consider the following two elements, one corresponding to each of the formulations. EXAMPLE 4.32: For the fournode plane strain element shown, assume that the displacements are interpolated using the four nodes and assume a constant pressure. Evaluate the matrix expressions used for the u/p formulation. y.v
Young's modulus £ Poisson's ratio v
Figure E4.32
A 4/1 element
Chap. 4
Formulation of the Finite Element Method
288
This element is referfed to as the u/p 4/1 element. In plane strain analysis we have €
—
l(e
€
)
""3"""
, €
+
zilf
29v
3(6’3‘
—
€&y
1 0u 
3—6—)’
=
EYY)
+
la_v
3ax3ay
1 =
_
∊⋁∶∊≍≍⋅⊢
∃
≣−∂−⋅⊊
€y
(a)
ody dvox
Yo 1
1/0u
—=(€&x 3 (e
+
ov = 6y>
+ 3 ( 3x
€y)
and S = C'e’, where 0
0
O
0
0
26 0
G
O O
0
0
0
2G
26 0 =
C,
E s
P
G
—
2(1 + »)
The displacement interpolation is as in Example 4.6,
u = Ha
with
u(x, y) = [u(x’ y)]; o{x, y)
H
—
[hl
0
h2 0
0 =[m
w
ws
ws
0 h
0 h
0
O
hy
hy
¢:+ ¢
h4 4 0
hs 3 0
¥
v
V2
V3
4] (b)
]
hy=3;0 1 +y) he = (1 + (1  y)
b = z(1 + 0)(1 + y); b =;(1 = 01  y);
Using (a) and (b), the straindisplacement interpolation matrices are %hl_x
%hz,x
_%hl,x BV
=
[hl'x
_%hz,x h2,x
...
_
F 0
..
hyy
hyy
−⋚⇂∥⊽≻∙
oo
—3hn
3 k B, =
and
e
el ⋮
−≣↥∣↕≵⋅⊃↗
by,
B,
hys
hy,x
_%hl'y ∣↕↿∙≻≀
∣↕≵∙≻≀
‘%hZ,y ⋅⋅
]
For a constant pressure assumption we have
H, = [1];
B = [po]
Since the degree of freedom P = [po] pertains only to this element and not to the adjacent elements, we can use static condensation to obtain from (4.147) the element stiffness matrix corresponding to the nodal point displacement degrees of freedom only;
K = Ku — K, K, 'K, The element is further discussed in Example 4.38.
Sec. 4.4
289
Incompatible and Mixed Finite Element Models
EXAMPLE 4.33: Consider the ninenode plane strain element shown in Fig. E4.33. Assume that the displacements are interpolated using the nine nodes and that the pressure is interpolated using only the four corner nodes. Refer to the information given in Example 4.32 and discuss the additional considerations for the evaluation of the matrix expressions of this element.
® Displacament noda @ Displacement and pressure node
Young's modulus £ Poisson's ratio v
Figure E4.33
A 9/4c element
This element was proposed by P. Hood and C. Taylor [A]. In the formulation the nodal pressures pertain to adjacent elements, and according to the above element nomination we refer to it as a u/pc element (it is the 9/4c element). The deviatoric and volumetric strains are as given in (a) in Example 4.32. The displacement interpolation corresponds to the nine nodes of the element,
[u(x,y)]_[h;*
oo
hE
i
0
L.
0]__,
L0
...
0
i
B
...
hllo,
ox, ]
@
where the interpolation functions A} are constructed as explained in Section 4.2.3 (or see Section 5.3 and Fig. 5.4). The deviatoric and volumetric straindisplacement matrices are obtained as in Example 4.32. The pressure interpolation is given by
P=[h1
hy
hs
h4]
pz
where the h; are those given in (b) in Example 4.32. A main computational difference between this element and the fournode element discussed in Example 4.32 is that the pressure degrees of freedom cannot be statically condensed out on the element level because the variables py, . . . , p4 pertain to the element we are considering here and to the adjacent elements, thus describing a continuous pressure field for the discretization.
Formulation of the Finite Element Method
290
Chap. 4
Let us now return to the discussion of what pressure and displacement interpolations should be used in order to have an effective element. For instance, in Example 4.32, we used four nodes to interpolate the displacements and assumed
a constant pressure,
and we may
ask whether
a constant pressure
is the
appropriate choice for the fournode element. Actually, for this element, it is a somewhat natural choice because the volumetric strain calculated from the displacements contains linear variations in x and y and our pressure assumption should be of lower order. When higherorder displacement interpolations are used, the choice of the appropriate pressure interpolation is not obvious and indeed much more difficult: the pressure should not be interpolated at too low a degree because then the pressure prediction could be of higher order and hence be more accurate, but the pressure should also not be interpolated at too high a degree because then the element would behave like the displacementbased elements and lock. Hence, we want to use the highest degree of pressure interpolation that does not introduce locking into the element. For example, considering the u/p formulation and biquadratic displacement interpolation (i.e., nine nodes for the description of the displacements), we may naturally try the following cases: 1. Constant pressure, p = po
(9/1 element)
2. Linear pressure, p = po + pix + p2y (9/3 element) 3. Bilinear pressure, p = po + pix + p2y + psxy (9/4 element) and so on, up to a quadratic pressure interpolation {which corresponds to the 9/9 element). These elements have been analyzed theoretically and by use of numerical experiments. Studies of the elements show that the 9/1 element does not lock, but the rate of convergence of pressure (and hence stresses) as the mesh is refined is only of o(h) because
a constant pressure is assumed in each ninenode element. The poor quality of the pressure prediction can of course also have a negative effect on the prediction of the displacements. Studies also show that the 9/3 element is most attractive because it does not lock and the stress convergence is of o(h?). Hence, the predictive capability is optimal since if a biquadratic displacement expansion is used, no higherorder convergence in stresses can be expected. Also, the 9/3 element is effective for any Poisson’s ratio up to 0.5 (but the static condensation of the pressure degrees of freedom is possible only for values of v < 0.5). Hence, we may be tempted to always use the 9/3 element (instead of the displacementbased ninenode element). However, we find in practice that the 9/3 element is computationally slightly more expensive than the ninenode displacementbased element, and when v is less than 0.48, the additional terms in the pressure expansion of the displacementbased element allow a slightly better prediction of stresses. The next u/p element of interest is the 9/4 element, and studies show that this element locks when v is close to 0.50; hence it cannot be recommended for almost incompressible analysis. In an analogous manner, other u/p elements can be constructed, and Table 4.6 summarizes some choices. Regarding these elements, we may note that the fournode twodimensional and eightnode threedimensional elements are extensively used in practice. However, the ninenode twodimensional and 27node threedimensional elements are frequently more powerful.
Sec. 4.4
Incompatible and Mixed Finite Element Models
As indicated in Table 4.6, the Q> — P, and P;
— P, elements are the first members
of two families of elements that may be used. That is, the quadrilateral elements @, — and the triangular elements P,
—
291
P,;,
P,_, n > 2, are also effective elements.
In Table 4.6 we refer to the infsup condition, which we will discuss in Section 4.5. From a computational point of view, the u/p elements are attractive because the element pressure degrees of freedom can be statically condensed out before the elements are assembled (assuming v < 0.5 but possibly very close to 0.5). Hence, the degrees of freedom for the assemblage of elements are only the same nodal point displacements that are also the degrees of freedom in the pure displacementbased solution. However, the u/pc formulation has the advantage that a continuous pressure field is always calculated. Table 4.7 lists some effective elements.
The Finite Element Solution of Totally Incompressible Conditions
If we want to consider the material to be totally incompressible, we can still use (4.140) and (4.143), but we then let k — . For this reason, we refer to this case as the limit problem. Then (4.143) becomes
∫ ∊⋁⊅−∠↿∨⇌∘
(4.149)
v
and (4.147) becomes, correspondingly,
[Il(i.m Koup] [:] _ [1;]
(4.150)
Hence, in the coefficient matrix, the diagonal elements corresponding to the pressure degrees of freedom are now zero. It follows that a static condensation of the element pressure degrees of freedom in the u/p formulation is no longer possible and that the solution of the equations of the complete assemblage of elements needs special considerations (beyond those required in the pure displacementbased solution) to avoid encountering a zero pivot element (see Section 8.2.5). Suitable elements for solution are listed in Tables 4.6 and 4.7. These elements are effective (except for the O, — P, elements) because they have good predictive capability irrespective of how close the behavior of the medium is to a situation of total incompressibility (but the procedure for solving the governing finite element equations must take into account that the elements in K,, become increasingly smaller as total incompressibility is approached). As already noted earlier, we refer to the infsup condition in Tables 4.6 and 4.7. This condition is the basic mathematical criterion that determines whether a mixed finite element discretization is stable and convergent (and hence will yield a reliable solution). The condition was introduced as the fundamental test for mixed finite element formulations by 1. Babuska [A] and F. Brezzi [A] and since then has been used extensively in the analysis of mixed finite element formulations. In addition to the infsup condition, there is also the ellipticity condition which has not received as much attention because frequently—as in the analysis of almost incompressible media—the ellipticity condition is automatically satisfied.
‘uonoipaid ssans 9JeIndoe I0J SUOTIBZIIDIOSIP aul a1mbal Aew vondumsse amssaid
JURISUOD I} NG “UOTIPUOD dnsjur a1y SoYSTIES JUSWIS]D Y,
(§°¢'y uondes Ul JUSWID[A JO UOISSNISIP 93s) uonpuod dnsjur oy £JSTIES 10U S0P JUW]D Y, ‘suorenyonyy amssaid aqissod pue uondwnsse amssard JUEISUOD S} JO IsneIq aleINdoeUl 9 ABUI SISSAIS nq ‘sjuouradeidsip poos Kjqruosear sjo1paxd jusuis[d Y,
syrewoy
WUAURR (I€
siwod [epoN
ln
A
0d =d
EoEo_oQ.N
1/0Z :a@¢ ut 1/8 @z ut
°d  30 UD =% EF
/8 agut
—
↔↻
/¥ acuw
cn&
↩≣⊟∘⊡
⊢≘≣⋮∾∼∾∼≣↘↘∙⋝⋅∼∼⋮⋅∼≋⇀⊾⋅≊∙⊏∾⇖↜∡∾≋∼↘⋅⋮⋝∾≅≌∾⊾⇖≣∙⊸≣≣⋮≃≣↪∾≅∙⊸⋮∾⋮∾∼∾≣∾⋛∾⊶∼ ⋅⊸⋮∾⋮∾↪↦↘∼⇖⋅⊸⋅∼≊⋅⊸⋮∾⋮∾∼∾⇖∖≖≣⋅≣∾⇤↖∾↯≣⋅⋮⋝⊛⋅∢⊔∎−⊠∢⋅−⋅
292
‘uonIpuod
dnsyut atp Ajsnes sjuswIpd ssopv'z=Eu'""d — /d siuoWord JenSueLn Jo Apurey 2 Jo Iaquusur Is1y YL "PIYSNIES 2q O} ST UOPIPUOI dnsjut ayy 1Ry uonIPUOd AP
01 193{qns pasn s sojqeLeA arnssaid Jo saquinu JsaSre]
ayJ, 9[qqnq 21qno 3y 4q payouud Iy sperwoukiod yo
aouds ap sajousp 4 “(s[qqnq 21qnd B IIM) uorsuedxa Juawade[dsip onespenb & Sursn Jusurdpd reyndueLn
reumndo ayy, {v] weney
"V'd PUE X19Zno1) " 99§
“UOTIBZNAIISIP suy armbar ew uondunsse amssaid jueisuod 9y Ing ‘uonipuod dnsJur oy soysmes
wewapd oy, [v] wemey "V'd PUE XI9ZNox) ‘W 99§
⊆∾↽∪∘⇁∽↼
⋅
↩⊏∾∈∾−∾ .. ⊝⋔→⋊↻−↶↔⋮∾⋈⊡↞∘∙−∾≅∾∘
∓↻∘∽∨≣⊶⊔⊶↔⋮∘∘ ↶↔≣⋅↬⋮
∘≡⋋⊔↔∽↝⊔∾∽∽⋮∘⊟↻−↻↻∽≗∶⋜ ⋅∾⋯≂↜−∣≂∣⋍↻∽⋮∾≣↻↔↻
−∾⊨∘⊔∾≔↼↔⋎≝∇↬∘⋋↔≣∾↬ ↻≘↬∘↼∘↶⊟∘⊟⊔∽⋮≗⊢ ⋅↤⋎≗∽⊶⊔∾∽∘⋣∘⊔∽⊶≣⋯⋮⋮∘∘↙ ↶↔≣∣↬⋮↻≘⊔∾≘≣∃↔⋮∘∘∘≘ ∘∶∘∾⋅⇁≣∽⊽∘≌∽⊶∽↻∃∾⋯↦∾⊳ ∽∽∘⊨↔∘ ∾⋅−∾
↬≔∷⋅⊟≣⊔∽∘⊔ oy∥≣ . uorsuedxo juswoserdsip ones . juouIse . penbiq. e Sursn −∾⋅⊶∘↼∾−≣−⋯∶∾⋮⋯↩⊾∘∘≓
z¢d + AWd
+ x\d +d
=d
\\
\\
Add
+x\d
+0d
=d
⋋≋⊹⋉↽⊾⊹∘⊾∥⊾
—_
Nm
£/6 ‘«@ut .ll
↽↝∖∖⋅∾⊶⋂↔⋅⋔∙−↝
cm
1/9 z ut
/o1 ‘g w
'd — i
¥/11 :q¢ wt €/L:gzut
‘€1°S PUB ‘T[S ‘'S ‘¥'S ‘C1'¥ 31 398 ‘suonouny uonejodiaw @ 104,

7Y L7 Yy
umouys
o.hum ⋅⋅ ∙⋯ ⊢∥⊶⊶∘⋈−⇈↽∽ ↼ ∞↭↽∘⊓⇅↼⇀↽−∙⊪↖−
∼≊⊹⋋≋⊹⋉↽⊾⊹∘⊾∥⊾
umoys eie Sepou 8qISIA ∶∘∶∞∾⋅↼↩≻∙∶∘
293
Wwawae 4O 19ju80 u1 apoN
zzu ‘SJUSWISID [eldre[upenb
Jo Apurey (9 — ‘G om
Jo 1equiomt 351y ay7, (V]
speway
10jAe], ") PuUE POOH °d 99§
‘
umoys
ale sapou aqisiA
uou g ay} Ajug woudp J¢
umoys ale sapou /Z [e10} 8yl
10 61 [qisia ayy AjuO
sjutod fepoN

WOWOP (7
:az ut
28/LT ‘A¢ W
/6
19— D
SIS
[y I1AVL
J(uoupuoo dnsfut
ay1 £f5ups SjUaWI]2 D pUD SNONULUOD 21D SIUWI]D UIIMIaq 2n5Sa4d pup SjuswaID]dsip) siuawala 9d [n aandaffs snouvp
294
SJUDWID INI Se
0) paLIgjol OsSTY "3jqqnq 2iqud ap £q patpuud
V)
14 srerwouAjod jo 2oeds ap saoup !
uniod ‘W pue ‘1zzaig J ‘powry "N 'd 99§
7 = U ‘sjuswop Jemueln
3o Aprurey '*d — “J Ap 3o roquow 51y AL, (V] 10[Ae] "D pue pooH ‘d 9°S
:gg ul
"€1°S PUR “[1°C ‘S°S ‘'S ‘€I'v "8I 998 ‘suonouny uonejodiaw
e
saqeneA ainssaid pue Juswadedsip y)im 9poN @
sajqeueA Juswedeidsip yum spoN
/S
>¢/p:gTwm 'qd  d
2$/01 :g¢ Wl o¢/9 gz w
dTg
) 104,
296
Formulation
of the Finite Element Method
Chap. 4
We may ask whether in practice it is really important to satisfy the infsup condition, that is, whether perhaps this condition is too strong and elements that do not satisfy it can still be used reliably. Our experience is that if the infsup condition is satisfied, the element will be, for the interpolations used, as effective as we can reasonably expect and in that
sense optimal. For example, the 9/3 element for plane strain analysis in Table 4.6 is based on a parabolic interpolation of displacements and a linear interpolation of pressure.The element does not lock, and the order of convergence of displacements is always o(h®), and
of stresses, o(h?), which is surely the best behavior we can obtain with the interpolations used. On the other hand, if the infsup condition is not satisfied, the element will not always display for all analysis problems (pertaining to the mathematical model considered) the convergence characteristics that we would expect and indeed require in practice. The element is therefore not robust and reliable. Since the infsup condition is of great fundamental importance, we present in the following section a derivation that although not mathematically complete does yield valuable insight. In this discussion we will also encounter and briefly exemplify the ellipticity condition. For a mathematically complete derivation of the ellipticity and infsup conditions and many more details, we refer the reader to the book by F. Brezzi and M. Fortin [A]. In the derivation in the next section we examine the problem of incompressible elasticity, but our considerations are also directly applicable to the problem of incompressible fluid flow, and as shown in Section 4.5.7, to the formulations of structural elements. 4.4.4
Exercises
4.33. Use the fournode and eightnode shell elements available in a finite element program and perform the patch tests in Fig. 4.17. 4.34. Consider the threedimensional eightnode element shown. Design the patch test and identify analytically whether it is passed for the element.
8
u=Y hu;+ o6y + axpy + ozps
i=1 8
v=Y hyv;+ oy
[l
wl=
~
w=y
I NAw
i=
1
+ osps + asdy
hw; + o101 + agpa + 0igds (1 + 221
+ yyM1 + z2)
Pr1=1xLga=1y% ¢3=12
Displacement interpolation functions
Sec. 4.4
Incompatible and Mixed Finite Element Models
297
4.35. Consider the HuWashizu functional IIuw in (4.114) and derive in detail the equations (4.116) to (4.121).
4.36. The following functional is referred to as the HellingerReissner functional'’ 1
Iyr(u,1) = ∫
—21'TC"1' av + ∫ ↑⊤∂⋮∣∥↿∨ ⋁
⋁
−∫
qu”dV—
∫ w7t ds ~ ∫ Su
Sy
V
fSuT(uSu—u,,)dS
where the prescribed (not to be varied) quantities are % in V, u, on S,, and f% on S;. Derive this functional from the HuWashizu functional by imposing € = C™'t. Then
invoke the stationarity of Ilyr and establish all remaining differential conditions for the volume and surface of the body.
4.37. Consider the functional I,
=1

∫
fS"T(llS"—
su
llp)ds
whereHisgivenin(4.109)anduparethedisplacementstobeprescribedonthesurfaceS,,. Hence,thevectorfsurepresentstheLagrangemultipliers(surfacetractions)usedtoenforcethe surfacedisplacementconditions.Invokethestationarityofl'llandshowthattheLagrange multipliertermwillenforcethedisplacementboundaryconditionsonSu. 4.38. Consider the threenode truss element in Fig. E4.29. Use the HuWashizu variational principle and establish the stiffness matrices for the following assumptions: (a)
Parabolic displacement, linear strain, and constant stress
(b) Parabolic displacement, constant strain, and constant stress
Discuss your results in terms of whether the choices of interpolations are sensible (see Exam
ple 4.29). 4.39. Show that the following stressstrain expressions of an isotropic material are equivalent. Ty
=
Tij
=
Kev6,»,
+
2GE.’,
(a)
ijrserx
(b)
T = Ce
©
where « is the bulk modulus, G is the shear modulus,
“T
E
E
30  29
2(1 + »)
E is Young’s modulus, » is Poisson’s ratio, ey is the volumetric strain, and €}; are the deviatoric strain components, €v = €u; AlSO,
Cijr:
=
/\60 6rs
€ =€ +
M«(airajs
€y
— =& 3
+
61':6jr)
7This functional is sometimes given in a different form by applying the divergence theorem to the second term.
298
Formulation of the Finite Element Method
Chap. 4
where A and p are the Lamé constants,
A

B U+ na 2
I K20 +»)
In (a) and (b) tensorial quantities are used, whereas in (c) the vector of strains contains the engineering shear strains (which are equal to twice the tensor components; €.g., ¥,y = €12 + €). Also, the stressstrain matrix C in (c) is given in Table 4.3.
4.40. Identify the order of pressure interpolation that should be used in the u/p formulation in order to obtain the same stiffness matrix as in the pure displacement formulation. Consider the following elements of 2 X 2 geometry. (a) Fournode element in plane strain
(b) Fournode element in axisymmetric conditions (c) Ninenode element in plane strain.
441. Consider the 4/1 element in Example 4.32 and assume that the displacement boundary condition to be imposed is u; = u#. Show formally that imposing this boundary condition prior to or after the static condensation of the pressure degree of freedom, yields the same element contribution to the stiffness matrix of the assemblage. 4.42. Consider the axisymmetric 4/1 u/p element shown. Construct the matrices B, By, C’, and H, for this element.
¢ 1
3
3


—2—
4.43. Consider the 4/3c element in plane strain conditions shown. Formulate all displacement and strain interpolation matrices for this element (see Table 4.7).
Sec. 4.4
Incompatible and Mixed Finite Element Models
299
4.44. Consider the 9/3 plane strain u/p element shown. Calculate the matrix Kp,.
⊺⊢−−−−−⊣
y
2
X
⊥⊡
Young's modulus £ Poisson's ratio v = 0.49
4.45. Consider the plate with the circular hole shown. Use a finite element program to solve for the stress distribution along section AA for the two cases of Poisson’s ratios » = 0.3 and » = 0.499. Assess the accuracy of your results by means of an error measure. (Hint: For the analysis with v = 0.499, the 9/3 element is effective.)
Plane strain condition Young's modulus £ = 200,000 MPa
300
Formulation of the Finite Element Method
Chap. 4
4.46. The static response of the thick cylinder shown is to be calculated with a finite element program.
f = force per unit length
30 mm
E = 200,000 MPa v=0.499
Use idealizations based on the following elements to analyze the cylinder. (a) Fournode displacementbased element (b) Ninenode displacementbased element
(¢) 4/1 u/p element. (d) 9/3 u/p element. In each case use a sequence of meshes and identify the convergence rate of the strain energy.
4.5 THE INFSUP CONDITION FOR ANALYSIS OF INCOMPRESSIBLE MEDIA AND STRUCTURAL PROBLEMS As we pointed out in the previous section, it is important that the finite element discretization for the analysis of almost, and of course totally, incompressible media satisfy the infsup condition. The objective in this section is to present this condition. We first consider the pure displacement formulation for the analysis of solids and then the displacement/pres
Sec. 4.5
The InfSup Condition
301
sure formulations. Finally, we also briefly discuss the infsup condition as applicable to structural elements. In our discussion we apply the displacement and displacement/pressure formulations to a solid medium. However, the basic observations and conclusions are also directly applicable to the solution of incompressible fluid flows if velocities are used instead of displacements (see Section 7.4). 4.5.1
The InfSup Condition Derived from Convergence Cpnsiderations
We want to solve a general linear elasticity problem (see Section 4.2.1) in which a body is
subjected to body forces f?, surface tractions f% on the surface S;, and displacement boundary conditions u’« on the surface S.. Without loss of generality we want to reach in this section, we can assume that the prescribed prescribed tractions £ are zero. Of course, we assume that the body so that no rigid body motions are possible. We can then write our
of the conclusions that displacements u®+ and is properly supported, analysis problem as a
problem of minimization, 1
min {Ea(v, v) + ⋮∫ vevy
(divv)2Wol—
2
Vol
∫
f8.v dVol}
(4.151)
Vol
where using indicial notation and tensor quantities (see Sections 4.3.4 and 4.4.3), a(u, v) = 2G
∑ ∊⋅⋅∣∫⊳≺∐⋝∊∣⋅∣∫⊳≼⋁⊃∠∣⋁∘↕ Voli,j
, » , 8 u v i d § ~ — ) u ( » , ( e = ) u ( , e{ ..
€;(u)
1 6u,» 6uj — — + — — =—
2(
%,
axi)
where x = E/[3(1 — 2»)] (bulk modulus), Young’s modulus, v = Poisson’s ratio.
5
(4.152) di
.
vV == 0o
G = E/[2(1 + v)] (shear modulus),
E =
V= {vlg—”)é € L(Vol), i,j = 1,2, 3;vils, = 0,i= 1, 2, 3} In these expressions we use the notation defined earlier (see Section 4.3) and we denote by
“Vol” the domain over which we integrate so as to avoid any confusion with the vector space V. Also, we use for the vector v and scalar g the norms
Ivlp =2i
ax;
;
I g8 = 1l g llZzven
(4.153)
L2(Vol)
where the vector norm   [lv is somewhat easier to work with but is equivalent to the Sobolev
norm    defined in (4.76) (by the PoincaréFriedrichs inequality).
302
Formulation of the Finite Element Method
Chap. 4
In the following discussion we will not explicitly give the subscripts on the norms but
always imply that a vector w has norm  wly and a scalar y has norm  ylo. Let u be the minimizer of (4.151) (i.e., the exact solution to the problem) and let V
be a space of a sequence of finite element spaces that we choose to solve the problem. These spaces are defined in (4.84). Of course, each discrete problem, 1
lim {—a(v,,, v,) + £
wevs
9
(2
↾∄⋅⋁⇂∎ dVol}
(div v,)* d Vol — ∫ Vol
(4.154)
Vol
has a unique finite element solution u,. We considered the properties of this solution in Section 4.3.4, and in particular we presented the properties (4.95) and (4.101). However, we also stated that the constants ¢ in these relations are dependent on the material properties. The important point now is that when the bulk modulus « is very large, the relations (4.95) and (4.101) are no longer useful because the constants are too large. Therefore, we
want our finite element space V, to satisfy another property, still of the form (4.95) but in which the constant ¢, in addition to being independent of h, is also independent of «. To state this new desired property, let us first define the “distance” between the exact solution u and the finite element space V; (see Fig. 4.22),
d(u, Vi) =
inf lu — vaf = lu —
(4.155)
VAEV,
where 1, is an element in V) but is in general not the finite element solution.
The Basic Requirements In engineering practice, the bulk modulus k may vary from values of the order of G to very large values, and indeed to infinity when complete incompressibility is considered. Our objective is to use finite elements that are uniformly effective irrespective of what value luupll
diu, V) = llu=d,ll
Figure 4.22
Schematic representation of solutions and distances; for optimal convergence
lu = w,) = c d(u, V,) with ¢ independent of h and .
Sec. 45
303
The InfSup Condition
takes. Mathematically, therefore, our purpose is to find conditions on V, such that
ull

uh”
=c
d(ll,
Vh)
with a constant ¢ independent of h and k.
156)
(4
’
These conditions shall guide us in our choice of effective finite elements and discretizations. The inequality (4.156) means that the distance between the continuous solution u and the finite element solution u, will be smaller than a (reasonably sized) constant ¢ times d(u, V3;) and that this relationship will be satisfied with the same constant ¢ irrespective of the bulk modulus used. Note that this independence of ¢ from the bulk modulus is the key property we did not have in Section 4.3.4 when we derived a relation such as (4.156) [see (4.95)]. Assume that the condition (4.156) holds (with a reasonably sized constant ¢). Then
if d(u, Vi) is o(h*¥), we know that ju — w, is also o(h*%), and since c is reasonably sized and
independent of &, we will in fact observe the same solution accuracy and improvement in accuracy as h is decreased irrespective of the bulk modulus in the problem. In this case the finite element spaces have good approximation properties for any value of k, and the finite element discretization is reliable (see Section 1.3). The relationship in (4.156) expresses our fundamental requirement for the finite element discretization, and finite element formulations that satisfy (4.156) do not lock (see Section 4.4.3). In the following discussion, we write (4.156) only in forms with which we can work more easily in choosing effective finite elements. One of these forms uses an infsup value and is the celebrated infsup condition. To proceed further, we define the spaces K and D,
Kig) ={vvE V,divv =g} D = {q q = div v for some v € V}
(4.157) (4.158)
and the corresponding spaces for our discretizations,
Ki(gn) = {vu  vi € Vi, div vi = gu} Dy = {gs  g» = div vs for some v, € Vi}
Hence the space K(gs), for div v, = gs. Also, the space reached by the elements v, element v, in V, such that and D.
(4.159) (4.160)
a given g,, corresponds to all the elements v in V;, that satisfy D, corresponds to all the elements g, with g, = div v, that are in Vi; that is, for any g, an element of D, there is at least one g, = div v,. Similar thoughts are applicable to the spaces K
We recall that when « is large, the quantity  div u, will be small; the larger k, the smaller  div ws, and it is difficult to obtain an accurate pressure prediction ps = ~K div u,. In the limit k — % we have div u, = 0, but the pressure p, is still finite (and
of course of order of the applied tractions) and therefore k(div u,)* = 0.
304
Formulation of the Finite Element Method
Chap. 4
Before developing the infsup condition, let us state the ellipticity condition for the problem of total incompressibility: there is a constant a greater than zero and independent of h such that a(Ve, V)
=
Vv,
“V;. "2
€
(4161)
Kh(O)
This condition in essence states that the deviatoric strain energy is to be bounded from below, a condition that is clearly satisfied. We further refer to and explain the ellipticity condition for the incompressible elasticity problem in Section 4.5.2. Let us emphasize that in this finite element formulation the only variables are the displacements.
Obtaining the InfSup Condition The
infsup condition—which
when
satisfied ensures that (4.156)
holds—can
now be
developed as follows. Since the condition of total incompressibility clearly represents the most severe constraint, we consider this case. Theng = 0, u belongs to K(g) for ¢ = 0 [that is, K(0)], and the continuous problem (4.151) becomes min
vek(©o)
f”v dVol}
{la(v, v) — ∫ 2
(4.162)
with the solution u, while the discrete problem is min
wek,y(0)
f”v;. dVol}
{la(vh, Vi) — ∫ 2
(4.163)
Vol
with the solution u,. Now consider condition (4.156). We notice that in this condition we compare distances. In the following discussion we characterize a distance as “small” if it remains of the
same order of magnitude as d(u, V,) as h decreases. Similarly, we will say that a vector is small if its length satisfies this definition and that a vector is “close” to another vector if the vector difference in the two vectors is small.
Sinceun, € K,(0), and therefore always ju — w, = ¢ d[u, K4(0)] (see Exercise 4.47),
we can also write condition (4.156) in the form
du, K.(0)] = c d(u, V,)
(4.164)
which means that we want the distance from u to X,,(0) to be small. This relation expresses the requirement that if the distance between u and V, (the complete finite element displacement space) decreases at a certain rate as h — 0, then the distance between u and the space in which the actual solution lies [because u, € K,(0)] decreases at the same rate. Figure 4.23 shows schematically the spaces and vectors that we use. Let u,o be a vector of our choice in K,(0) and let w, be the corresponding vector such that ﬁh
=
Uno
+
W,
(4165)
Sec. 45
305
The InfSup Condition
diu, Kx{0)]
Figure 4.23
Spaces and vectors considered in deriving the infsup condition
We can then prove that the condition in (4.164) is fulfilled provided that for all g, € Dy, there is a w, € Ki(gs) such that
(4.166)
wall = ¢ [ gnl where ¢’ is independent of h and the bulk modulus «.
First, we always have (see Exercise 4.48)
diviu — @)  = alju — @, and hence,
[ div @,
(4.167)
= a d(u, Vi)
(4.168)
where « is a constant and we used divua = 0. Second, we consider
lw = wol = Ju — @, + w
< u — Now
assume
that
(4.166)
holds
with
 +  wll
g, = div @i,..
Because
div u, = 0,
div @, = div w,, where we note that @i, is fixed by (4.155) and therefore g
we
have
is fixed, but by
choosing different values of u, different values of w, are also obtained. Then it follows that
lu—uwel
inf ‘Ihelgl(bh) viléeh
∥ Vi
B>0
(4.183)
∥ ∥ q" ∥
In other words, the infsup condition now corresponds to any element in V, and Pu(Dy).
Hence, when applying (4.166), (4.175), or (4.176) to the mixed interpolated u/p elements, we now need to consider the finite element spaces Vi, and P.(D»), where Py(Ds) is used instead of Dj.
EXAMPLE 4.36:
Prove that the infsup condition is satisfied by the 9/3 twodimensional «/p
element presented in Section 4.4.3.
For this proof we use the form of the infsup condition given in (4.176) (see F. Brezzi and K. J. Bathe [A]). Given u smooth we must find an interpolation, u, € V4, such that for each element m,
∫
. (div u — div u)g, dVol™ = 0
(a)
Volt™,
for all g polynomials of degree = Steens Sa
(u—w) nds
(d)
S;
We are left to use the two degrees of freedom at the element center node. We choose these in such a way that
∫
div(u—u,)deol"”)= ∫ Vol()
div(u — u)y dVol™ = 0
©
Vol(m)
We note now that (d) and (e) imply (a) and that u,, constructed element by element through (b) and (c), will be continuous from element to element. Finally, note that clearly if u is a (vector)
polynomial of degree 0
(4.204)
P, (div w,,) P,{div v,) dVol = ∫
(4.205)
Vol
P;.(divw,.)divv;.dVol Vol
The relation (4.204) is in matrix form
.
W{G,.V;.
(4.206)
0 7 P S 2 a v , ot sUp vTGiWs
inf sup
where W, and V; are vectors of the nodal displacement values corresponding to w; and vy,
and Gy, Sy are matrices corresponding to the operator b’ and the norm + lv, respectively. The matrices G, and S, are, respectively, positive semidefinite and positive definite (for the problem we consider, see Section 4.5.1). EXAMPLE 4.40: 1n Example 4.34 we calculated the matrix G, of a 4/1 element. Now also establish the matrix S, of this element.
To evaluate S, we recall that the norm of w is given by [see (4.153)]
pwip = =  22
aw,»
iJj
a.X.'j
12(Vol)
Hence, for our case
= [ LGS+ G+ (G = e where u, v are the components w;, i = 1, 2. Let us order the nodal point displacements in @ as in Example 4.34, @' =l
w
wus
ug
¥
o1
vy
v3;
?
= ∑ hi.yui
v4
By definition,  w, ]} = &’S,d. Also, we have du_
∑ ⇂↧∣⊳⋅≍∥∣⊳⊑
(b)
33)
and we write in (a)
©
Substituting from (c) and (b) into (a) we obtain
s 0= [ [ o+ o1 aray +1
Sh(l,
2)
=
∫
+1
∫ 1
2 =3 1
∣⋮∣↧↥⋅≍∣↧≵⋅≍
+
1
and so on,
hl_yhz_y]
dx
dy
=
_8
324
Formulation of the Finite Element Method
Chap. 4
Similarly, the terms corresponding to the v; degrees of freedom are calculated, and we obtain
4
§h]’
S"'[o
1 2
$=6l2
1
4
2 1
1
2
1
4
c
071
_S.
1
_1]1
4
1
Let us now consider the u/pc formulation. In this case the same expression as in (4.206) applies, but we need to use G» = (K,.)7T# '(K,.)4, where Tj, is the matrix of the L* norm of p, (see Exercise 4.59); i.e., for any vector of pressure nodal values P, we have
 pall = PET,Ps.
The form (4.206) of the infsup condition is effective because we can numerically evaluate the infsup value of the lefthand side and do so for a sequence of meshes. If the lefthandside infsup value approaches (asymptotically) a value greater than zero (and there are no spurious pressure modes, further discussed below), the infsup condition is satisfied. In practice, only a sequence of about three meshes needs to be considered (see examples given below). The key is the evaluation of the infsup value of the expression in (4.206). We can show that this value is given by the square root of the smallest nonzero eigenvalue of the problem Ghd)h
=
/\S}.d)}.
(4207)
Hence, if there are (k — 1) zero eigenvalues (because G is a positive semidefinite matrix)
and we order the eigenvalues in ascending order, we find that the infsup value of the expression in (4.206) is V. We prove this result in the following example. EXAMPLE 4.41:
Consider the function f(U, V) defined as
U'GV f@G, V) = (U"GU)VA(VTSV) /2
(@
where G is an n X n symmetric positive semidefinite matrix, S is an n X n positive definite matrix, and U, V are vectors of order n. Show that
inf sup f(U, V) = VA,
(®)
where A is the smallest nonzero eigenvalue of the problem
Gd
= ASd
(e
=0
dA — ∫ wpdA+boundaryterms(a)
= ∫
A
A
where k, Cs, 7, C, have been defined in (5.95) to (5.97) and y*® contains the assumed transverse shear strains ,YAS
oS = [ "s] = constant ¥z
The relation in (a) is a modified HellingerReissner functional. Substituting the interpolations for w, Bx, and B, into « and v, integrating over the element midsurface area A, and invoking the
stationarity of IIjiz with respect to the nodal point variables d,
478
Formulation and Calculation of Isoparametric Finite Element Matrices
where
Kb
=
B[T; CbBb
Chap. 5
dA
A
D=J‘C:dA=ACx A
G=C
∫ B, dA A
and B, and B, are straindisplacement matrices,
k = B,
v = B,a Using static condensation, we obtain the stiffness matrix of the element with respect to the nodal point variables only,
K=K,
+GD'G
As we discussed in Section 5.4.2, the pure displacementbased isoparametric plate element (i.e., using full numerical integration for the bending and transverse shear terms in the displacementbased stiffness matrix) is much too stiff (displays the shear locking phenomenon). The presentation in Example 5.44 shows that the onepoint integrated element has a variational basis quite analogous to the basis of the onepoint integrated isoparametric beam element. However, whereas the beam element is reliable and effective, the plate element stiffness matrix in Example 5.44 has a spurious zero eigenvalue and hence the element is unreliable and should not be used in practice (as was pointed out by J.L. Batoz, K. J. Bathe, and L. W. Ho [A]). The important point of this example is that a variational basis of an element might well exist, but whether the element is useful and effective can of course be determined only by a deeper analysis of the formulation. The equivalence between a certain isoparametric reduced or selectively integrated displacementbased element and a mixed formulation may also hold only for specific geometric element shapes and may also no longer be valid when nonisotropic material laws (or geometric nonlinearities) are introduced. An analysis of the effects of each of these conditions should then be performed. 5.5.7
Exercises
5.53. Evaluate the NewtonCotes constants when the interpolating polynomial is of order 3, i.e., ¥{r) is a cubic. 5.54. Derive the sampling points and weights for threepoint Gauss integration.
Sec. 5.5
479
Numerical Integration
5.55. Show that 3 X 3 Gauss numerical integration is sufficient to calculate the stiffness and mass matrices of a ninenode geometrically undistorted displacementbased element for axisymmetric analysis.
5.56. Show that 2 X 2 Gauss integration of the stiffness matrix of the eightnode plane stress displacementbased square element results in the spurious zero energy mode shown. (Hint: You
need to show that Bii = 0 for the given displacements.) u 2u
2u
(

Symmetric

deformations
5.57. Consider the 9/3 u/p element and show that 3 X 3 Gauss integration of a geometrically undistorted element gives the exact stiffness matrix. Also, show that 2 X 2 Gauss integration is not adequate. 5.58. Identify the requiréa integration order for full integration of the stiffness matrix of the sixnode displacementbased triangular element when using the Gauss integration in Table 5.8.
Plane stress element
5.59. Consider the ninenode plane stress element shown. All nodal point displacements are fixed except that u, is free. Calculate the displacement u; due to the load P. (a) Use analytical integration to evaluate the stiffness coefficient. (b) Use 1 X1, 2 X2, and 3 X 3 Gauss numerical integration
to evaluate
the
stiffness
coefficient. Compare your results.
E = 200,000 v=0.30
Thickness t=1.0
5.60. Consider the evaluation of lumped mass matrices for the elements shown in Table 5.9. Determine suitable Gauss integration orders for the evaluation of these matrices.
480
Formulation and Calculation of Isoparametric Finite Element Matrices
Chap. 5
5.61, Consider the plate bending element formulation in Example 5.29. Assume that the element stiffness matrix is evaluated with onepoint Gauss integration. Show that this element has
spurious zero energy modes."? 5.62. Consider the plate bending element formulation in Example 5.29 and assume that the bending strain energy is evaluated with 2 X 2 Gauss integration and the shear strain energy is evaluated with onepoint Gauss integration. Show that this element has spurious zero energy modes."
5.6 COMPUTER PROGRAM FINITE ELEMENTS
IMPLEMENTATION OF ISOPARAMETRIC
In Section 5.3 we discussed the isoparametric finite element formulation and gave the specific expressions needed in the calculation of fournode plane stress (or plane strain) elements (see Example 5.5). An important advantage of isoparametric element evaluations is the similarity between the calculations of different elements. For example, the calculation of threedimensional elements is a relatively simple extension from the calculation of twodimensional elements. Also, in one subroutine, elements with a variety of nodal point configurations can be calculated if an algorithm for selecting the appropriate interpolation functions is used (see Section 5.3). The purpose of this section is to provide an actual computer program for the calculation of the stiffness matrix of fournode isoparametric elements. In essence, SUBROUTINE QUADS is the computer program implementation of the procedures presented in Example 5.5. In addition to plane stress and plane strain, axisymmetric conditions can be considered. It is believed that by showing the actual program implementation of the element, the relative ease of implementing isoparametric elements is best demonstrated. The input and output variables and the flow of the program are described by means of comment lines.
SUBROUTINE
c C. C . C . C . C . C .
QUADS
QUA00001
(NEL,ITYPE,NINT,THIC,YM,PR,XX,S,I0UT)
PROGRAM TO CALCULATE ISOPARAMETRIC QUADRILATERAL ELEMENT PLANE STRESS, AND PLANE MATRIX FOR AXISYMMETRIC, CONDITIONS
QUA00002
STIFFNESS STRAIN
c . C C Cc Cc C Cc C Cc Cc C Cc C C Cc
. . . . . . . . . . . . .


INPUT VARIABLES  = NUMBER OF ELEMENT NEL = ELEMENT TYPE ITYPE EQ.0 = AXISYMMETRIC EQ.l = PLANE STRAIN EQ.2 = PLANE STRESS = GAUSS NUMERICAL INTEGRATION ORDER NINT = THICKNESS OF ELEMENT THIC = YOUNG'’S MODULUS YM = POISSON’S RATIO PR = ELEMENT NODE COORDINATES Xx(2,4) = STORAGE FOR STIFFNESS MATRIX 5(8,8) = UNIT NUMBER USED FOR OUTPUT I0UT
.
12 Note that these elements should therefore not be used in practice (see Section 5.5.5).
. . . . . .
QUA00003 QUA00004 QuUA00005 QUA00006 QUA00007 QUA00008
. . . . . . . . . . . . . .
QUA00010 QUA00O11 QUA00012 QUA00013 Quag0014 QUA00015 QUA00016 QUA00017 QUA00018 QUA00019 QUA00020 QUA00021 QUA00022 QUA00023
. QUA00009
Computer Program Implementation of Isoparametric Finite Elements

a0
oo
aoanooaao
Sec. 5.6
1 2 [sNeXe]
3
a0an
1 2 3
=
OUTPUT s(8,8)

. . . .
QUA00024 = CALCULATED STIFFNESS MATRIX QUA0002S QUA00026 QUA00027 IMPLICIT DOUBLE PRECISION (AH,02) QUA00028 . QUA00029 THIS PROGRAM Ié USED IN SINGLE PﬁEéIéIéN ARITHMETIC 6N CRAY . QUA00030 EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM MACHINES, . QUA00031 ENGINEERING WORKSTATIONS AND PCS. DEACTIVATE ABOVE LINE FOR QUA00032 SINGLE PRECISION ARITHMETIC. . QUA00033 . QUA00034 DIMENSION D(4,4),B(4,8),Xxx(2,4),5(8,8),XG(4,4),WGT(4,4),DB(4) QUA00035 QUA00036 MATRIX XG STORES GAUSS  LEGENDRE SAMPLING POINTS QUA00037 QUA00038 DATA XG/ 0.DO, 0.p0, 0.D0, 0.00, .5773502691896D0, QUA00039 .5773502691896D0, 0.00, 0.D0, .7745966692415D0, 0.D0, QUA00040 .774596669241500, 0.D00, .8611363115941D0, QuUA00041 .3399810435849D0, .3399810435849D0, .8611363115941p0 / QUA00042 QUA00043 MATRIX WGT STORES GAUSS  LEGENDRE WEIGHTING FACTORS QUA00044 QUA00045 DATA WGT / 2.DO, 0.D0, 0.p0, 0.00, 1.p0, 1.p00, QUA00046 0.00, 0.00, .5555555555556D0, .8888888888889D0, QuUA00047 .5555555555556D0, 0.00, .3478548451375D0, .6521451548625D0,QUA00048 .6521451548625D0,
OBTAIN
.347854845137500
STRESSSTRAIN
/
QuUA00049 QUA00050
LAW
QUA00051 QUA00052 QuUA00053
F=¥YM/(1.+PR) Gm=F*PR/(1.2.*PR)
ano
H=F
+
PLANE
QuUA00054 QUA00055
G
STRAIN
QUA00056 QuUA00057 QUADD058 QUA00059 QUA00060 QUA00061 QUA00062
ANALYSIS
D(1,1)=H D(1,2)=G D(1,3)=0.
D(2,1)=G D(2,2)=H D(2,3)=0.
QUA00063
QUA00064 QUA00065 QUADQG066 QuUA00067 QUA00068 QuUA00069 QUA00070
D(3,1)=0.
D(3,2)=0, D(3,3)=F/2.
ano
IF (ITYPE.EQ.1) THIC=1. GO TO 20 ENDIF AXISYMMETRIC
[sNeXe]
FOR
PLANE
THEN
Qua00071
ANALYSIS
D(1,4)=G D(2,4)=G D(3,4)=0. D(4,1)=G D(4,2)=G D(4,3)=0. D(4,4)=H IF (ITYPE.EQ.0)
c
481
STRESS
GO
TO
20
ANALYSIS
DO 10 I=1,3 A=D(1,4)/D(4,4) DO 10 J=I,3 D(I,J)=D(I,J)  D(4,J)*A 10 p(J,I)=D(I1,J)
CONDENSE
STRESSSTRAIN
MATRIX
QUA00072 QUA00073 QUAQ00Q074 QUA00075 QUA00076 Qua00077 QUA00078 QUA00079 QUA00080 QUA00081 QUA00082 QUA00083 QuUA00084 QUA00085S QUA00086 QUA00087 QUA00088 QUA00089 QUA00090 QUA00091
482
Formulation
and Calculation of Isoparametric Finite Element Matrices
QUA00092
STIFFNESS
ELEMENT
CALCULATE
C
Chap. 5
QUA00093 QUA00094 QUA00095 QUAD0096
C 20 30
c C
PO 30 I=1,8 DO 30 J=1,8 s(1,J)=0.
IST=3 IF (ITYPE.EQ.0) DO 80 LX=1,NINT RI=XG(LX,NINT) DO 80 Ly=1,NINT SI=XG(LY,NINT)
IST=4
DET
DETERMINANT
JACOBIAN
THE
B AND
OPERATOR
DERIVATIVE
EVALUATE
QUA00097 QUAD0098 QUA00099 QUA00100 QUA00101 QUA00102 QUA00103 QUAOOiM
QUA00105
C
QUA00106
(XX,B,DET,RI,SI,XBAR,NEL,ITYPE,IOUT)
STDM
CALL
QUADD107
C
QUAD0D108
STIFFNESS
ELEMENT
TO
CONTRIBUTION
ADD
c
QUAO00109
c
QUAO00110 QUA00111
IF (ITYPE.GT.0) XBAR=THIC WP=WGT( LX,NINT) *WGT (LY, NINT) *XBAR*DET
QUA00112 QUAO00113 QUADO114
DO 70 J=1,8 DO 40 K=1,IST DB(K)=0.0 40
50 60 70 80
DB(K)=DB(K)
+
QUADO116
D(K,L)*B(L,J)
QUA00117
I=J,8
60
DO
QUADD115
L=1,IST
40
DO
STIFF=0.0 DO 50 L=1,IST
QUA00118 QUAO00119
STIFF=STIFF + B(L,I)*DB(L) S(I,J)=S(I,J) + STIFF*WT CONTINUE CONTINUE
QUA00120 QUA00121 QUA00122 QUA00123
DO 90 J=1,8 DO 90 I=J,8 S(J,I)=S(I1,J)
QUAD0124 QUAD0125 QUAD0126 QUA00127
C 90
c
QUA00128 QUAD0129
RETURN
QUA00130 QUA00131
C END C C . (o
e .
PROGRAM TO EVALUATE THE STRAINDISPLACEMENT TRANSFORMATION AT POINT (R,S) FOR A QUADRILATERAL ELEMENT
C . c . c . C C
QUA00132
(XX,B,DET,R,S,XBAR,NEL,ITYPE,IOUT)
STDM
SUBROUTINE
. et
e
e
e
e
e
IMPLICIT
e
e
DOUBLE
e
e
e
e
e
PRECISION
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
(AH,02)
C
C
c C
C
C C
C
.
.
. QUA00139 . Qua0D140
QUADO141
QUA00142 QUA00143 QUADD144
XX(2,4),B(4,8),H(4),P(2,4),XJ3(2,2),XJ1(2,2)
DIMENSION
. QUADO0136 . QUA00137 . QUA00138
B
MATRIX
QUA0D0133 o117 X113 & 7 QUAO00135
R
RP
=
1.0
+
SP RM
= =
1.0 1.0
+ S  R
QUADD145 QUA00146
SM
=1.0
8
QUA00147 QUA00148
INTERPOLATION
QUAD0149%
FUNCTIONS
QUAO00150
H(l)
=
0.25%
RPx
SP
QUAD0151
H(2) H(3) H(4)
= = =
0.25* 0.25* 0.25*
RM* RM* RP*
SP SM SM
QUAO00152 QUA00153 QUA00154 QUA00155
COORDINATE
NATURAL
DERIVATIVES
OF
THE
INTERPOLATION
FUNCTIONS
QUA00155>$
QUA001
1.
WITH
RESPECT
C P(1,1) P(1,2) P(1,3) P(1,4)
= = = =
0.25* sp  p(1,1)  0.25*% SM  P(1,3)
TO
R
QUA00158 QUA0015% QUA00160 QUA00161 QUA00162 QUA00163
Sec. 5.6
c c c
c c c
2. WITH RESPECT TO S P(2,1) = 0,25* RP P(2,2) = 0.25% RM P(2,3) =  P(2,2) P(2,4) =  P(2,1) EVALUATE THE JACOBIAN MATRIX AT POINT 10 DO 30 I=1,2 DO
c C c
c c c
c c c
c c c c c c c c c c c c c c c c
Computer Program Implementation of Isoparametric Finite Elements
30
(R,S)
J=1,2
DUM = 0.0 DO 20 K=1,4 20 DUM=DUM + P(I,K)*XX(J,K) 30 XJ(I,J)=DUM COMPUTE THE DETERMINANT OF THE JACOBIAN MATRIX AT POINT (R,S) DET = XJ(1,1)* XJ(2,2)  XJ(2,1}* XJ(1,2) IF (DET.GT.0.00000001) GO TO 40 WRITE (IOUT,2000) NEL GO TO 800 COMPUTE INVERSE OF THE JACOBIAN MATRIX 40 DUM=1./DET XJI(1,1) = XJ(2,2)* DUM XJI(1,2) =XJ(1,2)* DUM XJI(2,1) =XJ(2,1)* DUM XJ1(2,2) = XJ(1,1)* DUM EVALUATE GLOBAL DERIVATIVE OPERATOR B K2=0 DO 60 K=1,4 K2=K2 + 2 B(1,Kk21) = 0. B(1,K2 )} = 0. B(2,K21) = 0. B(2,K2 ) = 0. DO 50 I=1,2 B(l,K21) = B(1,K21) + XJI(1,I) * P(I,K) 50 B(2,K2 ) = B(2,K2 ) + XJI(2,I) * P(I,K) B(3,K2 ) = B(1,K21) 60 B(3,K21) = B(2,K2 ) IN CASE OF PLANE STRAIN OR PLANE STRESS ANALYSIS DO NOT INCLUDE THE NORMAL STRAIN COMPONENT IF (ITYPE.GT.0) GO TO 900 COMPUTE THE RADIUS AT POINT (R,S) XBAR=0.0 DO 70 K=1,4 70 XBAR=XBAR + H(K)*XX(1,K) EVALUATE THE HOOP STRAINDISPLACEMENT RELATION IF (XBAR.GT.0.00000001) GO TO 90 FOR THE CASE OF ZERO RADIUS EQUATE RADIAL TO HOOP STRAIN DO 80 K=1,8 80 B(4,K)=B(1,K) GO TO 900 NON2ERO RADIUS
483
QUA00164 QUA00165 QUA00166 QUAD0167 QUA00168 QUAD0169 QUA00170 QUA00171 QUA00172 QUA00173 QUA00174
QUA00175
QUA00176 QUA00177 QUA00178 QUA00179 QUA00180 QUA00181 QUA00182 QUA00183 QUA00184 QUA00185 QUA00186 QUA00187 QUA00188 QUA00189 QUA00190 QUA00191 QUA00192 QUA00193 QUA00194 QUA00195 QUA00196 QUA00197 QUA00198 QUA00199 QUA00200 QUA00201 QUA00202 QUA00203 QUA00204 QUA00205 QUA00206 QUA00207 QUA00208 QUA00209 QUA00210 QUA(0211 QUA00212 QUA00213 QUA00214 QUA00215 QUA00216 QUA00217 QUA00218 QUAQ0219 QUA00220 QUA00221 QUA00222 QUA00223 QUA00224 QUA00225 QUA00226 QUA00227 QUA00228 QUA00229 QUA00230 QUA00231 QUA00232 QUA00233
484
Formulation and Calculation of Isoparametric Finite Element Matrices 90
100
DUM=1./XBAR K2=0 DO 100 k=1,4 K2=K2 + 2 B(4,K2 ) = 0. B(4,K21) = H(K)*DUM GO TO 900
C 800 900
sTOP RETURN
C 2000
FORMAT (//,' %*%* ERROR 1 ! ZERO OR NEGATIVE
C END
*** 7, JACOBIAN
DETERMINANT
FOR
ELEMENT
Chap. 5
QuAa00234 QUA00235 QuAa00236 QUA00237 Qua00238 Qua00239 QUA00240 Qua00241 Qua00242 Qua00243 QuUA00244 QUA00245 (’,IB,')')QUA88246 QUA00247 QUA00248
Il CHAPTER
SIX I
Finite Element
Nonlinear Analysis in Solid and Structural Mechanics
6.1 INTRODUCTION TO NONLINEAR ANALYSIS In the finite element formulation given in Section 4.2, we assumed that the displacements of the finite element assemblage are infinitesimally small and that the material is linearly elastic. In addition, we also assumed that the nature of the boundary conditions remains unchanged during the application of the loads on the finite element assemblage. With these assumptions, the finite element equilibrium equations derived were for static analysis
KU =R
(6.1)
These equations correspond to a linear analysis of a structural problem because the displacement response U is a linear function of the applied load vector R; i.e., if the loads are aR instead of R, where « is a constant, the corresponding displacements are aU. When this is not the case, we perform a nonlinear analysis. The linearity of a response prediction rests on the assumptions just stated, and it is instructive to identify in detail where these assumptions have entered the equilibrium equations in (6.1). The fact that the displacements
must be small has entered into the
evaluation of the matrix K and load vector R because all integrations have been performed over the original volume of the finite elements, and the straindisplacement matrix B of each element was assumed to be constant and independent of the element displacements. The assumption of a linear elastic material is implied in the use of a constant stressstrain matrix C, and, finally, the assumption that the boundary conditions remain unchanged is reflected in the use of constant constraint relations [see (4.43) to (4.46)] for the complete response.
If during loading a displacement boundary condition should change, e.g., a degree of freedom which was free becomes restrained at a certain load level, the response is linear only prior to the change in boundary condition. This situation arises, for example, in the analysis of a contact problem (see Example 6.2 and Section 6.7).
486
Finite Element Nonlinear Analysis in Solid and Structural Mechanics
Chap. 6
The above discussion of the basic assumptions used in a linear analysis defines what we mean by a nonlinear analysis and also suggests how to categorize different nonlinear analyses. Table 6.1 gives a classification that is used conveniently because it considers separately material nonlinear effects and kinematic nonlinear effects. The formulations listed in the table are those that we shall discuss in this chapter. TABLE 6.1
Classification of nonlinear analyses
Type of analysis
Materiallynonlinearonly
Description
Infinitesimal displacements and strains; the
Typical formulation used
Materiallynonlinear only (MNO)
Stress and strain measures
Engineering stress and strain
stressstrain relation is nonlinear Large displacements, large rotations, but small strains
Displacements and rotations of fibers are large, but fiber extensions and angle changes between fibers
Total Lagrangian (TL)
Updated Lagrangian (UL)
Second PiolaKirchhoff stress, GreenLagrange strain Cauchy stress, Almansi strain
are small; the stressstrain relation may be linear or nonlinear Large displacements, large rotations, and large strains
Fiber extensions and angle changes between fibers are large, fiber displacements and rotations may also be large; the stressstrain relation may be linear or nonlinear
Total Lagrangian (TL)
Updated Lagrangian (UL)
Second PiolaKirchhoff stress, GreenLagrange strain Cauchy stress, logarithmic strain
Figure 6.1 gives an illustration of the types of problems that are encountered, as listed in Table 6.1. We should note that in a materiallynonlinearonly analysis, the nonlinear effect lies only in the nonlinear stressstrain relation. The displacements and strains are infinitesimally small; therefore the usual engineering stress and strain measures can be employed in the response description. Considering the large displacement but small strain conditions, we note that in essence the material is subjected to infinitesimally small strains measured in a bodyattached coordinate frame x’, y’ while this frame undergoes large rigid body displacements and rotations. The stressstrain relation of the material can be linear o nonlinear, As shown in Fig. 6.1 and Table 6.1, the most general analysis case is the one in which the material is subjected to large displacements and large strains. In this case the stressstrain relation is also usually nonlinear. In addition to the analysis categories listed in Table 6.1, Fig. 6.1 illustrates another type of nonlinear analysis, namely, the analysis of problems in which the boundary conditions change during the motion of the body under consideration. This situation arises in
Sec. 6.1
487
Introduction to Nonlinear Analysis
£ 1073 cm
Au D= @i=1
®
=140
Ly
u
el =
.
+ Ay® = 6.6667 X 1072 cm
6.6667 X 1074
7
u® ∊≨↗∣⋟∶ L—=
—a
.. is elastic. < € — section a
5 .'. . 1.3333>