Search Finite Elements Website

Minimization of Potential Energy

 home > introductory mechanics > equilibrium
Prev Up Next

Introduction

This page takes us from the governing differential equation(s) of equilibrium to the minimization of potential energy principle that is the heart of finite element analysis.


Equilibrium Review

Recall that the governing differential equations of equilibrium are

\[ {\partial \sigma_{xx} \over \partial x} + {\partial \sigma_{xy} \over \partial y} + {\partial \sigma_{xz} \over \partial z} + f_x = \rho \, a_x \] \[ {\partial \sigma_{yx} \over \partial x} + {\partial \sigma_{yy} \over \partial y} + {\partial \sigma_{yz} \over \partial z} + f_y = \rho \, a_y \] \[ {\partial \sigma_{zx} \over \partial x} + {\partial \sigma_{zy} \over \partial y} + {\partial \sigma_{zz} \over \partial z} + f_z = \rho \, a_z \]
These three equations can be represented by a single matrix or tensor notation equation as

\[ \nabla \cdot \boldsymbol{\sigma} + {\bf f} = \rho \, {\bf a} \qquad \qquad \sigma_{ij,j} + f_i = \rho \, a_i \]
And for the case of static equilibrium, the acceleration vector on the right hand side is set equal to zero, leaving

\[ \nabla \cdot \boldsymbol{\sigma} + {\bf f} = 0 \qquad \qquad \sigma_{ij,j} + f_i = 0 \]
It is this form of the equilibrium equation that is the starting point for the development of the potential energy minimization principle that is at the heart of finite element theory.


Integrating the Equilibrium Equation

Beam
Consider the tapered beam to the right. It is being pulled in tension, so it is logical to expect that the internal force at any \(x\) is \(F\) and that the normal stress is therefore

\[ \sigma_{xx} \, = \, {F \over \pi \, R^2} \, = \, {F \over A_o} e^{2x/L} \]
This is just easy enough that a simple expression for \(\sigma_{xx}\) is possible. But what if the geometry were more complex and multiple loads were present. Then what?!

There are many possible analytical and numerical approaches that could be used to tackle such a problem. But over that past 60+ years, the finite element method has become by far the prefered choice. It is based on minimizing the potential energy of a mechanical system. And yes, minimizations do involve taking derivatives and setting them equal to zero.

In the mechanics world, it is clear that the minimization step should lead one back to the equilibrium equations, which we have already just set equal to zero for the case of static equilibrium. Therefore, we must find some appropriate form of an integral of the equilibrium equations that can then be differentiated to accomplish the minimization step.

Evidently, finding the most useful integral form took many decades. For example, it would be easy to write the following expression

\[ \text{Integral Form} = {1 \over 2} \left( \nabla \cdot \boldsymbol{\sigma} \right)^2 + {\bf f} \cdot \left( \nabla \cdot \boldsymbol{\sigma} \right) \]
and differentiate with respect to \((\nabla \cdot \boldsymbol{\sigma})\) to get back to the equilibrium equation (because it looks like \({1 \over 2} x^2 + f \, x\) ). Except that this is not very helpful, especially if you are also interested in the displacements of the structure. (It turns out that going from stresses and strains back to displacements is not a trivial matter. This is primarily due to the fact that one can have six independent strain components, but only three independent displacement components in a 3-D world.)

Instead, the following approach has been found to be the most useful. Granted, it will probably not be obvious at first why this is the case. Nevertheless, begin by first multiplying the equilibrium equation through by an arbitrary displacement (arbitrary, as in any direction) of differential length. Call this \(\delta {\bf u}\).

\[ ( \nabla \cdot \boldsymbol{\sigma} ) \cdot \delta {\bf u} + {\bf f} \cdot \delta {\bf u} = 0 \]

Clarification of \(\delta {\bf u}\)

It took me a long time to become comfortable with \(\delta {\bf u}\) when first confronted with it in college. It was primarily because I had never seen the \(\delta\) symbol used in such a way before. And that, combined with the idea that \(\delta {\bf u}\) was an arbitrary displacement really bothered me.

I finally realized that \(\delta {\bf u}\) was no different than \(d{\bf u}\), and that the \(\delta\) was simply being used to minimize confusion with other differentials, \(dV\) for volume and \(dS\) for surface area, that were coming up later in the derivation.

So if \(\delta {\bf u}\) bothers you, then replace it with \(d{\bf u}\), or even \(d{\bf u} \over dt\), which would represent a velocity vector. If you do use \(d{\bf u} \over dt\), the \(dt\) will float along throughout the entire derivation, and could then be canceled out at the end.
The next step is to recognize the following identity for differentiation of products

\[ \nabla \cdot ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) = \underbrace{( \nabla \cdot \boldsymbol{\sigma} ) \cdot \delta {\bf u}}_ {\matrix{\text{Term in} \\ \text{equilibrium} \\ \text{equation above}}} + \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) \]
produces a term that is in the equilibrium equation. Solve for this term

\[ ( \nabla \cdot \boldsymbol{\sigma} ) \cdot \delta {\bf u} = \nabla \cdot ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) - \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) \]
and replace it in the equation, giving

\[ \nabla \cdot ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) - \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) + {\bf f} \cdot \delta {\bf u} = 0 \]
It is common at this point to do two minor things: (i) shuffle the order of the terms, and(ii) multiply through by \(-1\). Doing so gives

\[ \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) - {\bf f} \cdot \delta {\bf u} - \nabla \cdot ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) = 0 \]
Next, integrate each term over the volume of the object

\[ \int \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) \, dV - \int {\bf f} \cdot \delta {\bf u} \, dV - \int \nabla \cdot ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) \, dV = 0 \]
and recognize that the divergence theorem can be applied to the third volume integral. Doing so will transform the term into a surface integral.

Looking at the divergence theorem for the third term:

\[ \int \nabla \cdot ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) \, dV = \int ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) \cdot {\bf n} \, dS \]
where \({\bf n}\) is the outward pointing unit normal for each differential element of surface area.

Applying the divergence theorem to the third term gives

\[ \int \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) \, dV - \int {\bf f} \cdot \delta {\bf u} \, dV - \int ( \boldsymbol{\sigma} \cdot \delta {\bf u} ) \cdot {\bf n} \, dS = 0 \]
But now, the third term has a \(\boldsymbol{\sigma}\) and an \({\bf n}\) in it. Recall that the dot product of the two gives the traction vector, \({\bf T}\). This is the case here, even though there is a \(\delta{\bf u}\) in between them. Replacing \(\boldsymbol{\sigma} \cdot {\bf n}\) with \({\bf T}\) gives

\[ \int \boldsymbol{\sigma} : \nabla (\delta {\bf u} ) \, dV - \int {\bf f} \cdot \delta {\bf u} \, dV - \int {\bf T} \cdot \delta {\bf u} \, dS = 0 \]
In the first term, the \(\nabla\) and \(\delta\) can be swapped to give

\[ \int \boldsymbol{\sigma} : \delta (\nabla {\bf u} ) \, dV - \int {\bf f} \cdot \delta {\bf u} \, dV - \int {\bf T} \cdot \delta {\bf u} \, dS = 0 \]
And even though \(\nabla {\bf u}\) is not \(\boldsymbol{\epsilon}\), it turns out that \( \boldsymbol{\sigma} : \delta (\nabla {\bf u}) \) is equal to \( \boldsymbol{\sigma} : \delta \boldsymbol{\epsilon} \) because of the symmetry of \( \boldsymbol{\sigma} \). So the equation becomes

\[ \int \boldsymbol{\sigma} : \delta \boldsymbol{\epsilon} \, dV - \int {\bf f} \cdot \delta {\bf u} \, dV - \int {\bf T} \cdot \delta {\bf u} \, dS = 0 \]
in which \(\delta \boldsymbol{\epsilon}\) can be replaced by \(\delta (\nabla {\bf u})\).

So what's the point? It is that the result here suggests what should be the quantity that, when minimized, leads to the equilibrium equation. The \( \boldsymbol{\sigma} : \delta \boldsymbol{\epsilon} \) term suggests that the term should be the strain energy density, \(w\), integrated over the volume. The \( {\bf f} \cdot \delta {\bf u} \) term suggests that this term should be the work done by internal forces, \({\bf f} \cdot {\bf u}\). And the \( {\bf T} \cdot \delta {\bf u} \) term suggests that this one should be the work done by external traction forces, \({\bf T} \cdot {\bf u}\).


Potential Energy Minimization

The quantity being minimized is called Potential Energy and is represented by, \(\pi({\bf u})\). It is equal to

\[ \pi({\bf u}) = \int w \, dV - \int {\bf f} \cdot {\bf u} \, dV - \int {\bf T} \cdot {\bf u} \, dS \]
and minimizing \(\pi({\bf u})\) with respect to the displacements, \({\bf u}\), leads to a solution that satisfies the equilibrium equation. For a linear elastic material, the strain energy density equals

\[ w = {1 \over 2} \boldsymbol{\sigma} : \boldsymbol{\epsilon} \]
so one could write

\[ \pi({\bf u}) = \int {1 \over 2} \boldsymbol{\sigma} : \boldsymbol{\epsilon} \, dV - \int {\bf f} \cdot {\bf u} \, dV - \int {\bf T} \cdot {\bf u} \, dS \]
Or one could replace \(\boldsymbol{\sigma}\) with \({\bf C} : \boldsymbol{\epsilon}\), giving yet another form

\[ \pi({\bf u}) = \int {1 \over 2} ({\bf C} : \boldsymbol{\epsilon}) : \boldsymbol{\epsilon} \, dV - \int {\bf f} \cdot {\bf u} \, dV - \int {\bf T} \cdot {\bf u} \, dS \]
And finally, due to the symmetry of the terms, one can insert \(\nabla {\bf u}\) for each of the \(\boldsymbol{\epsilon}\)'s.

\[ \pi({\bf u}) = \int {1 \over 2} ({\bf C} : \nabla {\bf u}) : \nabla {\bf u} \, dV - \int {\bf f} \cdot {\bf u} \, dV - \int {\bf T} \cdot {\bf u} \, dS \]
A great of advantage of all these expressions is that the quantity being minimized, the potential energy, is a function of displacements. So minimizing with respect to displacements will produce expressions for the displacements of the object in question. It is then possible, straight-forward in fact, to compute the strains from the displacements and then the stresses from the strains.

So lets return to the tapered beam and apply the potential energy minimization principle to the problem. We will start with the Rayleigh-Ritz approach, which is to guess a form of the displacement field solution, with one (or more) unknown constants. So lets guess that the displacements in the x-direction are

\[ u_x = C \, x \]
where \(C\) is an unknown constant, not to be confused with the stiffness tensor, \({\bf C}\), which is bold. The first step is to compute the strain energy density. The easiest approach this time is to take the derivative of the assumed displacement field with respect to \(x\) to get the strain. Doing so gives \(\epsilon = C\). And the stress is \(\sigma = E \, \epsilon = E \, C\). So the strain energy density is \({1 \over 2} E \, C^2\).

There are no internal forces, so the \(\int {\bf f} \cdot {\bf u} \, dV\) term equals zero.

And the only external force is \(F\) at \(x = L\), so \(\int {\bf T} \cdot {\bf u} \, dS\) is simply \(F \, C \, L\). Combining all this gives

\[ \pi({\bf u}) = \int {1 \over 2} E \, C^2 \, dV - F \, C \, L \]
Integrating over the volume is tricky because the cross-section varies. Fortunately, the \({1 \over 2} E \, C^2\) term is constant. The volume is \(A_o \left(L \over 2\right) \left[ 1 - e^{-2}\right]\).

So the total potential energy is

\[ \pi = {1 \over 4} A_o L \left[ 1 - e^{-2} \right] E \, C^2 - F \, C \, L \]
So now the undetermined constant is \(C\), and the way to minimize \(\pi\) is to differentiate the expression with respect to \(C\) and set the result equal to zero.