Finite Element Method: Physics and Solution Methods

By Sinan Muftu

Finite Element Method: Physics and Solution Methods aims to provide the reader a sound understanding of the physical systems and solution methods to enable effective use of the finite element method. This book focuses on one- and two-dimensional elasticity and heat transfer problems with detailed derivations of the governing equations. The connections between the classical variational techniques and the finite element method are carefully explained. Following the chapter addressing the classical variational methods, the finite element method is developed as a natural outcome of these methods where the governing partial differential equation is defined over a subsegment (element) of the solution domain. As well as being a guide to thorough and effective use of the finite element method, this book also functions as a reference on theory of elasticity, heat transfer, and mechanics of beams.

• Covers the detailed physics governing the physical systems and the computational methods that provide engineering solutions in one place, encouraging the reader to conduct fully informed finite element analysis
• Addresses the methodology for modeling heat transfer, elasticity, and structural mechanics problems
• Extensive worked examples are provided to help the reader to understand how to apply these methods in practice

• Language: English
• Publisher: Academic Press
• Release date: Jul 14, 2022
• ISBN: 9780128232002

Sinan Muftu

Sinan Müftü is a Professor of Mechanical Engineering at Northeastern University, Boston, USA. His research is in the general area of applied mechanics with applications in tribology and bioengineering, including mechanics of axially translating materials for roll-2-roll manufacturing systems, mechanics of high velocity particle impacts for cold spray additive manufacturing, and structure-function relationships in biological systems. He has taught the finite element method to undergraduate and graduate students in his institution since 2004 and developed customized programs for numerical analysis throughout his career. This book comes out of his experience and observations in teaching and conducting research in applied numerical analysis. Dr. Müftü is an elected fellow of the American Society of Mechanical Engineers for his contributions to mechanics of axially translating media.

Book preview

Chapter 1

Introduction

One of the major aims of natural sciences is to identify the cause–effect relationships in nature. Grasping these relationships is not only useful to satisfy inherent curiosity, but it is also helpful in building and designing systems that improve our daily lives.

Engineers and scientists interested in developing predictive capabilities for natural processes are generally involved in two major tasks: modeling and simulation. Our understanding of the term model refers to developing a mathematical representation of a physical process that is observed in nature. To this end, partial differential equations (PDEs) are typically used to formulate mathematical representations of physical processes. The term simulation, on the other hand, refers to solving the mathematical model and analyzing the results to gain insight into the cause–effects relationships.

1.1 Modeling and simulation

1.1.1 Boundary and initial value problems

Consider a system with dependent variables u, v, and w defined over a domain Ω, which itself occupies a subsection of space (Fig. 1.1). In general, each variable can take different values at different points in the domain and these values can also vary in time. Spatial position of a point P in the domain Ω can be identified with respect to a spatial reference system (e.g., (x, y, z)). If the position of point P also varies in time, the position of point P is said to be time dependent. Thus, for example, if u is a function of space and time u = u(x, y, z, t). In these notes, we will consider boundary value problems (BVPs) and initial value problems that are formulated by using PDEs. A very general representation of such a problem can be given as follows:

(1.1)

where is a differential operator of independent spatial variables x, y, z and time t, f = f(x, y, z, t) is typically a function that represents the internal effects that act on the system, and τ is the duration of interest.

Figure 1.1 A domain Ω with internal and external boundaries. Note that boundary conditions can be applied on both internal and boundaries.

The dependent variables interact with the outside of the domain Ω through the boundary Γ of the domain, and typically experience changes as a result of the external effects that are imposed on the boundary. These external effects are known as the boundary conditions which depend on the physics of the problem.

The Dirichlet boundary condition represents a prescribed value for a dependent variable,

(1.2)

Here the variable u of the solution domain is prescribed to ub on a segment of the boundary ΓE. In general, this prescribed variable can be a function of time t. The Dirichlet boundary condition is also known as the essential boundary condition.

The von Neumann boundary condition typically describes the external effects that cause a change in the system. Such effects include external forces, heat flow, etc. As we will demonstrate later in the notes, the von Neumann boundary conditions are typically represented as follows:

(1.3)

where is another differential operator, g is a given function, and ΓN represents the segment of the boundary over which the von Neumann boundary condition is applied. The von Neumann boundary condition, also known as the natural boundary condition or the nonessential boundary condition, can also vary in time.

The initial states of the dependent variables are also required to describe the problem properly. For example, the initial values of the variable u and its temporal rate of change u,t and u,tt can prescribed in the solution domain as initial conditions as follows:

(1.4)

where ,t and indicate differentiation with respect to time. Similar initial conditions have to be defined for v and w as well. Number of initial conditions that are sufficient to describe the problem properly depends on the physical system, and this topic is addressed in later chapters.

1.1.2 Boundary value problems

In some problems, only the steady state of the dependent variables is of interest and the temporal variation is neglected (or negligible). Thus, for example, u becomes only a function of the spatial dimensions u = u(x, y, z). A steady state boundary value problem can be formulated by dropping the time dependence as follows:

(1.5)

where for a boundary value problem is a differential operator of the independent spatial variables (x, y, z) and f = f(x, y, z). A steady state boundary value problem is also subject to the Dirichlet and/or von Neumann conditions on the boundary of the domain.

Example 1.1 Equation of motion of a solid bar

a) Derive the equation of motion of an elastic bar in terms of its deflection u(x,t). Initially, assume that the bar has a variable cross-sectional area A(x) and that it is subjected to distributed axial load q(x,t) and a concentrated force F at its free end as shown in Fig. 1.2. Also assume small deflections, linear elastic material behavior with constant elastic modulus E, and constant mass density ρ.

b) Obtain the steady state solution for the case of constant cross-section and zero distributed force.

Figure 1.2 A bar with variable cross-sectional area, subjected to an axially distributed load q ( x ), and the free body diagram of a small segment of this bar at position x .

Solution 1.1a: The solution domain Ω for this problem spans 0 < x < L. The boundaries Γ of the solution domain are located at x = 0 and x = L. Internal forces develop in the bar in response to external loading. The internal normal force N(x) at the cross-section x can be defined as follows:

(1.6)

where the average normal stress is defined as follows:

and where σ is the internal normal stress, A is the cross-sectional area of the bar. The equation of motion of the bar can be obtained by using Newton's second law on a small segment of the bar (Fig. 1.2). The balance of internal and inertial forces gives,

(1.7)

Hooke's law defines the constitutive relationship between the internal stress and strain for linear, elastic materials. For a slender bar, the Hooke's law can be given as follows:

(1.8)

where E is the elastic (Young's) modulus of the material. The strain–displacement, ε–u, relationship is given as follows:

(1.9)

Combining Eqs. (1.6–1.9), we find the internal force resultant as follows:

(1.10)

The equation of motion can then be found by combining Eqs. (1.7c) and (1.10),

(1.11)

This is a PDE that governs the dynamics of axial deflection u(x,t) along the bar. Its solution requires two boundary conditions and two initial conditions. The boundaries of this bar are located at x = 0, L. At the x = L boundary, the force resultant should be equal to the applied load, i.e., N(L) = F. By using Eq. (1.10), this condition can be expressed in terms of the bar deflection. The boundary conditions for this problem then become,

(1.12)

The initial conditions represent the state of deflection and velocity of the entire bar at t = 0. In general, these conditions can be represented as follows:

(1.13)

where u(0)(x) and are known functions.

Solution 1.1b: Let us find the steady state deflection of a bar with constant cross-sectional area, A0, subjected only to a concentrated force F at its free end. Eqs. (1.11) and (1.12) become,

(1.14)

(1.15)

Solution of Eq. (1.14) gives,

(1.16)

where C1 and C2 are integration constants that can be found by using the boundary conditions as follows:

And the variation of u(x) becomes,

(1.17)

Note that this relationship shows that the displacement varies linearly for a bar with cross-section that is fixed on one end and pulled by a force on the other.

1.2 Solution methods

In this work, modeling refers to mathematical formulation of a physical process. This requires background in the related subjects, certain mathematical tools, and experimental observations. In Chapter 2, we present the formulation of models for deformation of elastic solids and transfer and storage of thermal energy in solids and fluids. Solution of the mathematical model can be a challenging task and forms the general background of this work. Analytical solutions which can be expressed as relatively straight forward relationships between the dependent and independent variables exist only for a relatively small number of situations where the geometry and the physical nature of the problem can be simplified. Numerical methods are used otherwise. Among the numerical solution methods for solving PDEs are the finite difference, variational, and finite element methods.

The finite difference method (FDM) is implemented on the differential form of the BVP. The derivative operators of the PDE are approximated by finite difference operators. The solution domain is discretized in to a grid, and the unknowns are the values of the dependent variable at the nodes. The discretized version of the PDE is evaluated at each grid point. This results in a set of algebraic equations which can be represented in matrix form,

(1.18)

where [K] is the stiffness matrix representing the discretized form of the partial derivatives, {D} is the vector of unknown nodal values of the dependent variable, and {R} is the loading vector representing the external effects. The boundary conditions often require specialized treatment of the finite difference operators and modify the [K] matrix. The FDM is effective over relatively simple shapes such as rectangular and cylindrical domains in two-dimensional problems and parallelepiped or spherical domains in three-dimensional problems.

Solution of PDEs by variational methods involves use of a weighted-residual integral, e.g.,

(1.19)

where w is a weight function. One way to distinguish the well-known techniques such as the Rayleigh–Ritz, Galerkin, and collocation methods, which fall under the umbrella of the variational methods, is the choice of the weight functions. Nevertheless, in these traditional methods, an approximation function is chosen for the dependent variable over the entire solution domain. Therefore, traditional variational approaches can be impractical when the shape of the domain is not one of the simple shapes mentioned above. A more complete discussion of the variational methods is given in Chapter 3.

The finite element method (FEM) is based on the variational methods. However, the approximation function for the dependent variable is developed for small regions, known as elements, with regular shapes such as rectangles, triangles, tetrahedra, and hexahedra. As most irregular solution domains can be approximated as a collection of smaller regions (Figs. 1.3 and 1.4) the finite element method presents a powerful technique to solve BVPs over very complex shapes. In the FEM, the domain over which the problem is defined is divided into elements interconnected at nodal points and along edges. An element can take various shapes and can have a number of nodes associated with it. The dependent variable (i.e., the unknown quantity) is determined at the nodes.

Figure 1.3 Definitions of domain and boundary on a two-dimensional object subjected to distributed external forces, f . Finite element analysis requires discretization of the solution domain to subdomains called elements.

Figure 1.4 Elements communicate with the other elements on their boundaries. A line element with three nodes and a 2D element with four nodes are shown in this figure.

FEM overcomes the disadvantage of the traditional variational methods by providing a systematic procedure for the derivation of the approximation functions over subregions of the domain. As the solution domain is represented by a collection of elements, it is relatively straight forward to implement nonhomogeneous, and or discontinuous material properties.

As problems with complex geometry, material properties and boundary conditions cannot be effectively solved with analytical techniques, use of FEM provides must be accompanied with experimental verification. In addition, if possible, analytical forms of the problems should also be sought for simplified conditions as a check to the finite element solution. Finally, the importance of a systematic mesh convergence study for the solution cannot be understated.

Chapter 2

Mathematical modeling of physical systems

2.1 Introduction

The goal of this chapter is to give brief descriptions to modeling of deformation of linear elastic solids and thermal energy transfer and storage in physical systems. More detailed discussion of these topics can be found in the specialized references provided at the end of this chapter. Our goal is to demonstrate how to obtain mathematical models (representations) of physical systems by using the fundamental laws of physics. Thus, we will show that deformation of elastic solids can be described by using Newton's laws of motion. This will result in equations of motion represented as partial differential equations. Vibration of a long and slender bar (Section 2.1), deflection of a general deformable body (Section 2.2), and deflection of beams (Section 2.3) constitute examples of such systems. The principle of conservation of energy will be used to describe effects of heat transfer in a continuum (Section 2.4).

2.2 Governing equations of structural mechanics

When a deformable body is subjected to external effects such as external forces and/or imposed displacements on its boundary, its shape will change and internal forces will develop throughout its volume. The level of deformation for given external effects depends on the material of the deformable body. In this section, the equations of motion for small deflections of linear, elastic materials are presented. In particular, we are interested in small deformations of linear, elastic solids. To this end, following are discussed: i) concepts of external and internal forces and the concept of stress, ii) elastic deformations and the concept of small strain, iii) linear elastic constitutive relations, iv) balance laws, and v) total potential energy of a deformable body.

2.2.1 External forces, internal forces, and stress

External forces acting on continua fall under two categories: body force and surface traction as depicted in Fig. 2.1.

Figure 2.1 Body force , surface traction and surface normal .

The body force, acts on all material points contained in the deformable solid (domain) and it is typically due to presence of an external field that affects the solid at-a-distance. For example, gravitational and magnetic fields result in body forces. The unit of body force is force/volume [F/L³].¹

In contrast, the surface tractions act across the surface (boundary) of the solid body. External tractions are typically due to direct contact of one body with another. Examples include the contact stresses and wind pressure. The unit of surface traction is force/area [F/L²]. Concentrated force or distributed force (over a line) are special cases of surface tractions, with units [F] and force/length [F/L], respectively.

Internal forces develop in a deformable solid when it is subjected to external forces, and help keep the material intact. In general, the internal forces vary from point to point within the body. If we consider an infinitesimally small hexahedral volume ΔV( = dx.dy.dz) (Fig. 2.2) with the corresponding surface areas ΔAx = dy dz, ΔAy = dx dz, ΔAz = dx dy, we will see that each face of this small segment will be subjected to a resultant internal force, , , and . Note that the subscripts x, y, and z in the force vector designations indicate the direction of the unit outward normal of the face on which the force is acting. The same convention is also used in designating the faces of the