Engineering treatment of the finite element method

This is a draft of a new explanation as suggested on Talk: Finite element method.

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 and continuum mechanics.

Motivation

Template:Sec-stub FEM is related to linear algebra approaches for solving the forces and displacements of a truss. Each element can be thought of as a 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 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. 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)}$. [Insert explanation of notation.]

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 basis functions, or shape functions, vis.,

${\mathbf {u}}^{h}=\sum _{A=1}^{n}N_{A}d_{A}$

where the $N_{A}$s are the shape functions and the $d_{A}$s are the coefficients. Now we need only solve for the $d_{A}$s.

Projection asside

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.,

${\begin{bmatrix}1&0&0\\0&1&0\end{bmatrix}}{\begin{bmatrix}x\\y\\z\end{bmatrix}}={\begin{bmatrix}x\\y\end{bmatrix}}$

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:

$u\cdot 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)}\,d\Omega$

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 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)}\,d\Omega =\int _{\Omega }w_{i}f_{i}d\Omega +\sum _{i=1}^{n_{sd}}\left(\int _{\Gamma _{h_{i}}}w_{i}h_{i}\,d\Gamma \right)$

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

In general, FEM descretizes space to simplify the problem into a matrix inversion problem. 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.