Pixabay

Overview

This section will look at two additional ways of representing a real world system with a mathematical model.

• Shape functions can be combined, or superimposed, to create the complex shape of a dynamic system.
• A large structure can be represented as a number of small elements which can each individually be modelled as a simple component.

Topic 1 –  Shape Functions

In Section 1 of this module, during our discussion of the relationship between time domain and frequency domain responses, we demonstrated how a series of simple sinusoidal functions could be combined to create a complex waveform.

The logic behind the use of shape functions is similar.  Sometimes complicated parts take on complicated shapes……for example high speed rotors.

We need to find a way of representing these complex shapes.

Let us assume that we can represent the complex shape by some combination of simpler shapes (similar to the simple sinusoids that make up the complex wave form in our FFT discussion from Section 1).  In this case, let us try using the shapes of the sine wave, cosine wave, hyperbolic sine, and hyperbolic cosine.

Assume that the deflected shape, y(x), can be represented by a combination of these four functions, each with a different scale factor (A₁, A₂, A₃, A₄)

$\text{y(x)}={\text{A}}_1\text{cosλx}+{\text{A}}_2\text{sinλx}+{\text{A}}_3\text{coshλx}+{\text{A}}_4\text{sinhλx}$ 

Let us assume that our component is cantilevered at one end, so it has boundary conditions at x = 0 of deflection  y(0) = 0    and slope   y’(0) = 0    and it is simply supported at the other end, so it has boundary conditions at x = L of deflection y(L) = 0  and  y’’(L) = 0

Take the derivatives of our assumed shape functions

$\text{y(x)}={\text{A}}_1\text{cosλx}+{\text{A}}_2\text{sinλx}+{\text{A}}_3\text{coshλx}+{\text{A}}_4\text{sinhλx}$ 

$\text{y’(x)}=-{\text{A}}_1\text{λsinλx}+{\text{A}}_2\text{λcosλx}+{\text{A}}_3\text{λsinhλx}+{\text{A}}_4\text{λcoshλx}$ 

$\text{y’’(x)}=-{\text{A}}_1{\text{λ}}^2\text{cosλx}-{\text{A}}_2{\text{λ}}^2\text{sinλx}+{\text{A}}_3{\text{λ}}^2\text{coshλx}+{\text{A}}_4{\text{λ}}^2\text{sinhλx}$ 

Apply the boundary condition   y(0) = 0  

$\text{y(x)}=0={\text{A}}_1\text{cos0}+{\text{A}}_2\text{sin0}+{\text{A}}_3\text{cosh0}+{\text{A}}_4\text{sinh0}$  

$0={\text{A}}_1+{\text{A}}_3$ 

${\text{A}}_3=-{\text{A}}_1$ 

Apply the boundary condition  y’(0) = 0     (Recall cosh(0) = 1 and sinh(0) = 0)

$\text{y’(x)}=0=-{\text{A}}_1\text{λsin0}+{\text{A}}_2\text{λcos0}+{\text{A}}_3\text{λsinh0}+{\text{A}}_4\text{λcosh0}$ 

$0={\text{A}}_2+{\text{A}}_4$ 

${\text{A}}_4=-{\text{A}}_2$ 

Apply the boundary conditions     y(L) = 0        and     y’’(L) = 0

$\text{y(x)}=0={\text{A}}_1\text{cosλL}+{\text{A}}_2\text{sinλL}+{\text{A}}_3\text{coshλL}+{\text{A}}_4\text{sinhλL}$ 

$\text{y’’(x)}=0={\text{A}}_1{\text{λ}}^2\text{cosλL}-{\text{A}}_2{\text{λ}}^2\text{sinλL}+{\text{A}}_3{\text{λ}}^2\text{coshλL}+{\text{A}}_4{\text{λ}}^2\text{sinhλL}$ 

Use the previous results of   A₃ = -A₁      and    A₄ = -A₂    and get

$\text{(cosλL – coshλL)}{\text{A}}_1+\text{(sinλL – sinhλL)}{\text{A}}_2=0$ 

And

$-\text{(cosλL + coshλL)}{\text{A}}_1-\text{(sinλL + sinhλL)}{\text{A}}_2=0$ 

Put into matrix form

$\left(\begin{array}{cc}\text{(cosλL − coshλL)}& \text{(sinλL − sinhλL)}\\ \text{−(cosλL + coshλL)}& \text{−(sinλL + sinhλL)}\end{array}\right)\left(\begin{array}{c}{A}_1\\ A_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$ 

A₁ = 0 and A₂ = 0 would be a solution of this system, but it would be a trivial solution.  We could find a non-trivial solution from the determinant expression:

$\text{det}\left(\begin{array}{cc}\text{(cosλL − coshλL)}& \text{(sinλL − sinhλL)}\\ \text{−(cosλL + coshλL)}& \text{−(sinλL + sinhλL)}\end{array}\right)=0$ 

Expanding the determinant and using the trig identity cos²θ + sin²θ = 1

Yields:   $\text{cosλLsinhλL – sinλLcoshλL = 0}$ 

We cannot solve this explicitly.  However, MATLAB can solve it numerically, finding any and all values for λL that make the equation work.   Call these (λL)₁, (λL)₂, (λL)₃, etc.

We can now take each of these values (λL)_n  back to the earlier equation and solve for A₁ and A₂.  This would give us an expression for y(x) for each value of (λL)_n.   Think of these as mode shapes that the shaft might take.  Considering back to eigenvalues and eigenvectors, if these are mode shapes, then there should be a frequency associated with each mode shape.

Recognize that the (λL)_n  relate to the natural frequencies of the system.  However, since no actual geometric or material characteristics of the original beam were included in the problem, our result is only a ratio.

We must therefore dimensionalize our answer.  Recall ${\text{ω}}_{\text{n}}={\text{(k/m)}}^{\text{1/2}}$ 

Use the beam bending formula, $\text{δ}={\text{FL}}^3/\text{EI}$  and rearrange to solve for k, which is the spring rate  $\text{k}=\text{F}/\text{δ}=\text{EI}/{\text{L}}^3$ 

And the mass is the density times the volume, so for length L and Area A,  $\text{m = ρLA}$

So  ${\text{(k/m)}}^{\text{1/2}}={[\text{EI}/\text{ρ}{\text{AL}}^4]}^{\text{1/2}}$  and   ${\text{ω}}_2={\left({\text{λ}}_{\text{n}}\text{L}\right)}^2{\left[\text{EI/ρA}{\text{L}}^4\right]}^{\text{1/2}}$ 

Redemensionalizing our answer gives us frequencies of ${\text{ω}}_{\text{n}}={\left({\text{λ}}_{\text{n}}\text{L}\right)}^{\mathrm{2\; }}{\left[\text{EI/ρA}{\text{L}}^4\right]}^{\text{1/2}}$ in rad/sec.

Topic 2 – 1-D Finite Element Relationships

Finite Element Models (FEM) can be used to break extremely large structures down into a potentially very large number of simple elements.   The figure below demonstrates how this can result in a huge number of elements, as each intersection point on the web (called a node) is a single element and the lines of the web between the nodes represent the constraint relations between nodes.

Let us look at a much simpler model so that it will be much easier to follow how the computer builds these relationships.

Consider a small beam element, trapped between two rigid constraints.  Let P be an applied load and u₁ and u₂ be variables representing motion of the ends of the element.

Consider an element of length L with a fixed rigid constraint at the right end, being pushed with a force F₁ at the left end, so that there is an equal and opposite reaction F₁ at the other end, with an axial deflection u₁.

Recall the deflection of a beam in axial loading is δ = FL/AE

So for our situation above, F₁ = u₁AE/L

Now consider an element of length L with a fixed rigid constraint at the left end, being pulled with a force F₂ at the right end, so that there is an equal and opposite reaction F₂ at the other end, with an axial deflection u₂.

Similarly,   F₂ = u₂AE/L

The member is in static equilibrium so considering Figure 1 again, we could write the following two equations of equilibrium.

Put this in matrix form:

$\left(\begin{array}{c}\text{AE/L}\end{array}\right)\left(\begin{array}{cc}1& \text{−1}\\ \text{−1}& 1\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}P\\ P\end{array}\right)$ 

Now consider the beam from Figure 1 as if it were three elements, A, B, and C

Realise that the Load P is the same in each of the elements.  We will consider each element separately and apply half the load at each end.  Deflections at each end are also indicated.

This would allow us to create the three matrix equations shown below the figures.

Recognise that u_A2 = u_B1, let us call this simply u₂.

And u_B2 = u_C1, let us call this simply u₃.

For consistency, we will call u_A1 simply u₁   and  u_C3 simply u_4.

The three matrix equations would now be as follows:

Let us combine these three matrix equations into one.

Consider the boundary conditions.  At the left wall, the applied force will be reacted directly by the wall, and since it is a fixed boundary, u₁ will be zero.  Similarly, at the right wall, the force will be reacted directly by the wall and u₄ will be zero.  If we apply this boundary logic to the previous matrix equation:

This would represent the 3 element FEM that we created.  For a given P, A, E and L, we could solve it to determine the deflection at the element boundaries u₂ and u₃.

Using similar logic, if we were to divide the beam in Figure 1 into 9 elements, then we could create a similar matrix equation with a 10 by 10 matrix, from which we could predict the deflections u₁ through u₁₀ representing the deflection at the boundary between each of the elements.

Topic 2 – 2-D & 3-D Finite Element Relationships

Extending what was previously discussed for a 1D element, now consider a 2D element with four corner nodes and horizontal deflections (u) and vertical deflections (v) at each.

If we ponder the 1D example from the previous discussion, and consider that now we have 4 nodes, and 4 horizontal degrees of freedom (DoF) and 4 vertical DoF for a total of 8, our equation would look like:

Now imagine what a 3D cubic element with 8 nodes, and 24 DoF (3 at each node, one in the x direction, one in the y direction, and one in the z direction).  The square matrix would be 24 by 24, the DoF matrix would be a 24 character column, and there would be 24 possible applied forces. 

But we could use such an element to build a FEM of any structure.

There are many different types of finite elements, allowing the engineer performing an analysis of a structure to choose the most appropriate.   However, this discussion should be adequate to give you an understanding of the fundamentals. 

There are also a number of generally available FEA software routines available in the engineering world.   The one used by UHI is ANSYS, and the following are some videos showing how to build some relatively simple FEM models and performing stress and deflection FEA.

Following are links to videos demonstrating how to model a beam which is built-in at both ends (sometimes called a double-cantilever beam or an encastre beam)

Additional links to some useful videos will be included in the VLE.  However, there are hundreds available on the internet which have been made available for all types of modelling.

Summary

The advance in computing power has provided engineers with the ability to solve much larger and more complex problems.   Mathematical representations are required in order to do this.  This section has looked at the theory behind how this is done.  The existence of advanced modelling FEM and FEA routines makes it unnecessary for today’s engineer to build these mathematical representations themselves.   But it is fully appropriate that engineers still understand the concepts behind what the computer does on their behalf.

References and Bibliography

Ogata, K. (2004).  System Dynamics, 4th edition. Upper Saddle River, USA: Pearson Prentice Hall.