This information was made from notes from several sources.  The main one was from a Finite Element Course I took at the University of Iowa from Dr. Bhatti.

where U is the strain energy and  is the virtual work done by the system in deformation.  The strain energy is computed as a volume integral of the strain density .

        (compare to )

where

resulting in

where

For small deformations the strains are...

For large deformations, the strains are...

Later on when discussing interpolation functions, (ξ, ψ, ζ) are replaced with (u, v, w).  From the above equations, the following derivatives are needed:

for small deformations, the strains are


 

Where  for an eight node Hexahedral (Brick) element, assume a trial solution u(x,y,z), v(x,y,z), w(x,y,z) as a function of interpolation function and the degrees of freedom:

where  is the trial solution (u, v, and w), N is the shape functions as defined by the parent element, and  is the degree of freedom vector (ui, vi, wi).
(See Figure 4)
The shape functions are generalized as

where r, s, t is the parent element coordinate system that the element will be mapped too.  (See figure 5)  The coordinates of r, s, and t for nodal coordinates are...
TABLE  (based on figure 5)

1
2
3
4
5
6
7
8

                   resulting in 
 

For a twenty node isoparametric Solid Element, all definitions are just expanded to 20 instead of 8 with the following r,s and t coordinates and shapefunctions.
               resulting in 

For space, the rest of the building description will be done for the 8 node element.

Nodal coordinates before deformation are in the x, y, z coor. system.   In order to transform to the u, v, w coordinate system, we must define nodal coordinates in terms of r, s and t.  To do this we need trasformation functions or isoparametric mapping functions.

So, given an actuall element with coordinates (xi, yi, zi) for each node (for eight noded element)

then the interpolation  functions for the actuall element represented in the parent element coordinates (r, s, and t) is

Now because the interpolation functions are in terms of the parent element and eventually strains are needed in the original coordinate system, a transformation is needed to switch back and forth.

For the derivatives of u...

where J is the Jacobian matrix.  The derivatives of u (with respect to x, y, and z) can be computed by inverting the Jacobian matrix J, resulting in (with substitution from the approximate solution above),

;

Likewise for v and w,

;     and ;

with the reverse as

;
;

 

Now to get all of these derivatives for the jacobian and to check the determinant over the interval to see if the mapping was good.  Ideally, one would check the mapping by making sure that the determinant of the Jacobian is not zero anywhere over the element.  But for the 3-D case, this proves to be very difficult so common convention is to check whether the determinant of the Jacobian at all of the gauss points (of integration) are the same sign and if it is, then assuming it is not zero anywhere.

Remember that...

substituting with the definition gained above, results in

so

where K is the elemtne stiffness matrix

is the virtual work done by the system in the deformation.  The virtual work done by the system against the forces can be written as:

or 
or 
or 

where Fx, Fy and Fz are the body forces per unit volume and Px, Py, and Pz are surface or boundary forces per unit area.  Let the body forces per unit volume be the inertial forces:

, , .

Note that if there is no time dependence (static case) then there are no body forces per unit volume.  The boundy forces per unit area only happen at the boundary and will be primarily from the air flow and air pressure.  F is the body force vector and P is the equivalent nodal load vector.

;      and 
The  surface or boundary forces per area are external to the element and for the element here  would be of the form...

These componants are literally the stresses in each direction inposed by external forces.

So, the potential energy is

For FEM purposes, we want the variation of the potential energy, where each term is each element is zero with respect to its modal coordinates.   Let each integral represented by .

The derivatives with respect to individual nodal points in u, v, and w will then be taken (minimization) and set equal to zero as stated in the equations above.

If K is symmetric, then the minimization is easy with one large matrix equation:

Viscous damping effects can be included by replacing μ in the stiffness matrix K and replacing it with  μ+η ∂/∂t.  Let this new addition be the damping matrix D.  Where .

For static problems, this reduces to

For time dependent problems, the time derivatives can be split into descreat steps by a finite difference scheme (where n denotes the time step).  The time derivatives can be approximated by first order centeral differences as follows...

The stiffness can be replaced (Crank-Nickolson technique, Nickell, 1973) by the average at times n and n+1, or

With these substitutions, the equation above can be simplified to

 

???
To account for internal forces, let K be replaced by  where ...
I
f is the force vector from muscle contractions and I is the identity matrix.

so: 

These integrals will be approximated using Product-Gaussian integration.  So, first the integrals...

where ,  are Gause points, and l, m, and n are the number of Gause points in the r, s, and t directions respectively.  w represents Gause weights in each direction and  is the value of the integrand at the point .  (pg 9-21 has gauss weights and points).  For 2 point  Gaussian Quadrature, the weights are all unitary (1) and the positions between -1 and 1 in each direction comes from the zeros of the  2nd Legendre polynomial which is  or .

So now to evaluate the integrals


 
 
 
 
 
 
 
 
 

Convergence Requirements  (B5 page 126)

The patch test is a simple test that can be performed numerically, checking the validity of an elements formulation and implementation.  A small patch of elements is made with at least one node within the pathc so that the node is shared by two or more elements.  Also, one or more interelement boundaries must exist.  Execute a standard solution and examine the computered stresses.  computed stresses must agree with exact stresses for the physical probelm model allowing only for numerical noise being the difference.