Engineering treatment of the finite element method

:"This is a draft of a new explanation as suggested on ."

The finite element method (FEM) is a technique for finding approximate solutions to differential equations that is particularly useful in engineering. As of 2005, FEM is the primary analysis technique for computer modeling of mechanical systems as found in structural mechanics.

Motivation

FEM is related to linear algebra approaches for solving the forces and displacements of a truss. In solving a truss, each piece can be thought of as a linear spring. We then apply boundary conditions and a force and wish to find the displacement of each node. This can be solved by constructing a stiffness matrix which combines the stiffness of each element in the system and the nodes each connects. The equation "Kd" = "F", where "K" is the stiffness matrix, "d" is the displacement vector with an entry for each node in each direction, and "F" is a vector of forces at each node in each direction. Clearly "d" can be solved by inverting "K" to get "d" = "K" ⁻¹"F".

This solution is elegant, but in real-world applications, not all structures can be represented by pinned elastic beams. In a rough sense, the finite element method generalizes this approach to be applicable to more-general geometries and governing equations allowing us to use FEM to model everything form elastostatics to turbulent flow.

Approximation

Suppose we would like to solve a solid-mechanics problem using linear elastostatics. That is, we have a system for which

:$sigma\_\{ij,j\}+f\_i=0$ on the boundary, Ω

:$u\_i=g\_i$ on the region we have prescribed displacement boundary

conditions, $Gamma\_\{gi\}$

:$sigma\_\{ij\}n\_j=h\_i$ on the traction boundary, $Gamma\_\{h\_i\}$

where $sigma\_\{ij\}=c\_\{ijkl\}epsilon\{kl\}=c\_\{ijkl\}\; u\_\{(ij)\}$. (See notation section, below.)

We would like to find a displacement field, "u" that satisfieds these equations. Experience with differential equations on complicated geometry tells us that we cannot expect to find a closed form for "u"; the next best thing is to find some approximation, "u"^"h".

We would like "u"^"h" to be a "finite-dimensional function" (see vector space). In particular, we would like to construct it from a linear combination of some predefined basis functions, or shape functions, vis.,

:

where the $N\_A$s are the shape functions and the $d\_A$s are the coefficients. We can choose whatever shape functions we would like, however that choise will affect the quality of our solution.

Now, given a set of shape functions, we need only solve for the $d\_A$s.

Projection

Recall from linear algebra, that a dot product projects one vector onto another. Similarly, an "n"-dimensional vector can be projected onto an "n"-dimensional vector space using an "m"-by-"n" matrix, vis.,

:

That is, the matrix on the left projects the vector from three- to two-dimensional space. If you recall, the product of a matrix by a vector is computed by taking the inner product ("dot product") of each row of the matrix with that vector to produce the corresponding row (element) of the result. Consider the possibility that, if we were to define some other type of dot product, we could define matrix–vector multiplication in terms of that.

In our problem we would like to project the unknown "u" from an infinite-dimensional space onto our solution space. To this end, we would like to project it onto each of our basis functions; that dot product will give us the $d\_A$s.

Inner products

To define an inner product on a function, consider the definition of an inner product on vectors:

:$ucdot\; v\; =\; sum\_\{i=1\}^n\; u\_i\; v\_i.$

Note that a function, $f(x)$ is an infinite-dimensional vector in the sense that it is indexed by "x". So extrapolating the summation into a continuum, we could imagine an inner product, $a(f,g)$, of functions, $f(x)$ and $g(x)$ being defined as

:$a(f,g)\; =\; int\_\{-infty\}^\{infty\}\; f(x)g(x),dx.$

We will define other inner products, but they will all have a similar form.

Energy

So now we need to define an inner product on "u". While we have freedom in the matter, it will affect our final answer, so we would like to choose something sensible.

We will define the following inner product to project one function onto another.:$a(w,u)=int\_Omega\; w\_\{(i,j)\}c\_\{ijkl\}u\_\{(i,j)\},dOmega$The integral can be considered physically. The function "w" can be considered a variation of position, or "virtual displacement" and so $w\_\{(i,j)\}$ can be thought of as "virtual strain". In other words, $a(w,u)$ represents the strain energy done by the body were it to be deformed from "u" to "u" + "w".

Now, we would like to compute $d\_A\; =\; a(u,\; u^h)$; unfortunately, we have already given up on finding "u" itself. However, we "can" find this in terms of things we know. We noted that $a(w,u)$ defines a virtual strain energy. By conservation, the internal change in energy must be balanced by energy that went into and came out of the body, that is, the virtual work done by body forces (such as gravity) and the virtual work done by surface tractions. That is,

:$a(w,u)=int\_Omega\; w\_\{(i,j)\}c\_\{ijkl\}u\_\{(i,j)\},dOmega\; =\; int\_Omega\; w\_i\; f\_i\; dOmega\; +\; sum\_\{i=1\}^\{n\_\{sd\; left(\; int\_\{Gamma\_\{h\_i\; w\_i\; h\_i\; ,\; dGamma\; ight)$

For convenience of notation, we will define $(w,f)$ to be the body-force energy term and $(w,h)\_Gamma$ to be the traction energy term. That is, we have

:$a(w,u)\; =\; (w,f)\; +\; (w,\; h)\_Gamma.$

Now we can let "w" equal each of the shape functions in turn and use this relation to find the projection of "u" onto each of shape functions......

A single element

Thus far, we have not discussed "elements". While the practical effectiveness of this method depends on descretizing space into elements, we could still get "correct" answers by using shape functions which are non-zero across almost all of the domain.

In practical applications, the finite element method descretizes space to simplify the matrix inversion problem. By choosing sets of shape functions such that only a few are non-zero over any particular element of the domain, we will get a sparse matrix which is much less expensive to invert.

For the moment, consider the problem of constructing a single element which in some way will act like a spring to approximate the behavior of an elastic solid.:"Start a simple 2D elastostatics example here. (1D linear elasticity is boring and anything other than elasticity confuses the analogy with the truss analysis.)"Across the element, we will approximate the displacement field, "u", by a linear combination of a small number of shape functions.

:"etc."

Relation to variational methods for finding the optimal approximation...:"etc."

Combining elements

Now that we have a method for creating an element stiffness matrix we can create the global stiffness matrix in much the same way it was created for the case of the truss.

Notation

Plain subscript is an index. That is,:$u\_x$is the component of "u" in the "x" direction.

Subscript with a comma is for partial differentiation, so:$u\_\{x,y\}\; =\; frac\{partial\; u\_x\}\{partial\; y\}$

Parenthesis defines the symmetric derivative,:$f\_\{(x,y)\}=frac\{1\}\{2\}\; (f\_\{,x\}\; +\; f\_\{,y\}).$This is part of the definition of strain.

A domain is indicated by Ω and its boundary by Γ.