Finite element analysis for beams

So far we have looked at finite elements with interpolation functions that satisfied (C^0)continuity. However, beams are structural elements that require a different treatment. This topic though strictly speaking not essential for study in our introductory course, is important in the sense that it tells us that finite element analysis requires different approaches and is not always amenable through just one method. 

What are beam (or shell) elements:
The structures that are thin relative to their major dimensions can be modeled as beams or shells. 
1. beams and columns of high-rise structures
2. the sheet metal and frame of various vehicles are modeled using shell elements
3. the hull of a ship of fuselage of aircraft, where shell elements are used. 

Shell elements are more involved than beam elements. However, the idea behind beams and shell elements is fundamentally very similar.

Beam Kinematics

Nomenclature for the beam. The midplane for the beam is called as the neutral axis. The deformation of the Euler-Bernoulli beam happen as per this figure


Use the standard kinematics and equilibrium as below 

		Sign-convention for equilibrium
	
Note that the internal forces are:
	i. Internal moment: m(x)
	ii. Internal shear force: s(x)
and as described in the strength of material, e.g., Popov or Timoshenko and Gere, we get the following equations:

1. By summing all the vertical forces:

dividing this equation with Δx and taking the limit Δ x -> 0, we get

2. By considering moment equilibrium about the point x = y = 0, we get

and dividing  by Δ x and taking the limit Δ x -> 0 we get


3. Combining the two equations (1) and (2) we get


4. Now noting that using the basics of Euler-Bernoulli theory, where the internal moment m = E  I  (((d^2)u )/(d (x^2))), we get the following equation

Boundary Conditions

Differential equation is:

Boundary Conditions:
This governing equation is a fourth order differential equation and the boundary conditions for this equation are:
	
Note: 
	1. the quantity n is external normal to the beam along the length of the beam. For example, in the Figure of the beam above with the loading, at the left end of the beam n = -1 and at the right end of the beam n = +1. This is not the same n as shown in the Figure above regarding nomenclature of the beam. 
	2. Our sign convention is that (Overscript[m, -]) is positive anticlockwise. Hence the inclusion of the normal n as discussed above ensures that moment is correctly defined  as in the sign-convention figure. 
	3. Also the shear force (Overscript[s, -]) is positive upwards, hence by multiplying with n we ensure that we obey the sign convention for shear force as described in the sign-convention figure. 
	
At any end of the beam you could have boundary conditions in the following common combinations

Like in the case of elasticity at any point we cannot have energy conjugates: m and θ  or s and (Subscript[u, y])being specified simultaneously. Mathematically, this means:

Strong form to weak form

In the case of the beams, since the governing equation is fourth order, the weak form will look at a bit different. However, we re-write the fourth order equation in terms of the second order equation in terms of the bending moment. 

and from here:

Integrating this equation by parts we get

where like in the previous cases, we note that the weight w should vanish on the essential boundaries (Subscript[Γ, u]) to give

We now do another integration parts on the term involving first derivatives of both m and w

In the above equation, we note that ((d w)/(d x)) should vanish on the essential boundary (Subscript[Γ, θ]) -- here we are encountering something completely new than anything before that. 

Combining all the relations above together, we get the following final expression for the weak form

This is actually the statement of principle of virtual work. The left hand side of the equation tells us the internal virtual work done and the right side is the external virtual work. In this case w = (Subscript[δu, y]) can be thought of to be the virtual displacement of the beam. Now, by noting that:
m = E I (((d^2)(Subscript[u, y]))/(d (x^2)))
the weak for will look as

As before we have symmetry in the way w and (Subscript[u, y])appear in the weak form derivatives. 

We now have to carefully look into the structure of U and (Subscript[U, 0]), respectively, the spaces corresponding to (Subscript[u, y]) and w respectively. Since we have second derivative of both these functions appearing in the integral we can no longer have (C^0) functions approximating them. We need (C^(1 ))functions, i.e. functions without jumps and kinks. At a more physical level, the presence of kink will make the deformations incompatible (something that was desirable in linear elasticity). Since the integrals need to finite, (C^1)functions that are also integrable are called as (H^2)functions. 

	Kink makes the displacement loose compatibility. 
	
	So we define U and (Subscript[U, 0])are spaces such that
	
	and our goal is to 
	
	find (Subscript[u, y])∈ U such that the weak form
	
	holds for ∀ w(x) ∈ (Subscript[U, 0]). 

It could also be shown that weak form implies strong form.

Finite element discretization

The Lagrange polynomials that we had discussed earlier are no longer sufficient since they have kinks in them. The simplest functions with (C^1) are the Hermite polynomials. Now, since we need to maintain both slope and displacement continuity when we more from one element to the other we need to express the interpolation in terms of both displacements and rotations. Hence the degree of freedom at every node of the beam element are (Subscript[u, y]) and θ as shown below


and the corresponding nodal forces are work conjugates. 

Here f and m are the internal shear force and the moments (very similar to structural mechanics.) The Hermite polynomials for interpolation are

where 

and the displacement field at any point is
	(Subscript[u, y]^h)(x) = (Subscript[N, u1])(Subscript[u, y1])+ (Subscript[N, θ1])(Subscript[θ, 1])+ (Subscript[N, u2])(Subscript[u, y2])+ (Subscript[N, θ2])(Subscript[θ, 2])
For this to be valid interpolation field, it should satisfy the Kroeneker delta property at the nodes, which it indeed satisfies

This ensures that the continuities at the nodes both for the displacement  (Subscript[u, y]) and the slope (Subscript[du, y])/dx are indeed satisfied. 


The Hermite functions for two node beams. 

Note that the derivative transformation between x and ξ is:

Also note that any integral of the form: (∫_(Ω^e)) f(x) dx = (l/2)(∫_(-1)^(+1))f(ξ) dξ.
When we use the Galerkin scheme, the same interpolation functions are used for (Subscript[u, y]) and w.

The second derivatives are expressed as

Now, we check the properties of the shape functions . For simplicity (without loosing generality) we choose the value of the length of the beam as  l = 1.

Discrete equations

As before we can the residual (or the reactions) as

The element stiffness matrix is:

The external force matrix is:

where (Subscript[f, Ω]^e) and (Subscript[f, Γ]^e)  are the element and body forces, respectively. Everything else now is similar to what we have been doing before. 

If the stiffness E I is constant over the element, the element stiffness is given as

For constant pressure

It can be seen that uniform load results in both nodal moments and vertical nodal forces, corresponding to fixed end moments for a uniformly loaded beam. 

It itself such beam elements are particularly different that what we do in structural mechanics. However, they provide some understanding of what one does when encountered with higher order differential equations. 

Also, as we see below, these ideas could be extended for finite element analysis of shells.

Example problem:

Consider a beam problem shown in Figure below. The beam ANS is clamped at the left side and is free at the right side. Spatial dimensions are in meters forces in N and distributed loading p in N (m^(-1)). The beam bending stiffness is EI = (10^4)N (m^2). The natural boundary conditions at x = 12 m are (Overscript[s, -]) = -20 N and (Overscript[m, -]) = 20 N m.

The beam is subdivided into two finite elements as shown in this figure. 

Step-1: Degrees of freedom and preliminary steps

The global displacement matrix is defined as
: 
Element stiffness matrix

For element-1: E I = (10^4) , L = 8:

and similarly for element-2 : E I = (10^4), L = 4

Global stiffness matrix:

The global stiffness matrix is computed using direct assembly:

Boundary force matrix:

For element 1:  (Subscript[f, Γ]^((1))) = (([0, 0, 0, 0 ])^T) because it does not have boundary (Subscript[Γ, s]) or (Subscript[Γ, m]). 
For element 2: 

quite obviously because of the properties of the shape functions and their derivatives. 

The global boundary force matrix is obtained by direct assembly
     
     Body force matrix:
     
     
  We also have a point force acting at the center of element-1. Since the nodal forces transmitted due to a body force of form p(x) = (Subscript[F, 0])δ(x - (Subscript[x, 0])) is:
  	(Subscript[f, Ω]^e) = (∫_(Ω^e))(N^eT)(Subscript[F, 0])δ(x - (Subscript[x, 0])) dx = (Subscript[F, 0])(N^eT)((Subscript[x, 0]))=(Subscript[F, 0])(N^eT)((Subscript[ξ, 0])), 
  where (Subscript[ξ, 0]) corresponds to the point (Subscript[x, 0]) through the equation:
  	ξ = ((2 x)/l)- 1.
  Note, that here we do not have the additional factor of (l/2) because δ(x-(Subscript[x, 0])) has units of (1/length) (units of (1/x)). However, ξ does not have any units. 
For element 1:  In this case  a distributed force of p(x) = -1 acts over the beam an a point force -10 N acts at  x = (l/2) which implies that ξ = 0.  So the total body force on this element is:
	
For element 2: The point force in element 2, acts on the first node where ξ = -1 giving

The direct assembly of the force matrix gives the global force vector as:

While accounting for the essential boundary conditions ((Subscript[Γ, u]) and (Subscript[Γ, θ])) as the partitioned degrees of freedom, the stiffness matrix for the system is:

As usual solving for the displacements and the unknown reactions at the essential boundary point we get

Post-processing: Bending moment and shear force



How does this answer (displacement, shear force and bending moment) compare with the exact answer? The displacement is quite reasonable as compared with the exact solution. Consider the origin of the shape functions for the beam.

u/.DSolve[{EI u''''[x]==0, u[0]==u1, u'[0] == 0,u'[l]==0,u[l]==0}, u, x][[1]];
%[x]//FullSimplify;
N1= %/.x->(1 + ξ) (l/2)//FullSimplify
u/.DSolve[{EI u''''[x]==0, u[0]==0, u'[0] == θ1,u'[l]==0,u[l]==0}, u, x][[1]];
%[x]//FullSimplify;
N2= %/.x->(1 + ξ) (l/2)//FullSimplify
u/.DSolve[{EI u''''[x]==0, u[0]==0, u'[0] == 0,u'[l]==0,u[l]==u2}, u, x][[1]];
%[x]//FullSimplify;
N3= %/.x->(1 + ξ) (l/2)//FullSimplify
u/.DSolve[{EI u''''[x]==0, u[0]==0, u'[0] == 0,u'[l]==θ2,u[l]==0}, u, x][[1]];
%[x]//FullSimplify;
N4= %/.x->(1 + ξ) (l/2)//FullSimplify

These shape functions are the same as what we have used. Thus shape functions that we have used earlier infact corresponding to solving the basic beam equation with boundary displacements at nodes 1 and 2 (Subscript[u, 1]), (Subscript[θ, 1]), (Subscript[u, 2]), (Subscript[θ, 2]).  As a result any problem with loading only at the nodes will be exactly solved using the FEA formulation discussed. However, bending moment will be approximation because with the FEA shape functions bending moments will be linear in x and will have kinks or jumps when we move from one element to the other if there are external force or bending moment, respectively, at the common node. The worst, however, will be shear force since wherever we have a point force on the beam, we will end up creating jump in the shear force. Hence, it is  a good idea to use a new node wherever point forces are present. The actual comparison between the current FEA solution with two elements and the exact solution is shown below.

What next?

The polynomial space of (C^1)elements that was used above of Hermite polynomials is not available in FEniCS. An alternative procedure which could be implemented would be to utilize (C^0) elements, but put a constraint on the kinks between the elements using the so called Discontinous Galerkin elements (https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html). 

The simple Euler-Bernoulli formulation can be extended to deep beams where even shear can be dominant. A most commonly used model for this is what is called as Timoshenko beam theory (https://en.wikipedia.org/wiki/Timoshenko_beam_theory). The corresponding finite element formulation can given rise to what is called as shear locking phenomena, whose effect can be softened using the so called reduced integration formulation (http://14.139.134.16/cmmacs/pdf/ch06.pdf).   

There are many, many other interesting features and abnormalities that can arise in finite element formulation which are topics for a more advance course. Further, the natural transition from beam elements is to shear elements. We will not go into such details in our course. Instead, I will provide these as some of the topics on which you will do a mini-course project and do a 5 minute presentation.