FEM 2024 Notes - Finite Element Methods - Aalborg University PDF

Notes for course on Finite Element Methods, 2024 Introduction Erik Lund & Esben Lindgaard Department of Materials and Production, Aalborg University, {el, elo}@mp.aau.dk. Motivation for the course The Finite Element Method (FEM) is the most universal numerical method for solution of boundary value problems expressed by partial differential equations, and it is therefore by far the most frequently used tool for stress, strain and displacement analysis of real-life structures, machines, and parts thereof. Thus, the central concepts, theories and methods taught in this course do not only constitute an impor- tant basis for the subsequent part of the curriculum, but also for the students’ professional work after graduation. Department of Materials and Production Introduction 1 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 1 Objective of the course 1/2 The objective of the course is stated in the Study Curriculum and repeated here for convenience. Goal: Students who complete the module are expected to: Knowledge: Be able to use the finite element method in static stress analysis. Have knowledge of element technology, such as bar, beam, solid and shell elements. Be able to solve structural dynamics and vibrations problems using methods such as free vibrations, modal methods and direct time integration methods. Be able to apply nonlinear finite element methods including solution of systems of nonlinear equations, geometrically nonlinear problems, contact problems, and nonlinear material models. Be able to perform linearised buckling analysis. Department of Materials and Production Introduction 2 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 2 Objective of the course 2/2 and Evaluation Skills: Demonstrate a basic understanding of concepts, theory and applications of finite element analysis from a mechanical engineering view point. Be able to perform linear and nonlinear static and dynamic stress analysis using commercial finite element software. Be able to conduct a systematic and critical assessment of analysis choices and obtained finite ele- ment results. Competences: Be able to apply the concepts, theories and techniques covered in the area of linear and nonlinear finite element analysis on engineering problems using commercial software programs. Be able to judge the opportunities and limitations of finite element simulations with regard to engi- neering problems. Form of examination: Internal, oral examination. Department of Materials and Production Introduction 3 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 3 Prerequisites Calculus, mechanics, solid mechanics and the courses Finite Element Methods 1 and Finite Element Methods 2 (or similar). However, several slides covering part of the two above mentioned courses have been added to these notes (pp. 19 - 82) in order to support the students that have not attended the above mentioned courses. The first lecture gives a quick summary of this material. Literature Cook et al. (2002): R.D. Cook, D.S. Malkus, M.E. Plesha & R.J. Witt: “Concepts and Applications of Finite Element Analysis”, Fourth Edition, John Wiley & Sons, New York, 2002. Lund & Lindgaard (2024): Erik Lund & Esben Lindgaard: “Notes for course on Finite Element Methods, 2024”. Furthermore, some figures are taken from Cook et al. (1989): Cook, R.D., Malkus, D.S. & Plesha, M.E.: “Concepts and Applications of Finite Element Analysis”, Third Edition, John Wiley & Sons, New York, 1989. For each lecture it is assumed, that the participants have read the course material, such that focus can be put on the difficult parts of the material. The reading guidelines are given on the following page. Department of Materials and Production Introduction 4 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 4 Outline of course and reading guidelines 1/7 Lecture 1: 4 hours lecture Lecturer: Esben Lindgaard Literature: Cook et al. (2002): Chapters 1, 2, 3, and 4. Lund & Lindgaard (2024): Notes for lecture 1. Contents: Introduction to FEM. The characteristic element matrix. Strong and weak form. Admissible configuration. Interpolation. 1D and 2D shape functions. Different approaches to FEM, e.g. variational methods and Galerkin FEM. Element stiffness matrices (1D bar element, 2D beam element and frame element, 2D bilinear solid element, and 3D trilinear solid element). Consistent load vectors. Solution of systems of linear equations. Equilibrium and compatibility of a FE solution. Convergence requirements. Properties of 2D bilinear element. Home work assignments: Approx. 4 hours. Important note 1: It is assumed that the course participants have passed finite element courses similar to the bachelor courses in this area given at AAU. Important note 2: The first lecture is a condensed summary of the material covered in the bachelor courses on finite element methods given at AAU and will serve as background for the remaining of the course. Only parts of the material will be presented while the rest is for self-study. This lecture is optional for those who have attended the bachelor courses, however we recommend you to participate anyway. Important note 3: The participants should install the ANSYS programme on their pc’s after this lecture and complete the “ANSYS Workbench Tutorial: Introduction to Static Structural” available at the course home page. Department of Materials and Production Introduction 5 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 5 Outline of course and reading guidelines 2/7 Lecture 2: 4 hours lecture Lecturer: Erik Lund Literature: Cook et al. (2002): Sections 6.1 - 6.13 and 7.1 - 7.3. Lund & Lindgaard (2024): Notes for lecture 2. Contents: Element technology. 2D isoparametric solid elements. Numerical integration. Static conden- sation and sub-structuring. Elements with interior degrees of freedom, for example the improved 4-node isoparametric 2D element (QM6). 3D solid elements. Shell elements. Home work assignments: Approx. 2 hours. Lecture 3: 2 hour lecture and 2 hours exercises Lecturer: Erik Lund Literature: Lund & Lindgaard (2024): Notes for lecture 3. Contents: Element technology. Symmetry conditions. Coordinate transformation. Basic use of ANSYS. Exercises to be solved using ANSYS Workbench. Home work assignments: Approx. 4 hour. Important note: The participants should have installed the ANSYS programme on their pc’s before the lecture and they should have completed the “ANSYS Workbench Tutorial: Introduction to Static Structural” available at the course home page. Department of Materials and Production Introduction 6 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 6 Outline of course and reading guidelines 3/7 Lecture 4: 1/2 hour student presentation of exercises from lecture 3, 2 hours lecture and 1 1/2 hours exercises Lecturer: Erik Lund Literature: Cook et al. (2002): Sections 2.10 - 2.11 and 12.1 - 12.2. Lund & Lindgaard (2024): Notes for lecture 4. Contents: Solution of multi-physics problems, exemplified by thermoelastic problems. Macros and pa- rameters in ANSYS. Combining ANSYS with other simulation environments. ANSYS Workbench. Home work assignments: Approx. 1 hour. Lecture 5: 1/2 hour student presentation of exercises from lecture 4, 2 hours lecture and 1 1/2 hours exercises Lecturer: Erik Lund Literature: Cook et al. (2002): Sections 9.1 and 9.8 - 9.11. Lund & Lindgaard (2024): Notes for lecture 5. Contents: Methods for error estimation and adaptive mesh generation. Home work assignments: Approx. 2 hours. Department of Materials and Production Introduction 7 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 7 Outline of course and reading guidelines 4/7 Lecture 6: 1/2 hour student presentation of exercises from lecture 5, 2 1/2 hours lecture and 1 hours exercises Lecturer: Erik Lund Literature: Cook et al. (2002): Sections 11.1 - 11.5 and 11.7. Lund & Lindgaard (2024): Notes for lecture 6. Contents: Structural dynamics and vibrations. Classification of dynamic problems. Dynamic equilibrium in finite element form derived using principle of virtual work. Problems of free vibration (eigenvalue problems). Mass and damping matrices. Free, undamped vibrations. Solution of eigenvalue problems. Methods of modal analysis. Home work assignments: Approx. 4 hours. Lecture 7: 1/2 hour student presentation of exercises from lecture 6, 2 hours lecture and 1 1/2 hours exercises Lecturer: Erik Lund Literature: Cook et al. (2002): Sections 11.11 - 11.13. Lund & Lindgaard (2024): Notes for lecture 7. Contents: Structural dynamics. Explicit and implicit direct methods for time integration. Home work assignments: Approx. 4 hours. Department of Materials and Production Introduction 8 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 8 Outline of course and reading guidelines 5/7 Lecture 8: 1/2 hour student presentation of exercises from lecture 7 and 3 1/2 hours lecture Lecturer: Esben Lindgaard Literature: Cook et al. (2002): Sections 17.1, 17.2, 17.9, 18.1 and 18.4 - 18.7. Lund & Lindgaard (2024): Notes for lecture 8. Contents: Introduction to nonlinear finite element methods. Classification of nonlinearities. Solution of systems of nonlinear equations using Newton-Raphson schemes or arc-length methods. General finite element formulation of geometrically nonlinear problems. The relation between geometrically nonlinear analysis and linearized buckling analysis. Home work assignments: Approx. 2 hours. Lecture 9: 1/2 hour lecture and 3 1/2 hours exercises Lecturer: Esben Lindgaard Literature: Lund & Lindgaard (2024): Notes for lecture 8. Contents: Exercises to be solved using ANSYS Workbench. Home work assignments: Approx. 4 hours. Department of Materials and Production Introduction 9 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 9 Outline of course and reading guidelines 6/7 Lecture 10: 1/2 hour student presentation of exercises from lecture 9, 2 hours lecture and 1 1/2 hours exercises Lecturer: Esben Lindgaard Literature: Cook et al. (2002): Section 17.8. Lund & Lindgaard (2024): Notes for lecture 10. Contents: Contact problems and their implementation in implicit and explicit finite element programs. Home work assignments: Approx. 4 hours. Lecture 11: 1/2 hour student presentation of exercises from lecture 10, 2 hours lecture and 1 1/2 hours exercises Lecturer: Esben Lindgaard Literature: Cook et al. (2002): Sections 17.3 - 17.6. Lund & Lindgaard (2024): Notes for lecture 11. Contents: Introduction to finite element modelling of nonlinear material models. Home work assignments: Approx. 4 hours. Department of Materials and Production Introduction 10 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 10 Outline of course and reading guidelines 7/7 Lecture 12: 1/2 hour student presentation of exercises from lecture 11 and 1 hour lecture Lecturer: Esben Lindgaard Literature: Lund & Lindgaard (2024): Notes for lecture 12. Contents: Introduction to workshop and presentation of oral examination questions. Home work assignments: Work on workshop projects - approx. 20 hours per student. Lecture 13: Workshop Lecturers: Esben Lindgaard and Erik Lund Contents: Workshop - student presentations. Feedback on the chosen problems and underlying theory. Important note 1: The workshop projects deal with finite element analysis of a load carrying mechanical component or system defined by the student. Important note 2: Note that the workshop can be considered as an informal seminar and therefore not directly a part of the evaluation of the course. However, it is a prerequisite for passing the course that the workshop project is completed, handed-in, and presented at the workshop! Important note 3: The oral examination is individually and will be based partly on the handed-in workshop project and one of the oral examination questions provided in lecture 12. Department of Materials and Production Introduction 11 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 11 Introduction to Workshop Projects (1/2) Motivation for the workshop projects is to apply the concepts, theories and techniques covered in the course of finite element methods on a practical problem of interest. Thus establishing a link between theory and practical problem solving skills. The project work should be accomplished by using the com- mercial finite element program ANSYS Workbench (or ANSYS Mechanical APDL). Approach for workshop projects Done in semester groups Selection of structural problem(s) of interest (component or system) Selection and definition of analysis studies Perform analyses on the structural problem(s) Document project work, i.e. problem setting, procedure, modelling issues, results etc. (slide show format - max 15-20 slides - please include names and slide numbers in slide show) Present results at the workshop (max. 15-20 minutes pr. group) The amount of work load for the workshop project is approximately 20 hours per student. Department of Materials and Production Lecture 1: Introduction to Workshop Projects 12 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 12 Introduction to Workshop Projects (2/2) Note that the workshop can be considered as an informal seminar and therefore not directly a part of the evaluation of the course. However, it is a prerequisite for passing the course that the workshop project is completed, handed-in, and presented at the workshop! All participants are expected to give feedback to the other presentations during the workshop. The oral examination is individual and will be based partly on the handed-in workshop project and one of the oral examination questions provided in lecture 12. Description of the oral examination The oral examination is individual and consists of two parts. Part one consists of a general overhearing and questioning about the handed-in course workshop project. Thus, remember to bring a printed copy of your handed-in workshop project. Part two includes draw of one exam question. A number of examination questions will be distributed in lecture 12 and thus known in advance. There is no preparation time. However, you may quickly look at your own notes for approx. 30 seconds in order to refresh your memory. After that no material, e.g. books, notes etc. may be used under the examination except for the handed-in workshop project. Under the oral examination, you are expected to demonstrate a basic understanding of concepts, theory, and applications of the finite element method and the chosen examination question. The examiners may interrupt your presentation with questions. The examination has 13 minutes duration. Department of Materials and Production Lecture 1: Description of the oral examination 13 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 13 Exercises using ANSYS Exercises for the course will be placed on the course webpage. We will use the commercial finite element (FE) code ANSYS for the exercises - please install ANSYS on your pc. For installation go to: \\m-tech.aau.dk\studies\software\FEM\ANSYS (open link using File Explorer). Recommended web pages: https://www.ansys.com/academic/students and https://courses.ansys.com/. At \\m-tech.aau.dk\studies\software\FEM\ANSYS\FEA_Best_Pract_2019R1_EN_Student (open link using File Explorer) you can find an ANSYS “FEA Best Practices” virtual classroom module. Reading guide for these notes/slides For each lecture there is a series of notes/slides as a help for the students during these lectures. Each lecture begins with an agenda and list of content to create an overview. These notes/slides do not give a complete and comprehensive description of the subject, and as such they should be viewed as supplementary information to the text book Cook et al. (2002). There is a notation list on page 16, which is consistent with Cook et al. (2002). There is an index list at the end of the slides. There is a reference list at the end of the slides. Department of Materials and Production Lecture 1: Description of the oral examination 14 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 14 A list of excellent FE text books In this course Cook et al. (2002) is used as text book, however, if you would like to know more about this subject it is recommended to explore the following books on the FEM: O.C. Zienkiewicz & R.L. Taylor: The Finite Element Method, Vol. 1, 2 & 3, Sixth Edition, Butterworth- Heinemann, Oxford, ISBN 0-7506-6431-2 (2005). T.J.R. Hughes: The Finite Element Method, Dover, New York, ISBN 0 4864 1181 8 (2000). K.-J. Bathe: Finite Element Procedures, Prentice-Hall (1996). C. A. Felippa: Introduction to Finite Element Methods (2013 ed.), Aerospace Engineering Sciences, Uni- versity of Colorado, Boulder, previously available at http://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/Home.html. C. A. Felippa: Advanced Finite Element Methods (2013 ed.), Aerospace Engineering Sciences, Uni- versity of Colorado, Boulder, previously available at http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/Home.html. Notice that the notes/book chapters written by Felippa are freely available and they are excellent! The versions of these notes listed above can be downloaded at the course home page. Department of Materials and Production Lecture 1: Description of the oral examination 15 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 15 List of nomenclature (only the most important symbols) 1/3 [] a matrix (rectangular or quadratic) ⌊⌋ a matrix row {} a vector ⌈⌋ a diagonal matrix φ a function for element interpolation {φe} nodal value of state variable φ [N ] matrix with interpolation functions Ni, i = 1,... , n, where n is the number of nodes in the element, i.e. φ = [N ]{φe} x, y, z global Cartesian coordinates xi, yi, zi global Cartesian coordinates for node i ξ, η, ζ, ξ1, ξ2, ξ3 local coordinates u, v, w translational degrees of freedom u,x compact notation for the partial deriverative ∂u ∂x ui , v i , w i translational degrees of freedom at node i θx , θy , θz rotation degrees of freedom Department of Materials and Production List of nomenclature 16 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 16 List of nomenclature (continued) 2/3 {d} element displacement vector {ǫ} strain vector [B] strain-displacement matrix [k] element stiffness matrix [E] constitutive matrix {re} consistent element load vector L, S, V length, surface and volume t thickness or time r radius (axial symmetry) [K] global stiffness matrix {D} global displacement vector {R} global load vector {P } global vector with concentrated nodal loads {F } nodal body forces per unit volume {Φ} nodal surface tractions [J] Jacobian matrix |J| determinant of the Jacobian matrix Department of Materials and Production List of nomenclature 17 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 17 List of nomenclature (continued) 3/3 {σ} stress vector U strain energy Ue element strain energy [m] element mass matrix [M ] global mass matrix [c] element damping matrix [C] global damping matrix ρ mass density κd viscous damping factor ˙ {d} element velocity vector or nodal velocities by element {Ḋ} global velocity vector or nodal velocities ¨ {d} element acceleration vector or nodal accelerations by element {D̈} global acceleration vector or nodal accelerations ω eigenvalue, eigenfrequency or circular frequency {D} nodal amplitudes, eigenvector or eigenmode {Z} modal coordinates or modal displacements [φ] modal matrix whose columns are the eigenvectors Department of Materials and Production List of nomenclature 18 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 18 Course on Finite Element Methods, 2024 Lecture 1 Lecturer: Esben Lindgaard Esben Lindgaard and Erik Lund, Department of Materials and Production, Aalborg University The following slides summarize much of the material covered in the previous courses on finite element methods. These slides serve as background material for the remaining of this course! Introduction. Basic idea of FEM. The characteristic element matrix. Strong form and weak form. Admissible configuration. Interpolation. 1D and 2D shape functions. Variational methods: Principle of stationary potential energy. Principle of virtual work. Galerkin FEM. Element stiffness matrices (1D bar element, 2D beam and frame elements, 2D bilinear solid element, and 3D trilinear solid element). Consistent load vectors. Solution methods for systems of linear equations. Equilibrium and compatibility of FE solution. Convergence requirement. Properties of 2D bilinear element. Department of Materials and Production Lecture 1: Introduction 19 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 19 Basic idea of FEM The FEM is the most versatile and applied numerical method for solving partial differential equations, which can describe physical problems such as: Heat conduction problems Structural mechanic problems Fluid problems (however, finite volume methods are generally preferred for these problems) Electro-static and magnetic-static problems etc. In this course we will mainly focus on static and dynamic stress/strain analyses, i.e. the governing equations are as follows: Equilibrium equations Geometric conditions Physical (constitutive) conditions and the compatibility equations. Additionally, the boundary conditions must be fulfilled. Essentially, the FEM introduces numerical models that replaces the governing partial differential equa- tions of a given problem with algebraic equations, which can be solved by standard equation solvers. Department of Materials and Production Lecture 1: Basic Idea of FEM 20 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 20 Very briefly, the idea behind the FEM is that we approximate the complicated governing equation of the problem with a model, which consists of simple and piecewise continuous solutions, i.e. the FEM consists of the following: Discretization. The structure is divided into a series of simple elements, which are defined by n nodes. One element is defined as a special domain, where the function φ is interpolated across this domain by the use of nodal values {φe} = {φ1, φ2,... , φn}. n is the number of nodes associated to the element. Interpolation. Simple piecewise continuous functions (Lagrange polynomials) are typically used for the interpolation in each element. However, in lecture 5 we will see the use of advanced hierarchical shape function (Legendre polynomials) for p-FEM. Cook et al. (2002) p. 5, fig. 1.2-1: Steps in modelling and example of FE discretization. Department of Materials and Production Lecture 1: Basic Idea of FEM 21 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 21 Cook et al. (2002) p. 5, fig. 1.3-1: FE analysis of a tapered bar Department of Materials and Production Lecture 1: Basic Idea of FEM 22 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 22 Cook et al. (1989) pp. 2, fig. 1.1-2 Cook et al. (1989) pp. 3, fig. 1.1-3 The smooth and continuous function φ is interpolated across a three-node triangular element by a linear polynomial φ̃ = a1 +a2 x+a3 y. The smooth and continuous function φ is interpolated across a four-node quadrilateral element by a bilinear polynomial φ̃ = a1 + a2x + a3y + a4xy. Department of Materials and Production Lecture 1: Basic Idea of FEM 23 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 23 The characteristic element matrix The characteristic element matrix depends on the physical problem that we would like to solve. In the field of mechanics: the element stiffness matrix relates global nodal displacement to nodal forces. Consequently, the state variables of the problem are global displacements. In the field of heat conduction: the element conductivity matrix relates nodal temperatures to nodal heat fluxes. The characteristic element matrix can be derived in three ways: 1. The direct method requires and is based on physical insight into the problem. This can only be used in a few simple cases. 2. Variational methods as the Rayleigh-Ritz-method, which is based on the principle of stationary po- tential energy or in other words the minimum of the total elastic potential energy. If the governing differential equation of the problem only contains differentials of equal-order, then there exists a func- tional to the problem. A functional is an integral expression describing the problem, which implicitly contains the governing partial differential equations and non-essential boundary conditions. An ex- ample would be the potential energy known from structural mechanics. 3. Weighted residual methods can be applied to all types of partial differential equations. An example is the Galerkin method, which results in the same discretized governing equations as achieved by the variational methods. Department of Materials and Production Lecture 1: The characteristic element matrix 24 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 24 Strong form and weak form It is essential to recognize that we want to solve a physical problem described by partial differential equa- tions written in ’strong form’. However, the FEM is based on integral expressions, i.e. the functional of the total potential energy depends not on the state variables and its derivatives at a particular point, but upon their integrated effect over the domain, i.e. the problem is written and solved in the ’weak form’. Consequently, the conditions that must be fulfilled such as equilibrium conditions and stress boundary conditions are only satisfied in an average sense, not at very point! Department of Materials and Production Lecture 1: Strong form and weak form 25 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 25 Strong form and weak form - explained For a static linear elastic structural problem we can describe it as follows Strong Form: Weak Form: Problem described by a set of partial differential Problem described by a functional giving the inte- equations that must be fulfilled at every point in gral effect over the structure. the structure. Z Z Z Π= f (x, y, z, u, v, w, ux, uy ,...)dxdydz σij,j + Fi = 0 1 P1 εij = (ui,j + uj,i) 2 E ν σij = (εij + δij εkk ) R 1+ν 1 − 2ν Π= σij (εij ) dV P1 z Pn x y σij (εij ) z Pn x y Department of Materials and Production Lecture 1: Strong form and weak form 26 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 26 Admissible configuration Every applied interpolation function in a FE model has to result in an admissible (deformed) configuration. The term admissible configuration is associated with the system that we want to describe. We define the system as being a physical structure with applied loads. The configuration of the system is the set of all positions of all particles in the structure. If the order 2m is the highest order of derivative of the state variable in the governing PDE, then the boundary conditions can be divided into the following classes: Essential boundary conditions, which also are called kinematic/geometric/displacement/primary Bound- ary Conditions (BC), involve derivatives of order 0 through m − 1, the zeroth derivative being the dependent variable itself. Nonessential boundary conditions, which also are called natural BC, involve derivatives of order m through 2m − 1. Problem bar beam bending 2D heat conduction PDE AEu,xx +q = 0 EIw,xxxx −q = 0 k∇2T + Q = 0 2m, m-1, 2m-1 2,0,1 4,1,3 2,0,1 essential BC on u only on w and w,x on T only nonessential BC on σx = Eu,x on M = EIw,xx on q = −k(T,x lB + T,y mB ) and V = EIw,xxx (i.e. q = −kT,n) Department of Materials and Production Lecture 1: Admissible configuration 27 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 27 An admissible configuration is any configuration that fulfills both internal compatibility and the essential boundary conditions, see Cook et al. (2002) pp. 83-87 and pp.137-138. Cook et al. (2002) p. 138, fig. 4.2-2. The upper configuration of the bending beam is not admissible, since it has four flaws, which ones? However, the two lower configurations are admissible, since the nonessential BC need not to be fulfilled. Based on the definition of the admissible configuration, we can establish two requirements for the inter- polation: Bars, 2D/3D solid models, etc., where 2m = 2, can be interpolated by piecewise continuous polyno- mials. This is called C0 continuity. Beams, plates and shells, where 2m = 4, need to be interpolated such that tangent continuity is satisfied. This is called C1 continuity. Department of Materials and Production Lecture 1: Admissible configuration 28 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 28 A field or polynomial function is said to have Cm continuity, if the field or polynomial function has continu- ous state variables up to m-order derivatives. dφ C0 continuity: φ is continuous, however φ,x = dx is not. C1 continuity: φ and φ,x are continuous. This is required for Bernoulli-Euler beam elements and Kirchhoff plate and shell elements. Cook et al. (2002) p. 85, fig. 3.2-1, the function φ1 is C0 continuous, while φ2 is C1 continuous. Department of Materials and Production Lecture 1: Admissible configuration 29 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 29 Interpolation Interpolation is demonstrated with the help of piecewise continuous polynomials for a bar element, see Cook et al. (2002) pp. 85-87 (Cook et al. (1989) fig. 3.8-1) The axial displacement u across the length from x = 0 to x = LT can be ap- proximated by three linear displacement fields (the e indicates approximated value) ũ = a1 + a2x for 0 ≤ x ≤ x2 ũ = a3 + a4x for x2 ≤ x ≤ x3 ũ = a5 + a6x for x3 ≤ x ≤ x4 Based on the requirement of an admissible configuration, the following displacement field is obtained: ũ = a2x for 0 ≤ x ≤ x2 ũ = a2x2 + a4(x − x2) for x2 ≤ x ≤ x3 ũ = a2x2 + a4(x3 − x2) + a6(x − x3) for x3 ≤ x ≤ x4 The adopted approach for determining the displacement field is not practical, because the constants ai give no physical meaning and since the approach is problem dependent. Instead we will describe the displacements in terms of nodal displacements (degree of freedoms) ui. Department of Materials and Production Lecture 1: Interpolation 30 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 30 The axial displacement field across the generic element shown below is assumed to be linear in the axial coordinate s. The assumed displacement field ũ must fulfill that ũ = ui at the left node and ũ = uj at the right node. Consequently, we can now interpolate ũ linearly along the element by: n X L−s s ũ = ui + uj or ũ = ⌊N ⌋{d} = Ni d i L L i=1 where we have defined the shape function matrix [N ] and the element nodal displacement vector {d}: L−s s ui ⌊N ⌋ = ⌊ ⌋ and {d} = L L uj The elements Ni in the shape function matrix are called shape functions, basis functions or interpolation functions. The degrees of freedom u1 and u2 are coefficients for the shape functions Ni in the approxi- mated expression for the displacements, i.e. u1 and u2 play the same role as constants a1 and a2 in the interpolation function ũ = a1 + a2x. Department of Materials and Production Lecture 1: Interpolation 31 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 31 Determination of the 1D shape function by the a-basis-method The shape functions have intuitively been derived on the previous slides, however, this can be done in a systematic way by the a-basis-method. Based on the generic bar element on the previous slide, the linear displacement field can be written as follows: a1 ũ = a1 + a2s or ũ = ⌊1 s⌋ a2 As previously, the requirement is that ũ = ui for s = 0 and ũ = uj for s = L: ui 1 0 a1 = or {d} = [A]{a} uj 1 L a2 This equation for {a} is then solved and substituted into the upper equation: ui {a} = [A]−1{d} and ũ = ⌊1 s⌋[A]−1 uj n X By using the definition of the interpolation, ũ = ⌊N ⌋{d} = Nidi, it can be written as: −1 i=1 " # 1 0 1 0 L−s s −1 ⌊N ⌋ = ⌊1 s⌋[A] = ⌊1 s⌋ = ⌊1 s⌋ 1 1 =⌊ ⌋ 1 L − L L L L Department of Materials and Production Lecture 1: 1D and 2D shape functions 32 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 32 The expression found by the a-basis-method is equivalent to the intuitively derived expression for the shape function matrix ⌊N ⌋. The same approach can be used if the interpolation functions are chosen to be, e.g., quadratic polynomials. Alternatively, one can use Lagrange’s Pn interpolation formula, see Cook et al. (2002) pp. 86-87 and fig. 3.2-4, where φ̃ = ⌊N ⌋{φe} = i=1 Ni φi The characteristics for all C0 polynomial shape functions: 1. All shape functions Ni, i = 1,... , n and the function φ are polynomials of the same order. 2. Ni = 1 when x = xi and Ni = 0 when x = xj. n X 3. Ni = 1. i=1 Department of Materials and Production Lecture 1: 1D and 2D shape functions 33 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 33 Determination of the 2D shape functions for a 4 node bilinear element We want to determine the shape functions for a 2D solid 4 node bilinear element. Here φ is assumed on the form φ̃ = a1 + a2x + a3y + a4xy. If we introduce values at the edge of the element between node 1 and 4 as φ̃14 and similarly with φ̃23 (y -direction), then: b−y b+y b−y b+y φ̃14 = φ1 + φ4 and φ̃23 = φ2 + φ3 2b 2b 2b 2b Cook et al. (2002) p. 97, fig. 3.6-1. Next we interpolate in the x-direction using φ̃14 and φ̃23: a−x a+x a−x b−y b+y a+x b−y b+y φ̃ = φ̃14 + φ̃23 = φ1 + φ4 + φ2 + φ3 2a 2a 2a 2b 2b 2a 2b 2b = N1 φ 1 + N2 φ 2 + N3 φ 3 + N4 φ 4 where we have obtained (a − x)(b − y) (a + x)(b − y) (a + x)(b + y) (a − x)(b + y) N1 = , N2 = , N3 = , N4 = 4ab 4ab 4ab 4ab Department of Materials and Production Lecture 1: 1D and 2D shape functions 34 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 34 Short introduction to variational methods As previously mentioned, the characteristic element matrix can be derived by three different approaches: 1. Direct methods 2. Variational methods, e.g. the principle of stationary potential energy or the principle of virtual work 3. Weighted residual methods, e.g. Galerkin FEM. Variational methods, such as the principle of stationary potential energy and the principle of virtual work, are based on integral expressions describing the problem, which implicitly contains the governing equa- tions and non-essential boundary conditions. If no approximation is introduced, the variational methods will produce exact solutions. However, varia- tional methods are often combined with the Rayleigh-Ritz methods, thus generally producing approximate solutions where equilibrium conditions and non-essential boundary conditions only are satisfied in an av- erage or integral sense, not at every point! The theoretical background for the calculus of variations and energy methods can be found in e.g. C.L. Dym & I.H. Shames. Solid Mechanics. A Variational Approach, Augmented Edition. Springer New York 2013, Chapter 2 and 3. Pdf version available through www.aub.aau.dk. Let us start out by taking a closer look at the principle of stationary potential energy. Department of Materials and Production Lecture 1: Short introduction to variational methods 35 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 35 Short introduction to the principle of stationary potential energy If the governing differential equation of the problem only contains derivatives of equal order, then there exists a functional to the problem, i.e. we can apply the calculus of variations. A functional is an integral expression describing the problem, which implicitly contains the governing partial differential equations and non-essential boundary conditions. An example would be the potential energy known from structural mechanics. The definition of a potential within the field of structural mechanics: The ability to perform work. Requirement: conservative, elastic system, i.e. the work is path-independent and there is energy con- servation. The total elastic potential Πp consists of two parts: The internal potential Πip , which is the stored elastic energy, i.e. the strain energy U The external potential Πep = Ω, i.e. the potential of work of the external forces i.e. Πp = Πip + Πep = U + Ω Department of Materials and Production Lecture 1: Principle of Stationary Potential Energy 36 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 36 The principle of stationary potential energy (the minimum of the total potential energy) is given as (see Cook et al. (2002) pp. 137-138): Among all the admissible configurations of a conservative system, those that satisfy the equations of equilibrium make the potential energy stationary w.r.t. small admissible variations of displacements. The potential Πp of a discrete system described by the global displacement vector {D} with the elements {D1, D2,... , Dn}T is a function of Di, i.e. Πp = Πp (D1, D2,... , Dn). Consequently, the principle of stationary potential energy can be written as: ∂Πp ∂Πp ∂Πp dΠp = dD1 + dD2 +... + dDn = 0 ∂D1 ∂D2 ∂Dn This expression has to be valid for any small admissible variation of the configuration {D}, such that the following is required: ∂Πp ∂Πp = 0 or = {0} ∂Di ∂{D} In the following, the principle of stationary potential energy will be demonstrated for a 1D spring element in equilibrium. Department of Materials and Production Lecture 1: Principle of Stationary Potential Energy 37 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 37 1D example: The principle of stationary potential energy of a spring The spring has a stiffness coefficient k and a length of L in the undeformed configuration. The force P is applied and the spring is elongated by D. We will use an imaginary wall a distance H away from the right end of the spring. The distance H is arbitrarily chosen, however it is chosen such that the spring does not deform the distance D beyond the wall distance H. The potential energy Πp consists of to contributions: Πp = Πip + Πep = U + Ω where the stored elastic energy is 1 Πip = U = k D2 2 and the external potential of the external loads is given by: Πep = Ω = P (H − D) This is reasoned by the argument that the external force P has done the work P D by elongating the spring with D , i.e. the force has lost the potential P D of the total available potential P H. Modified version of Cook et al. (2002) p. 139, fig. 4.2-3 Department of Materials and Production Lecture 1: Principle of Stationary Potential Energy 38 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 38 We are now able to set up the potential energy Πp: 1 Πp = Πip + Πep = U + Ω = k D2 + P (H − D) 2 The equilibrium configuration Deq is found by the principle of stationary potential energy: 1 2 ∂Πp dΠp d 2 k D + P (H − D) = = = k D−P = 0 ∂D dD dD i.e. we have obtained the governing equation of the system: kD = P Cook et al. (2002) p. 140, fig. 4.2-4: Notice that Cook et al. (2002) chooses H to be equal to 0, which can be confusing. Department of Materials and Production Lecture 1: Principle of Stationary Potential Energy 39 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 39 The potential energy of an elastic body The potential energy of an elastic body is described in Cook et al. (2002) pp. 142-146. Here, a summation of the expression is given, so please read the pages in Cook et al. (2002) or in another book such as C.L. Dym & I.H. Shames. Solid Mechanics. A Variational Approach, Augmented Edition. Springer New York 2013. Here the derivation of the strain energy U for the elastic body can be found: Z Z Z ǫij Z 1 U = u0 dV = σij dǫij dV = σij ǫij dV V V 0 V 2 where u0 is strain energy density. The potential energy is the sum of the stored strain energy U and the potential of the external forces Ω, when neglecting initial strains and stresses: Z Z Z 1 T Πp = {ǫ} [E]{ǫ}dV − {u}T {F } dV − {u}T {Φ} dS − {D}T {P } V 2 V S Here, {ǫ} is the strain vector, {σ} the stress vector, [E] the constitutive matrix, {F } the body forces, {φ} the surface tractions, V the volume and S the surface. Additionally, Hooke’s law of linear elasticity is given by {σ} = [E]{ǫ}, which has been used in the expression for the strain energy U. Let us next consider how to apply the Rayleigh-Ritz method in order to determine approximate solutions for our structural problems. Department of Materials and Production Lecture 1: Principle of Stationary Potential Energy 40 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 40 The Rayleigh-Ritz method Basic idea behind Rayleigh-Ritz is to avoid solving very complicated differential equations by applying the method to a functional that describes the mathematical model. The result is a substitute model with finite degrees of freedom which is described by algebraic equations rather than by differential equations. From the principle of stationary potential energy we know that: Among all the admissible configurations of a conservative system, those that satisfy the equations of equilibrium make the potential energy stationary w.r.t. small admissible variations of displacements. In the Rayleigh-Ritz method we determine a solution based upon a subset of all admissible configurations. This subset we refer to as the approximating field. Polynomials are usually applied as approximating fields since they automatically satisfy compatibility conditions, e.g. ũ = a1 + a2x + a3x2 +... + anxn−1. The e indicates that we have introduced an approximation via the approximating field. e p, the equilibrium configuration is By the stationary value of the potential energy of the substitute model Π defined by the n algebraic equations ep ∂Π = 0 for i = 1,... , n ∂ai e p is minimum. which solution gives the unknown parameters ai such that the model potential Π Department of Materials and Production Lecture 1: The Rayleigh-Ritz method 41 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 41 The Rayleigh-Ritz method ep Πp , Π The Rayleigh-Ritz method yields the exact solution if the exact field is included in the approximating field {ũ}. In such case Π e p,min = Πp,min. The Rayleigh-Ritz method: {u} ⇒ {ũ} The Rayleigh-Ritz method has a classical form and ⇓ a modern finite element form. In the classical form a ep Πp ⇒ Π single approximating field is used to model the entire e p,min Π e p , the proposed displacement By a minimization of Π structure whereas in the FE form the structure is di- field {ũ} is approximated as close as possible to the vided into many subdomains (finite elements) each Πp,min exact equilibrium configuration. having a simple approximating field. Let us now consider the modern version of the Rayleigh-Ritz method where the approximating field is the interpolation functions for all finite elements in the structure. We will insert this approximating displacement field {ũ} in the expressions of the potential energy, make it stationary and from that derive the governing equilibrium equations on finite element form. Consequently, the elastic body is now described by a finite number of parameters (element properties). Department of Materials and Production Lecture 1: The Rayleigh-Ritz method 42 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 42 Element stiffness matrix and consistent load vector based on the principle of station- ary potential energy Finite element discretization of a 2D structure: The basis of displacement based FEA: displacements at any point within the elastic body/continuum can be interpolated on the basis of the nodal displacements. We will denote the approximate displacements as {ũ}. In case of a 2D solid problem the displacements are given as: where [N ] contains the shape functions (interpolation functions) ũ {d} = ⌊u1 v1 u2 v2... un vn⌋T {ũ} = = [N ]{d} ṽ = ⌊{d1}T {d2}T... {dn}T ⌋T n is the number of nodes for the element considered Department of Materials and Production Lecture 1: The Rayleigh-Ritz method on FE form 43 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 43 Consequently, the generalized expression for the element strain vector is then written as: {ǫ̃} = [∂]{ũ} = [∂][N ]{d} ≡ [B]{d} where [B] ≡ [∂][N ] is the strain-displacement matrix This is inserted into the potential energy, where ne is the number of elements in the elastic body: X X ne ne ep = Π ep + Π i ep = U e e+Ω e= 1 T Π {d} [k]{d} − {d}T {re} − {D}T {P } 2 j=1 j=1 where the symmetric element stiffness matrix [k] and the consistent load vector {re } is introduced as: Z Z Z [k] = [B]T [E][B]dV {re} = [N ]T {F }dV + [N ]T {Φ}dS Ve Ve Se This is expanded to global level for the entire elastic body: ne X ne X [K] = [k], {R} = {re} + {P } j=1 j=1 The principle of stationary potential energy gives the governing equations: 1 ep ∂Π ∂ 2 T {D} [K] {D} − {D} {R} T = = [K] {D} − {R} = {0} , i.e., [K]{D} = {R} ∂{D} ∂{D} Department of Materials and Production Lecture 1: The Rayleigh-Ritz method on FE form 44 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 44 The principle of virtual work Another simple method to derive the governing equations is the principle of virtual work also called the principle of virtual displacements. This method is even more general than the principle of stationary potential energy, since it is valid for any constitutive law and thereby not, as the principle of stationary potential energy, restricted to elastic bodies. In fact stationary potential energy can be considered as a special case of the principle of virtual work. The theoretical background for this principle can be found in e.g. C.L. Dym & I.H. Shames. Solid Mechanics. A Variational Approach, Augmented Edition. Springer New York 2013, pp. 7-9, pp. 117-120, pp. 373-376, pp. 559-566, pdf version available through www.aub.aau.dk. The principle of virtual work: It is required that the work done by external forces, body forces, and internal stresses is in equilibrium when the deformable body is exposed to a virtual, kinematically admissible deformed configuration (i.e., a virtual/imaginary displacement which fulfills compatibility and essential boundary conditions). Department of Materials and Production Lecture 1: Principle of Virtual Work (self-study) 45 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 45 The symbol δ has almost the same Index list for the finite element formulation: meaning as d for differential, but {F } : body forces {Φ} : surface tractions the symbol δ is used when consid- {p}i : concentrated nodal loads {u} : displacements ering a virtual displacement. {δu} : virtual displacement {δǫ} : virtual strain due to {δu} The principle of virtual work for a single element with the volume Ve can be expressed as follows: Z Z Z n X {δǫ}T {σ}dV = {δu}T {F }dV + {δu}T {Φ}dS + {δu}Ti {p}i | Ve {z } | Ve Se {z i=1 } work of the internal stresses = work done by body forces, surface forces, change in elastic energy and nodal forces This equation expresses that the incremental storage of strain energy is equal to the work done by the external forces on the deformable body in equilibrium, when it is exposed to a static and virtual kinematic admissible displacement {δu}. It should be noticed that this equation also can be derived by multiplying the equilibrium equation σij,j + Fi = 0 with δui, and perform integration over the entire volume, which is illustrated by a bar element in Cook et al. (2002) pp. 155-156. Now, we will introduce the interpolation {ũ} = [N ]{d}, where {ũ} is the approximated displacements at an arbitrary point within the element volume, [N ] is the shape function matrix and {d} is the nodal displacement vector for the element. Department of Materials and Production Lecture 1: Principle of Virtual Work (self-study) 46 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 46 The linearized geometric equations state that {ǫ̃} = [∂]{ũ}, i.e. {ǫ̃} = [∂][N ]{d} = [B]{d}, where [B] ≡ [∂][N ] The matrix [B] is called the strain-displacement matrix, since it gives the relationship between the strains {ǫ̃} at a point and the nodal element displacement vector {d}. Applying a virtual displacement yields: {δ ũ} = δ([N ]{d}) = [N ]{δd} and {δǫ̃} = δ([B]{d}) = [B]{δd} i.e. {δ ũ}T = {δd}T [N ]T and {δǫ̃}T = {δd}T [B]T Inserting this into the equilibrium equation for a single element yields: Z Z Z n ! X {δd}T [B]T {σ̃}dV = {δd}T [N ]T {F }dV + [N ]T {Φ}dS + {p}i | Ve {z } | Ve Se {z i=1 } { rint } external load vector {rext} Pn Here i=1 {p}i assembles the concentrated point loads, assumed to be located at nodes, into an ele- ment nodal load vector. The constitutive equation (Hooke’s law for a linear elastic material) states that {σ̃} = [E]{ǫ̃} (neglecting initial stresses or strains), which is inserted as well: Z Z Z n ! X {δd}T [B]T [E][B]dV {d} = {δd}T [N ]T {F }dV + [N ]T {Φ}dS + {p}i | Ve {z } | Ve Se {z i=1 } { rint } external load vector{rext} Department of Materials and Production Lecture 1: Principle of Virtual Work (self-study) 47 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 47 This expression is valid for any arbitrary, kinematically admissible nodal displacement {δd} ⇒ [k]{d} = {rext} where the element stiffness matrix [k] is given by Z [k] = [B]T [E][B]dV Ve and the external load vector is given by Z Z n X {rext} = [N ]T {F }dV + [N ]T {Φ}dS + {p}i Ve Se i=1 For a model with ne elements, the stiffness contributions are added together to form the global stiffness matrix, i.e. ne X ne X [K]{D} = {R} where [K] = [k] , {R} = {rext} n=1 n=1 The global nodal displacement vector {D} consists of the unknown nodal displacements in the entire model. Consequently, we have now obtained a linear algebraic system of equations, which describes the static equilibrium equations. Department of Materials and Production Lecture 1: Principle of Virtual Work (self-study) 48 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 48 Example of the weighted residual method: Galerkin FEM The basic idea with weighted residual methods is that if we introduce an interpolated expression for the state variable in the governing equations, then there will be introduced a residual R , i.e. an error due to the interpolated state variable ũ used as approximation for the real state variable u. This residual is then weighted in an integral sense thereby obtaining a FE formulation, i.e. a weak formulation of the governing equations. The Galerkin method is one among many weighted residual methods: “Collocation method”: The residual R is required to be zero in a series of points. “Subdomain method”: The integral of R is required to be zero over a series of subdomains. “Least squares method”: By using the least squares method the second power of R is minimized over a domain: ∂I = 0, i = 1,... , n ∂ai Z where I = R2dV is the squared residual over the volume V and ai contains the constants in the V interpolation expressions. Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 49 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 49 “The Galerkin method:” The integral of the residual is weighted by the shape functions and is required to be equal to zero: Z Ri = Wi R dV = 0, i = 1,... , n ∂ ũ where the weight functions Wi are given by: Wi = = Ni , i = 1,... , n ∂ai Galerkin FEM yields identical results with the variational methods (if there exists a functional for the prob- lem), when the same approximation ũ (the same “trial function”) is used. As previously mentioned, the benefit of a weighted residual method is that it can be used for problems for which there does not exist a functional, for example if the governing equation contains derivatives of odd order. An example is the Navier-Stokes equations. Remember that the approximated ũ must be an admissible configuration (see p. 27). Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 50 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 50 We will in the following derive the Galerkin FEM for a simple bar loaded by an axial tip force P and an axial linearly varying force q(x) = cx, see Cook et al. (2002) p. 180, fig. 5.1-1: A is the cross-sectional area and E is the elastic modulus. These properties are assumed to be constant along the length of the bar. The governing equations for this problem are then written as: Equilibrium equation: Aσx,x + q = 0 Geometric condition: ǫx = u,x Physical (constitutive) condition: σx = Eǫx together with the requirement of compatibility, which is fulfilled if the approximated state variable ũ is an admissible configuration. Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 51 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 51 The governing equation is obtained by using the equilibrium, geometric and constitutive conditions: AEu,xx +q = 0 An axial tip force P acts at x = LT , such that the nonessential boundary condition is (2m=2, 2m-1=1): AEu,x −P = 0 The essential boundary condition is only on u (2m=2, m-1=0). Thus we can find the exact analytical solution by integration. However, we will imagine that it is not possible and divide the bar into ne elements. Every element is associated with an approximated and linearly varying variable ũ, i.e.:   L−x x  ⌊N ⌋ = ⌊ ⌋ ũ = ⌊N ⌋{d} where L L  u  {d} = u1 2 The degrees of freedom u1 and u2 are the coefficients for the shape functions Ni in the approximated displacement ũ, i.e. u1 and u2 act just as the constants a1 and a2 in the expression: ũ = a1 + a2x Consequently, the weight functions in the Galerkin method are as follows: ∂ ũ L−x x Wi = = Ni where N1 = , N2 = ∂di L L Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 52 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 52 We are now ready to set up the residual equation (Galerkin) by inserting the approximation ũ in the governing equation. This will lead to an error, also called the residual R. The Galerkin finite element method is obtained by weighting the residual with the shape functions and require that the integral of these values is zero. Z ne X L Ni (AE ũ,xx +q) dx = 0 j=1 0 Cook et al. (2002) p. 188, fig. 5.3-2 The figure displays the shape functions on the structural level after the bar elements have been assem- bled. The number of shape functions and degrees of freedom is ne + 1, and the number of Galerkin residual equations is ne + 1. It is seen that in order to fulfill compatibility between the elements, we must have continuity in ũ and ũ,x. This requirement can be reduced only to require continuity in ũ by using integration by parts. Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 53 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 53 The previous statement can be seen from the fact that if the integral contains derivatives of order 2m, then the integral will only be defined if all derivatives up to order 2m − 1 of the integrand are continuous. In case of the bar element 2m − 1 = 2 − 1 = 1 but by integration by parts this can be reduced to 2m − 1 = 1 − 1 = 0! R R Using integration by parts u dv = u v − v du for each element yields the natural boundary conditions for the Galerkin method: Z L Z L NiAE ũ,xx dx = [NiAE ũ,x ]L0 − Ni,xAE ũ,x dx 0 | {z } 0 nonesssentiel BC At arbitrary x, axial force in the bar is: AE ũ,x −F = 0 ne Z X L and the Galerkin residual equation is: Ni (AE ũ,xx +q) dx = 0 j=1 0 Xne Z L ne X We obtain by inserting: (−Ni,xAE ũ,x +Niq) dx + [NiF ]L0 = 0 j=1 0 j=1 We will then introduce the strain-displacement matrix ⌊B⌋ −1 1 u1 ⌊B⌋ = ⌊N, x⌋ and ǫ̃x = ũ,x = ⌊B⌋{d} = L L u2 Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 54 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 54 By introducing and inserting ⌊B⌋: ne Z X L ne Z X L ne X ⌊B⌋T AE⌊B⌋dx {d} = ⌊N ⌋T q dx + [⌊N ⌋T P ]L0 j=1 |0 {z } j=1 |0 {z } j=1 [k]j {re }j ne X Since ⌊N ⌋0 = ⌊1 0⌋ and ⌊N ⌋L = ⌊0 1⌋ the last term [⌊N ⌋T F ]L0 gives a global vector consisting of the j=1 concentrated nodal loads. This vector is called {P }. Additionally, the element stiffness matrix [k]j and the consistent element load vector {re }j for the j th element are defined as Z L Z L [k]j = ⌊B⌋T AE⌊B⌋dx {re}j = ⌊N ⌋T q dx 0 0 By summation of all contributions from the element stiffness matrices the global stiffness matrix [K] is introduced: Xne Xne Z L [K] = [k]j = ⌊B⌋T AE⌊B⌋dx j=1 j=1 0 Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 55 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 55 By summation of all contributions from the consistent element load vectors the global load vector {R} is introduced: Z ne X ne X L {R} = {re}j + {P } = ⌊N ⌋T q dx + {P } j=1 j=1 0 where the concentrated nodal loads {P } were introduced as: ne X {P } = [⌊N ⌋T F ]L0 j=1 Consequently, the system of equations is as follows: [K]{D} = {R} which contains all the governing equations for the problem. The Galerkin finite element formulation allows us to set up the governing system of algebraic equations for an axially loaded bar. The approach is the same for more advanced PDE’s, see Cook et al. (2002) pp. 191-195 for a 2D plane elasticity theoretical formulation. The governing equations obtained by the principle of stationary potential energy are the same as obtained with the Galerkin FEM. See further information in Cook et al. (2002) chapter 5. Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 56 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 56 Comparisons between the residual methods Let us consider the simple bar element previously introduced by different weighted residual methods (“Collocation method”, “Subdomain method”, “Least squares method” and “Galerkin method”). The ex- ample is taken from Cook et al. (2002) section 5.1-5.2, pp. 180-186. The solution for this problem by the different methods is shown in Cook et al. (2002) p. 186, table 5.2- 1. Firstly, we see that none of the methods gives the exact solution through-out the bar on both the displacement and strain. However, the Galerkin method yields the exact values of the state variable u (this is only valid for this example, and this is not a generality!). None of the methods results in the correct strain through-out the bar element. Department of Materials and Production Lecture 1: Galerkin FEM (self-study) 57 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 57 Element stiffness matrices When deriving the element stiffness matrix one needs to start with the strain-displacement relations (geometric conditions) for the considered element. This will enable one to set up the strain-displacement matrix, when the interpolation (shape) functions have been chosen. Next the element stiffness matrix [k] and the consistent load vector {re} can be established on the basis of the general equations derived by either the principle of virtual work, the principle of stationary potential energy or Galerkin FEM.  L−x x  1D bar element  ⌊N ⌋ = ⌊ ⌋ ũ = ⌊N ⌋{d} where L L  u Cook et al. (1989) p. 113, fig. 4.2-1:  {d} = u1 2 The strain-displacement relations are as follows: d (⌊N ⌋{d}) d ⌊N ⌋ ǫ̃x = dũ/dx = = {d} = ⌊B⌋{d} dx dx since {d} are fixed values used for interpolations. d ⌊N ⌋ −1 1 The strain-displacement matrix is then given by: ⌊B⌋ = = ⌊ ⌋ dx Z L L L T EA 1 −1 This is inserted into the expression for the stiffness matrix: [k] = ⌊B⌋ AE⌊B⌋dx = 0 L −1 1 Department of Materials and Production Lecture 1: Element stiffness matrices 58 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 58 ”Plane” Bernoulli-Euler beam element The stiffness matrix of the Bernoulli-Euler beam element is derived by the standard beam theory, where the transverse shear deformation is ignored, see Cook et al. (1989) p. 16, fig. 1.5-2: The rotations w,x are assumed to be small and the cross-sectional areas (transverse normals) are as- sumed to remain straight and normal to the deformed neutral line after deformation. The axial displacement u is given by u = −z w,x. Since ǫx = u,x, the strain can be computed as ǫx = −z w,xx where w,xx is the curvature of the beam. Department of Materials and Production Lecture 1: Element stiffness matrices 59 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 59 In case of bending problems for beams, plates and shells, we would like to have a compact notation for the strain displacement relations as previously, therefore we define the [B ∗] matrix such that [B ∗]{d} represents the curvature. [B ∗] is used in these notes to emphasize its difference to [B]. Consequently, the Bernoulli-Euler beam element has the following interpolation functions: 2 d w̃ = ⌊N ⌋{d} and w̃,xx = ⌊B ∗⌋{d} where ⌊B⌋ = −z⌊B ∗⌋ = −z 2 ⌊N ⌋ dx In order to determine the matrix ⌊B ∗⌋ we need to derive an expression for the shape function [N ] for the beam element. An admissible configuration for a beam element requires continuity in both the interpolated value φ̃ and its derivative φ̃,x. This is called Hermitian interpolation. Cook et al. (2002) p. 87, fig. 3.2-3 We use a cubic polynomial for interpolation: φ̃ = a1 + a2x + a3x2 + a4x3 or φ̃ = ⌊X⌋{a} where ⌊X⌋ = ⌊1 x x2 x3⌋ and {a} = ⌊a1 a2 a3 a4⌋T We then require that: φ̃ = φ1 and φ̃,x = θ1 for x = 0 φ̃ = φ2 and φ̃,x = θ2 for x = L Department of Materials and Production Lecture 1: Element stiffness matrices 60 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 60 These requirements are inserted into the cubic expression for φ̃:      1 0 0 0  φ1     a1    θ1  0 1 0 0  a2 =  2 3  or {d} = [A]{a}  φ2    1 L L L 2   a3   θ2 0 1 2L 3L a4 Thus {a} = [A]−1{d}, and we obtain the following: φ̃ = ⌊X⌋{a} = ⌊X⌋[A]−1{d} = ⌊N ⌋{d} where the shape functions ⌊N ⌋ are given by ⌊N ⌋ = ⌊X⌋[A]−1 The four shape functions are dis- played on Cook et al. (2002) p. 87, fig. 3.2-4. Department of Materials and Production Lecture 1: Element stiffness matrices 61 / 82 Aalborg University, Denmark Lund & Lindgaard: Notes for course on Finite Element Methods, 2024 p. 61 Now we can readily derive the matrix ⌊B ∗⌋, which is used in the element stiffness matrix. The displacement w̃ is interpolated by a cubic expression: Beam element with 4 de- grees of freedom: Cook et al. w̃ = ⌊N ⌋{d} = ⌊N ⌋ {w1 θ1 w2 θ2}T (1989) p. 113, fig. 4.2-2: where 3x2 2x3 2x2 x3 3x2 2x3 x2 x3 ⌊N ⌋ = ⌊1 − 2 + 3 x− + 2 − 3 − + 2⌋ L L L L L2 L L L Thus ⌊B ∗⌋ is determined as ∗ d2 6 12x 4 6x 6 12x 2 6x ⌊B ⌋ = ⌊N ⌋ = ⌊− + − + 2 − − + 2⌋ dx2 L2 L3 L L L2 L3 L L Thus the element stiffness matrix is calculated as:   Z L 12 6L −12 6L EI  6L 4L2 −6L 2L2  [k] = ⌊B ∗⌋T EI⌊B ∗⌋dx = 3  0 L −12 −6L 12 −6L  6L 2L2 −6L 4L2 where the constitutive element matrix [E] is replaced with EI due to the z 2-term, which arises from the chosen definition of [B ∗], see also Cook et al. (2002) pp. 24-29 and pp. 90-91. Department of Materials and Production

FEM 2024 Notes - Finite Element Methods - Aalborg University PDF

Document Details

Tags

Related

Summary

Full Transcript

Upgrade to continue