Finite element method

2D FEM solution for a magnetostatic configuration (lines denote the direction of calculated flux density and colour - its magnitude)

2D mesh for the image above (mesh is denser around the object of interest)

Finite element methods (FEM) (sometimes referred to as finite element analysis) are numerical techniques used for finding approximate solutions of partial differential equations (PDE) as well as of integral equations. The solution approach is based either on reducing the differential equations to a set of linear equations (steady state problems) or into an approximating system of ordinary differential equations (time-dependent problems).

FEM use a more general approach than the finite difference method or the finite volume method, both of which are types of finite element methods, by using a larger number of function spaces. Its generality makes it suitable to use with complex domains, such as cars and oil pipelines, or moving domains, such as melting wax. Extension to the basic method allow for other challenging problems, such as modeling discontinuous solutions with XFEM.

FEM are used for computer simulations in both research science and engineering applications. There are many finite element software packages, both free and proprietary.

History

The finite-element method originated from the needs for solving complex elasticity, structural analysis problems in civil and aeronautical engineering. Its development can be traced back to the work by Alexander Hrennikoff (1941) and Richard Courant (1942). While the approaches used by these pioneers are dramatically different, they share one essential characteristic: mesh discretization of a continuous domain into a set of discrete sub-domains, usually called elements.

Hrennikoff's work discretizes the domain by using a lattice analogy while Courant's approach divides the domain into finite triangular subregions for solution of second order elliptic partial differential equations (PDEs) that arise from the problem of torsion of a cylinder . Courant's contribution was evolutionary, drawing on a large body of earlier results for PDEs developed by Rayleigh, Ritz, and Galerkin.

Development of the finite element method began in earnest in the middle to late 1950s for airframe and structural analysis and gathered momentum at the University of Stuttgart through the work of John Argyris and at Berkeley through the work of Ray W. Clough in the 1960s for use in civil engineering. By late 1950s, the key concepts of stiffness matrix and element assembly existed essentially in the form used today and NASA issued request for proposals for the development of the finite element software NASTRAN in 1965. The method was provided with a rigorous mathematical foundation in 1973 with the publication of Strang and Fix's An Analysis of The Finite Element Method, and has since been generalized into a branch of applied mathematics for numerical modeling of physical systems in a wide variety of engineering disciplines, e.g., electromagnetism and fluid dynamics.

Application

Visualization of how a car deforms in an asymmetrical crash using finite element analysis.

A variety of specializations under the umbrella of the mechanical engineering discipline (such as aeronautical, biomechanical, and automotive industries) commonly use integrated FEM in design and development of their products. Several modern FEM packages include specific components such as thermal, electromagnetic, fluid, and structural working environments. In a structural simulation, FEM helps tremendously in producing stiffness and strength visualizations and also in minimizing weight, materials, and costs. FEM allows detailed visualization of where structures bend or twist, and indicates the distribution of stresses and displacements. FEM software provides a wide range of simulation options for controlling the complexity of both modeling and analysis of a system. Similarly, the desired level of accuracy required and associated computational time requirements can be managed simultaneously to address most engineering applications. FEM allows entire designs to be constructed, refined, and optimized before the design is manufactured. This powerful design tool has significantly improved both the standard of engineering designs and the methodology of the design process in many industrial applications. The introduction of FEM has substantially decreased the time to take products from concept to the production line. It is primarily through improved initial prototype designs using FEM that testing and development have been accelerated. In summary, benefits of FEM include increased accuracy, enhanced design and better insight into critical design parameters, virtual prototyping, fewer hardware prototypes, a faster and less expensive design cycle, increased productivity, and increased revenue.

Outline of general mathematical framework

A function in H¹₀, with zero values at the endpoints (blue), and a piecewise linear approximation (red).

Basis functions v_k (blue) and a linear combination of them, which is piecewise linear (red).

A piecewise linear function in two dimensions.

The basic idea is to replace the infinite dimensional linear problem:

Find $u\in H_{0}^{1}$ such that

$\forall v\in H_{0}^{1},\;a(u,v)=\int fv$

with a finite dimensional version:

Find $u\in V$ such that

$\forall v\in V,\;a(u,v)=\int fv$

where $V$ is a finite dimensional subspace of the Sobolev space $H_{0}^{1}$. There are many possible choices for $V$, FEM limits these choices by using a finite element defined below.

In general, the finite element method is characterized by the following process.

• Formulate a weak or variational form of the PDE.
• Choose a finite element.
• Create a mesh for the continuous domain ($\Omega$) using the reference shape of choosen finite element. (not described here)
• Assemble a set of linear equations by applying the weak form on each element in the mesh.

Weak or variational form

The variational form, sometimes referred to as the weak form, is obtained by integration of the strong form of PDE with any given function in the approximating subspace. If able it is usual to integrate by parts using a Green's theorem. For example, finding $u$ in the Poisson equation:

(1) $-\Delta u=f$

where $f$ is given. One defines the bilinear form:

(2) $a(u,v):=\int _{\Omega }f\cdot v\,dx=-\int _{\Omega }\Delta u\cdot v\,dx=\int _{\Omega }\nabla u\cdot \nabla v\,dx$

where $v$ is any function in the approximating subspace, $\nabla$ denotes the gradient and $\cdot$ denotes the dot product. Using this bilinear form one will create an equation for each basis function of $v$ by substituting each basis function for $u$ and set equal to the linear form:

(3) $L(v):=\int _{\Omega }f\cdot v\,dx.$

Choosing a finite element

The finite element is defined as):

• a reference shape (such as an interval, triangle, square, tetrahedron, or other suitable polygon),
• a set of basis functions for representing $V$(piecewise polynomial basis functions are the most common), and
• a set of evaluation node on the reference shape where the basis functions are evaluated.

It is common to use Lagrange polynomials with the evaluation nodes placed at the vertices of the polygon. In the one-dimensional case, one uses the evaluation nodes at $x_{k}$, using piecewise linear functions $v_{k}$ in $V$ whose value is $1$ at $x_{k}$ and zero at every $x_{j},\;j\neq k$, i.e.,

$v_{k}(x)={\begin{cases}{x-x_{k-1} \over x_{k}\,-x_{k-1}}&{\mbox{ if }}x\in [x_{k-1},x_{k}],\\{x_{k+1}\,-x \over x_{k+1}\,-x_{k}}&{\mbox{ if }}x\in [x_{k},x_{k+1}],\\0&{\mbox{ otherwise}},\end{cases}}$

for $k=1,...,n$; this basis is a shifted and scaled tent function. Similarly, for the two-dimensional case, one uses one basis function $v_{k}$ per vertex $x_{k}$ of the discretization of the planar region $\Omega$. The function $v_{k}$ is the unique function of $V$ whose value is $1$ at $x_{k}$ and zero at every $x_{j},\;j\neq k$.

The primary advantage of this choice of basis is that the inner products

$(v_{j},v_{k})=\int _{0}^{1}v_{j}v_{k}\,dx$

and

$a(v_{j},v_{k})=\int _{0}^{1}v_{j}'v_{k}'\,dx$

will be zero for almost all $j,k$. In the one dimensional case, the support of $v_{k}$ is the interval $[x_{k-1},x_{k+1}]$. Hence, the integrands of $(v_{j},v_{k})$ and $a(v_{j},v_{k})$ are identically zero whenever $|j-k|>1$.

Similarly, in the planar case, if $x_{j}$ and $x_{k}$ do not share an edge of the triangulation, then the integrals

$\int _{\Omega }v_{j}v_{k}\,ds$

and

$\int _{\Omega }\nabla v_{j}\cdot \nabla v_{k}\,ds$

are both zero.

Linear equation assembly

Solving the two-dimensional problem $u_{xx}+u_{yy}=-4$ in the disk centered at the origin and radius 1, with zero boundary conditions.
(a) The triangulation.

(b) The sparse matrix L of the discretized linear system.

(c) The computed solution, $u(x,y)=1-x^{2}-y^{2}.$

The process of formulating the linear equations is called assembly. Some methods will explicitly form a matrix while others will only assemble elements as each cell of the mesh is iterated. The different choices in assembly are implementation details that often depend on the specific problem being solved.

To assemble the linear equations, one defines $u(x)=\sum _{k=1}^{n}u_{k}\cdot v_{k}(x_{k})$ and $f(x)=\sum _{k=1}^{n}f_{k}\cdot v_{k}(x_{k})$, where $u_{k}$ and $f_{k}$ are the coefficients of $u$ and $f$ at the points $x_{k}$, and $v_{k}$ is the basis function that equals one at $x_{k}$. The problem becomes

$a(u,v)=a(\sum _{k=1}^{n}u_{k}v_{k},\sum _{j=1}^{n}v_{j})=\int _{\Omega }(\sum _{k=1}^{n}f_{k}v_{k})\cdot (\sum _{j=1}^{n}v_{j})\,dx$

which can be written on column vector form:

$\sum _{k=1}^{n}u_{k}a(v_{k},v_{j})=\sum _{k=1}^{n}f_{k}\int v_{k}v_{j}$ for $j=1,...,n$.

Denote by $\mathbf {u}$ and $\mathbf {f}$ the column vectors $(u_{1},...,u_{n})^{t}$ and $(f_{1},...,f_{n})^{t}$, and $K=(K_{ij})$ (known as the stiffness matrix) and $M=(M_{ij})$ (known as the mass matrix) be matrices whose entries are $K_{ij}=a(v_{i},v_{j})$ and $M_{ij}=\int v_{i}v_{j}$ then we may rephrase as

$K\mathbf {u} =M\mathbf {f}$.

As we have discussed before, most of the entries of $L$ and $M$ are zero because the basis functions $v_{k}$ have small support. So we now have to solve a linear system in the unknown $\mathbf {u}$ where most of the entries of the matrix $L$, which we need to invert, are zero.

Such matrices are known as sparse matrices, and there are efficient solvers for such problems (much more efficient than actually inverting the matrix.) In addition, $L$ is symmetric and positive definite, so a technique such as the conjugate gradient method is favored. For problems that are not too large, sparse LU decompositions and Cholesky decompositions still work well. For instance, Matlab's backslash operator (which uses sparse LU, sparse Cholesky, and other factorization methods) can be sufficient for meshes with a hundred thousand vertices.

Additional concerns

More advanced implementations (adaptive finite element methods) utilize a method to assess the quality of the results (based on error estimation theory) and modify the mesh during the solution aiming to achieve approximate solution within some bounds from the 'exact' solution of the continuum problem. Mesh adaptivity may utilize various techniques, the most popular are:

• moving nodes (r-adaptivity)
• refining (and unrefining) elements (h-adaptivity)
• changing order of base functions (p-adaptivity)
• combinations of the above (e.g. hp-adaptivity)

Various types of finite element methods

Generalized finite element method

The Generalized Finite Element Method (GFEM) uses local spaces consisting of functions, not necessarily polynomials, that reflect the available information on the unknown solution and thus ensure good local approximation. Then a partition of unity is used to “paste” these spaces together to form the approximating subspace. The effectiveness of GFEM has been shown when applied to problems with domains having complicated boundaries, problems with micro-scales, and problems with boundary layers.

XFEM

Main article: XFEM

Spectral methods

Main article: Spectral methods

Meshfree methods

Main article: Meshfree methods

See also

• Boundary element method
• Discrete element method
• Finite element machine
• Finite element method in structural mechanics
• Galerkin method
• Multiphysics
• Patch test
• Rayleigh-Ritz method
• List of finite element software packages

External links

• NAFEMS -- The International Association for the Engineering Analysis Community
• IFER -- Internet Finite Element Resources - an annotated list of FEA links and programs
• Workshop "The Finite Element Method in Biomedical Engineering, Biomechanics and Related Fields"
• Finite Element Analysis Resources- Finite Element news, articles and tips
• Finite Element Books- books bibliography
• Mathematics of the Finite Element Method
• Finite Element Methods for Partial Differential Equations
• FEM AVI-gallery at CompMechLab site, St.Petersburg State Polytechnical University, Russia
• Illustrative Examples of Implementing the Finite Element Method in Software CFDLab, University of Texas at Austin, USA