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