8.2 The Finite Element Method for dynamic linear elasticity

 

In this section, we show how to extend the finite element method to solve dynamic problems. Specifically:

 

1. The principle of virtual work is used to derive a discrete system of equations for the (time varying) nodal displacements

 

2. Three methods for integrating the equation of motion with respect to time are presented: (i) explicit time-stepping; (ii) implicit time stepping; and (iii) modal analysis.  The properties of each scheme are illustrated by solving simple problems.

 

3. As always example codes are provided so you can see how the method is actually implemented, and explore its predictions for yourself.

 

 

 

8.2.1 Review of the governing equations of dynamic linear elasticity

 

 

As before, we begin by summarizing the governing equations for a standard dynamic linear elasticity problem.   A generic linear elasticity problem is shown in the figure. We are given:

 

1. The shape of the solid in its unloaded condition $R_0\approx R$

 

2. The initial stress field in the solid (we will take this to be zero in setting up our FEM code)

 

3. The elastic constants for the solid $C_{ijkl}$

 

4. The thermal expansion coefficients for the solid, and temperature distribution (we will take this to be zero for our FEM code, for simplicity)

 

5. The initial displacement field in the solid $\stackrel{o}{u}(x)$, and the initial velocity field $\stackrel{o}{v}(x)$

 

6. A body force distribution $b(x,t)$ acting on the solid (Note that in this section we will use b to denote force per unit volume rather than force per unit mass, to avoid having to write out the mass density all the time)

 

7. Boundary conditions, specifying displacements $u^{*}(x,t)$ on a portion ${\partial }_{1}R$ or tractions $t^{*}(x,t)$ on a portion ${\partial }_{2}R$ of the boundary of R

 

We then wish to calculate (time varying) displacements, strains and stresses $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ satisfying the governing equations of dynamic linear elasticity

 

1. The strain-displacement equation ${\epsilon }_{ij}={\displaystyle \frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)}$

 

2. The elastic stress-strain law ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$

 

3. The equation of motion for stresses $\partial {\sigma }_{ij}/\partial x_i+b_j=\rho {\partial }^2{u}_j/\partial t^2$

 

4. The boundary conditions on displacement and stress

$u_i=u_i^{*}\text{on}{\partial }_{1}R{\sigma }_{ij}{n}_i=t_j^{*}\text{on }{\partial }_{2}R$

 

 

 

8.2.2 Expressing the governing equations using the principle of virtual work

 

Just as for static problems, the principle of virtual work can be used to write the governing equation for the displacement field in a linear elastic solid in an integral form (called the `weak form’).  Instead of solving the governing equations listed in the preceding section, the displacements, strains and stresses $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ are calculated as follows. 

 

1. Find a displacement field $u_i(x_j)$ satisfying

$\begin{array}{l}{\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial u_k}{\partial x_l}{\displaystyle \frac{\partial \delta v_i}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}dV+{\displaystyle \underset{R}{\int }\rho {\displaystyle \frac{{\partial }^2{u}_i}{\partial t^2}\delta v_{i}dV}}-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}\delta v_{i}dA=0}}}}}\\ u_i=u_i^{*}\text{on}{\partial }_{1}R\end{array}$

for all virtual velocity fields $\delta v_i$ satisfying $\delta v_i=0$ on ${\partial }_{1}R$.

 

2. Compute the strains from the definition ${\epsilon }_{ij}={\displaystyle \frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)}$

 

3. Compute the stresses from the stress-strain law ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$

 

 

The stress will automatically satisfy the equation of motion and boundary conditions, so all the field equations and boundary conditions will be satisfied.

 

You can derive this result by following the method outlined in Section 8.1.2.

 

 

 

8.2.3 The finite element equations of motion for linear elastic solids

 

We now introduce a finite element approximate following exactly the same procedure that we used for static problems. We choose to calculate the displacement field (now a function of time) at a set of n discrete nodes, as sketched in the figure. We denote the coordinates of these special points by $x_i^a$, where the superscript a ranges from 1 to n.  The unknown displacement vector at each nodal point will be denoted by $u_i^a$.  We then set up the FEA equations for dynamic problems as follows:

 

1. As before, the displacement and virtual velocity fields are interpolated between nodal values as

$u_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)u_i^a}$       $\delta v_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)\delta v_i^a}$

 

2. Substitute into the virtual work equation to see that

${\displaystyle \underset{R}{\int }\rho N^b{N}^a{\displaystyle \frac{{\partial }^2{u}_i^b}{\partial t^2}\delta v_i^{a}dV}}+{\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N_{}^b}{\partial x_l}{\displaystyle \frac{\partial N_{}^a}{\partial x_j}{u}_k^b\delta v_i^{a}dV}}}-{\displaystyle \underset{R}{\int }{b}_i{N}^a\delta v_i^{a}dV}-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}{N}^a\delta v_i^{a}dA}=0$

 

3. This is once again a linear system which may be expressed in the form

$\left(M_{ab}{\ddot{u}}_i^b+K_{aibk}{u}_k^b-F_i^a\right)\delta v_i^a=0$

where

$K_{aibk}={\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^a(x)}{\partial x_j}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}dV}}}{F}_i^a={\displaystyle \underset{R}{\int }{b}_i{N}^a(x)}dV+{\displaystyle \underset{\partial 2R}{\int }{t}_i^{*}{N}^a(x)}dA$

are the finite element stiffness matrix and force vector introduced in Section 7.2 (the matrix form of the stiffness in Section 8.1.13 is usually most convenient), and

$M_{ab}={\displaystyle \underset{R}{\int }\rho N^a{N}^{b}dV}$

is known as the finite element mass matrix.  Procedures for computing the mass matrix are given below.

 

4. Again, since the principal of virtual power holds for all virtual velocity fields such that $\delta v_i^a=0$ for $x_i^a$ on ${\partial }_{1}R$ we conclude that this requires

$\begin{array}{l}\left(M_{ab}{\displaystyle \frac{{\partial }^2{u}_i^b}{\partial t^2}+K_{aibk}{u}_k^b-F_i^a}\right)=0\forall (i,a):x_k^a\text{not on}{\partial }_{\text{1}}R\\ u_i^a=u_i^{*}\forall (i,a):x_k^a\text{ on}{\partial }_{\text{1}}R\end{array}$

Nodal velocities must also satisfy initial conditions. This is a set of coupled second-order linear differential equations that may be integrated with respect to time to determine $u_i^a$ as a function of time.  Observe that $M_{ab}$ and $K_{aibk}$ are constant matrices for linear elastic problems, but $F_i^a$ will in general be a function of time.

 

 

To find the solution, we now need to integrate these equations with respect to time.  For linear problems, there are two choices

 

1. Brute-force time stepping methods (e.g. Newmark time integration).  This technique can be used for both linear and nonlinear problems.

 

2. Modal methods, which integrate the equations of motion exactly.  This method only works for linear problems.

 

Both methods are described in detail below.

 

 

 

8.2.4 Newmark time integration for elastodynamics

 

The general idea is simple.  Given values of $u_i^a$ and $\partial u_i^a/\partial t$ at some time t, we wish to determine values at a slightly later time $t+\Delta t$, using the governing equations of motion.  There are of course many ways to do this.  Here, we will present the popular ‘Newmark’ time integration scheme.

 

Newmark method for a 1 DOF system: It is simplest to illustrate the procedure using a 1-DOF system (e.g. a forced,  undamped spring-mass system) first.   For this case the equation of motion is

 

$m\ddot{u}+ku-f=0$

 

where f(t) and the initial conditions $u(0),\dot{u}(0)$ are given.

 

Suppose that we are able to get estimates for the acceleration $\ddot{u}(t)$ and $\ddot{u}(t+\Delta t)$ both at the start and end of a general time step $\Delta t$.  We could then use a Taylor series expansion to obtain estimates of displacement and velocity at time $t+\Delta t$

$\begin{array}{l}u(t+\Delta t)\approx u(t)+\Delta t\dot{u}(t)+{\displaystyle \frac{\Delta t^2}{2}\left[(1-{\beta }_2)\ddot{u}(t)+{\beta }_2\ddot{u}(t+\Delta t)\right]}\\ \dot{u}(t+\Delta t)\approx \dot{u}(t)+\Delta t\left[(1-{\beta }_1)\ddot{u}(t)+{\beta }_1\ddot{u}(t+\Delta t)\right]\end{array}$

Here, ${\beta }_1$ and ${\beta }_2$ are two adjustable parameters that determine the nature of the time integration scheme.  If we set ${\beta }_1={\beta }_2=0$ the acceleration is estimated based on its value at time t.  This is known as an explicit time integration scheme.  Alternatively, if we put ${\beta }_1={\beta }_2=1$, the acceleration is estimated from its value at time $t+\Delta t$.  This is known as implicit time integration.

 

The acceleration at the start and end of the time step is computed using the equation of motion.  At time t we have

$\ddot{u}(t)={\displaystyle \frac{1}{m}\left[-ku(t)+f(t)\right]}$

while at time $t+\Delta t$

$\begin{array}{l}m\ddot{u}(t+\Delta t)+ku(t+\Delta t)-f(t+\Delta t)=0\\ \Rightarrow m\ddot{u}(t+\Delta t)+k\left\{u(t)+\Delta t\dot{u}(t)+{\displaystyle \frac{\Delta t^2}{2}\left[(1-{\beta }_2)\ddot{u}(t)+{\beta }_2\ddot{u}(t+\Delta t)\right]}\right\}-f(t+\Delta t)=0\\ \Rightarrow \ddot{u}(t+\Delta t)={\displaystyle \frac{1}{m+k{\beta }_2\Delta t^2/2}\left\{-k\left[u(t)+\Delta t\dot{u}(t)+{\displaystyle \frac{\Delta t^2}{2}(1-{\beta }_2)\ddot{u}(t)}\right]+f(t+\Delta t)\right\}}\end{array}$

 

This then gives us the following time stepping scheme: Given: k, m, f(t), $u(0),\dot{u}(0)$

1. At t=0 compute $\ddot{u}(0)={\displaystyle \frac{1}{m}\left[-ku(0)+f(0)\right]}$.

 

2. Then for successive time steps, compute

$\ddot{u}(t+\Delta t)={\displaystyle \frac{1}{m+k{\beta }_2\Delta t^2/2}\left\{-k\left[u(t)+\Delta t\dot{u}(t)+{\displaystyle \frac{\Delta t^2}{2}(1-{\beta }_2)\ddot{u}(t)}\right]+f(t+\Delta t)\right\}}$

$\begin{array}{l}u(t+\Delta t)\approx u(t)+\Delta t\dot{u}(t)+{\displaystyle \frac{\Delta t^2}{2}\left[(1-{\beta }_2)\ddot{u}(t)+{\beta }_2\ddot{u}(t+\Delta t)\right]}\\ \dot{u}(t+\Delta t)\approx \dot{u}(t)+\Delta t\left[(1-{\beta }_1)\ddot{u}(t)+{\beta }_1\ddot{u}(t+\Delta t)\right]\end{array}$

 

 

Newmark applied to FEM equations: This can immediately be extended to the general n DOF case as follows.

 

Given: $K_{aibk}$, $M_{ab}$, $F_i^a(t)$, $u_i^a(0),{\dot{u}}_i^a(0)$

Then:

 

1. At  t=0, compute ${\ddot{u}}_i^a(0)=M_{ab}^{-1}\left[-K_{bick}{u}_k^c(0)+F_i^b(0)\right]$

 

2. Then for successive time steps, solve

$\begin{array}{l}{M}_{ab}{\ddot{u}}_i^b(t+\Delta t)+{\displaystyle \frac{1}{2}{\beta }_2\Delta t^2{K}_{aibk}{\ddot{u}}_k^b(t+\Delta t)=}\\ -K_{aibk}\left[u_k^b(t)+\Delta t{\dot{u}}_k^b(t)+{\displaystyle \frac{\Delta t^2}{2}(1-{\beta }_2){\ddot{u}}_k^b(t)}\right]+F_i^a(t+\Delta t)\end{array}$

for ${\ddot{u}}_i^b(t+\Delta t)$

 

3. Then compute

$\begin{array}{l}{u}_i^a(t+\Delta t)\approx u_i^a(t)+\Delta t{\dot{u}}_i^a(t)+{\displaystyle \frac{\Delta t^2}{2}\left[(1-{\beta }_2){\ddot{u}}_i^a(t)+{\beta }_2{\ddot{u}}_i^a(t+\Delta t)\right]}\\ {\dot{u}}_i^a(t+\Delta t)\approx {\dot{u}}_i^a(t)+\Delta t\left[(1-{\beta }_1){\ddot{u}}_i^a(t)+{\beta }_1{\ddot{u}}_i^a(t+\Delta t)\right]\end{array}$

 

 

 

 

8.2.5 1-D implementation of a Newmark scheme.

 

It is straightforward to extend a static FEM code to dynamics.  As before, we illustrate the method in 1D before developing a 3D code.

 

Consider a long linear elastic bar as shown in the figure.  Assume

 

1. The bar has shear modulus $\mu$ and Poisson’s ratio $\nu$  and mass density $\rho$;

 

2. The bar has cross section $h\times h$ and length L;

 

3. The bar is at rest and strain free at t=0 with $u_1=d{u}_1/dt=0$ 

 

4. It is constrained on all its sides so that $u_2=u_3=0$

 

5. The bar is subjected to body force $b=b(x_1,t)e_1$,

 

6. The bar is either loaded or constrained at its ends, so that the boundary conditions are either  $t_1=t^{*}(0,t),t_1=t^{*}(L,t)$ or displacement $u_1(0)=u^{*}(0)u_1(L)=u^{*}(L)$ at x=0 and x=L. 

 

 

For this 1-D example, then, the finite element equations reduce to

$M_{ab}{u}_1^b+K_{ab}{u}_1^b=F_{}^a$

where

$\begin{array}{l}{M}_{ab}=h^2{\displaystyle \underset{0}{\overset{L}{\int }}\rho N^a(x_1)N^b(x_1)d{x}_1}{K}_{ab}=h^2{\displaystyle \underset{0}{\overset{L}{\int }}{\displaystyle \frac{2\mu (1-\nu )}{1-2\nu }{\displaystyle \frac{\partial N^a(x_1)}{\partial x_1}{\displaystyle \frac{\partial N^b(x_1)}{\partial x_1}d{x}_1}}}}\\ F_{}^a=h^2{\displaystyle \underset{0}{\overset{L}{\int }}{b}_{}{N}^a(x_1)d{x}_1}+h^2{t}_1^{*}(0)N^a(0)+h^2{t}_{}^{*}(L)N^a(L)\end{array}$

 

We adopt exactly the same interpolation scheme as for the static solution (a 1D finite element mesh is shown in the figure)  and calculate stiffness matrix and force vector accordingly.  In addition, we need to determine the mass matrix.  This can also be done by evaluating integrals for each element by numerical quadrature, and then assembling element mass matrices into a global matrix.  Specifically,

 

· For each element, compute the element mass matrix. The mass matrix can either be evaluated by numerical quadrature, by evaluating the formula

$m_{ab}=h^2{\displaystyle \sum _{I=1}^{n_I}{w}_I\rho N^a\left({\xi }_I\right)N^b\left({\xi }_I\right)J({\xi }_I)}$

 

where

$J=\left|{\displaystyle \frac{\partial x}{\partial \xi }}\right|=\left|{\displaystyle \sum _{a=1}^{N_e}{\displaystyle \frac{\partial N^a({\xi }_I)}{\partial \xi }{x}_{}^a}}\right|$

and the integration points ${\xi }_I,w_{I}I=\mathrm{1...}{n}_I$ , and shape functions $N^{(a)}(\xi ),a=\mathrm{1...}{N}_e$ were listed in Section 7.2.5. Note that that a higher order integration scheme is required to integrate the mass matrix than we needed for the stiffness matrix.  Specifically, for quadratic elements we need a cubic  integration scheme (with 3 integration points) while for linear elements we need a quadratic scheme (2 integration points).  If the mass matrix is under-integrated it will be singular (you can check this using the example code provided below).  Alternatively, the mass matrix can be integrated analytically, which yields the results in Table 8.10)

 

· Assemble the contribution from each element to the global mass matrix $M_{ab}={\displaystyle \sum _{l=1}^L{m}_{ab}}$

 

· Modify the stiffness and mass matrices  to enforce the constraints ${\ddot{u}}^1={\ddot{u}}^{*}(0,t){\ddot{u}}^{(L)}={\ddot{u}}^{*}(L,t)$

 

· Run the Newmark time-stepping scheme outlined above.

 

 

 

8.2.6 Example 1D dynamic FEM code and solution

 

A simple example Matlab code for this 1D example is given in the file FEM1D_newmark.m.  It is set up to solve for displacements for a bar with the following properties:

 

1. Length 5 and unit unit x-sect area,

 

2. Mass density 100; shear modulus 50, Poisson’s ratio 0., (so the wave speed is 1);

 

3.  No body force,  displacement u=0 at x=0, and traction t=10 applied at x=L suddenly at t=0 and then held constant for t>0

 

The code can be downloaded from

https://github.com/albower/Applied_Mechanics_of_Solids

 

As an example, the figure on the right shows the variation of displacement at the end of the bar (x=L) with time  You should be able to verify for yourself that the exact solution is a saw-tooth, amplitude 0.5 and period 20 (the time for a plane wave to make two trips from one end of the bar to the other)

 

Notice how the peaks of the saw-tooth are blunted $–$ this is because the FEM interpolation functions cannot resolve the velocity discontinuity associated with a propagating plane wave.

 

The Matlab file produces animations of the displacement and velocity field in the bar.  In the exact solution, a plane wave propagates down the bar, and repeatedly reflects off the free and fixed ends of the bar.  You can see this quite clearly in the simulation results, but the wave-front is not sharp, as it is supposed to be.  This gives rise to spurious high frequency oscillations in the velocity fields.                   

 

It is interesting to investigate the effects of the two adjustable parameters in the time integration scheme. First, we look at the solution with ${\beta }_1=0.5,{\beta }_2=0$ (displacement update is fully explicit; the velocity update is a mid-point rule $–$ this scheme is used frequently in explicit dynamic simulations, because each time-step can be computed very quickly).  The figures below show the predicted displacement at the end of the bar for two values of time step $\Delta t$.

 

        

 

The result with smaller time step is good, but with a larger time step the solution is oscillatory.  This is an example of a numerical instability, which is a common problem in explicit time stepping schemes.  Eventually the solution blows up completely, because the oscillations grow exponentially.  With a small enough time step an explicit scheme will usually work, but the critical time step for stability is proportional to the time required for a wave to propagate through the smallest element in the mesh.

 

With larger values of ${\beta }_2$ the problem disappears.  In fact, one can show that for the equation of motion considered here, setting

${\beta }_2\ge {\beta }_1{\beta }_1\ge 1/2$

makes the time stepping scheme is unconditionally stable $–$ no oscillations will occur even for very large time steps.

 

Stability does not necessarily mean accuracy, however.  Results with a fully implicit integration scheme ( ${\beta }_1={\beta }_2=1$ ) and a large time step are shown below.  This shows a different problem $–$ energy is dissipated due to the numerical time integration scheme.   So, larger values of $\beta$ buy  stability by introducing artificial damping, at the expense of a loss of accuracy.  A good compromise is to set ${\beta }_1={\beta }_2=1/2$.

 

 

 

8.2.7 Lumped mass matrices

 

In the Newmark time integration scheme, accelerations are computed by solving a set of linear equations

$\begin{array}{l}{M}_{ab}{\ddot{u}}_i^b(t+\Delta t)+{\displaystyle \frac{1}{2}{\beta }_2\Delta t^2{K}_{aibk}{\ddot{u}}_k^b(t+\Delta t)=}\\ -K_{aibk}\left[u_k^b(t)+\Delta t{\dot{u}}_k^b(t)+{\displaystyle \frac{\Delta t^2}{2}(1-{\beta }_2){\ddot{u}}_k^b(t)}\right]+F_i^a(t+\Delta t)\end{array}$

This is the most time-consuming part of the procedure.  But notice that if we set ${\beta }_2=0$ and somehow find a way to make the mass matrix diagonal, then equation solution is trivial.

 

The so-called `lumped’ mass matrix is a way to achieve this.  Instead of using the correct element mass matrix

$M_{ab}={\displaystyle \underset{R}{\int }\rho N^a{N}^{b}dV}$

we replace it with a diagonal approximation.  Various procedures can be used to do this. 

 

1. Integrate the mass matrix using integration points located at the element nodes

 

2. Diagonal scaling procedure: Set $M_{aa}=c{M}_{aa}$ and $M_{ab}=0a\ne b$, with the constant c selected to satisfy

${\displaystyle \sum _a^{}{M}_{aa}}={\displaystyle \underset{V_e}{\int }\rho dV}$

3. Row-sum method: Set

$M_{aa}={\displaystyle \sum _b^{}{M}_{ab}}$ and $M_{ab}=0a\ne b$.

 

In general, the three approaches give different answers.  For some element types, it turns out that some of these procedures give zero or negative masses, which is clearly undesirable.   Use whichever method gives the best answers.

 

Explicit time integration with a lumped mass matrix has been implemented in the 1D FEM code FEM1D_newmark.m. To use a lumped mass matrix, set the lumpedmass variable to true near the top of the code.

 

 

 

As an example, the figure above  shows the predicted displacement at the end of the bar with relatively large timestep (set dt=0.1, beta1=1 and lumpedmass=true in the code to reproduce the result)

 

In general lumped mass matrices introduce some additional damping which can help stabilize the time stepping scheme, but lead to some accuracy loss.  You can see this in the result shown in Figure 8.20.  Using a smaller timestep reduces the energy loss.

 

Explicit time integration is cheap, easy to implement and is therefore a very popular technique.  Its disadvantages are that it is conditionally stable and can require very small timesteps (the critical step size for stability is proportional to the size of the smallest element in the mesh and inversely proportional to the wave speed). 

 

 

 

8.2.8 Example 2D and 3D dynamic linear elastic code and solution

 

As always, a sample code can be downloaded from

https://github.com/albower/Applied_Mechanics_of_Solids

 

1. The code is in a file called FEM_2Dor3D_linelast_dynamic.m.  You will need to modify the line of the code that reads the input file to point to wherever you store the input file, as described in 8.1.14.

 

2. An input file that sets up a simulation of the vibration of a cantilever beam that is subjected to a transverse force at its end is provided in the file Linear_elastic_dynamic_beam.txt.

 

The format for the input file is described in detail in Section 8.1.14.  The mesh is assumed to be at rest at time t=0, and any loads are applied as a step and held constant for t>0.  You can always extend the code to do more complex loadings if you need to.  The motion of the beam is animated at the end of the file.

 

 

 

8.2.9 Modal method of time integration

 

The finite element equations

$\left(M_{ab}{\ddot{u}}_i^b+K_{aibk}{u}_k^b-F_i^a\right)=0$

are a standard set of coupled, second order linear differential equations such as one would meet in analyzing forced vibrations in discrete systems.  They can therefore be solved using the usual modal techniques.  The procedure is to rearrange the system of equations to yield n decoupled 2^nd order differential equations, by means of the following substitutions.  We will drop index notation in favor of matrix notation.

 

1. Write the governing equation as $M\ddot{u}+Ku=F(t)$

 

2. Let

$q=M^{1/2}u$              $H=M^{-1/2}K{M}^{-1/2}$

and substitute for u

$\ddot{q}+Hq=M^{-1/2}F$

 

3. Note that H is a symmetric, positive definite matrix.  Consequently, we may perform a spectral decomposition and write

$H=Q\Lambda Q^T$

where Q is an orthogonal matrix ( $Q^{T}Q=I$ ) and $\Lambda$ is diagonal.  This spectral decomposition is accomplished as follows: compute the eigenvalues ${\lambda }_i$ and corresponding normalized eigenvectors $r^{(i)}$ of H (normalized to satisfy $r^{(i)}\cdot r^{(i)}=1$ ), then set

$\Lambda =\left[\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right]Q=\left[\begin{array}{cccc}{r}_1^{(1)}& r_1^{(2)}& & \\ r_2^{(1)}& r_2^{(2)}& & \\ ⋮& & \ddots & \\ r_n^{(1)}& & & r_n^{(n)}\end{array}\right]$

 

4. Next let $w=Q^{T}q$ and rearrange the equation of motion a second time to get

$\ddot{w}+\Lambda w=Q^T{M}^{-1/2}F$

Since $\Lambda$ is diagonal this is now a set of n uncoupled ODEs for w.

 

5. The initial conditions for w follow as

$w=Q^T{M}^{1/2}u(t=0)\dot{w}=Q^T{M}^{1/2}\dot{u}(t=0)$

 

6. Solve the decoupled ODEs for w, then  recover the displacements through

$u=M^{-1/2}Qw$

You can use any method you like to solve the ODEs for w $–$ for harmonic forcing you could use the exact solution.  For general F(t) you could use Fourier transforms (the FFT works for periodic F).  You could even use the Newmark algorithm if you want, although this would be rather perverse!

 

 

Effectively, this procedure reduces the general forced FEM problem to an equivalent one that involves solving the equations of motion for n forced, spring-mass systems.

 

 

 

8.2.10 Natural frequencies and mode shapes

 

The procedure outlined in 8.2.9 can also be used to extract the natural frequencies and mode shapes for a vibrating, linear elastic solid.  To see how to do this, note that

 

1. By definition, the natural frequencies and mode shapes of a continuous solid or structure are special deflections and frequencies for which the solid will vibrate harmonically.

 

2. The governing equations for w can be regarded as describing the vibration of n uncoupled spring-mass systems;

 

3. If we solve the problem with F=0, and excite only the ith spring-mass system (e.g. by choosing appropriate initial conditions), the solution will be a harmonic vibration at the natural frequency of the ith spring-mass system. 

 

 

4. The natural frequencies of the n decoupled systems are therefore the natural frequencies of vibration of the structure.  The corresponding eigenvectors $r^{(i)}$ are related to the deflections $u^{(i)}$ associated with the ith vibration mode through $u^{(i)}=M^{-1/2}{r}^{(i)}$.

 

 

 

The following procedure can be used to determine the natural frequencies and mode shapes of a linear elastic solid:

 

1. Compute the finite element mass and stiffness matrices M and K

 

2. Find $H=M^{-1/2}K{M}^{-1/2}$

 

3. Compute the eigenvalues ${\lambda }_i$ and corresponding eigenvectors $r^{(i)}$ of H

 

4. The natural frequencies ${\omega }_i$, are related to the eigenvalues of H by ${\lambda }_i={\omega }_i^2$ .

 

5. The deflections $u^{(i)}$ associated with the ith vibration mode follow as $u^{(i)}=M^{-1/2}{r}^{(i)}$. The deflections calculated in this way will have some random magnitude $–$ if you like, you can normalize them appropriately.

 

The procedure is simple on paper (as always) but there are two problems in practice.  First, we need to find the square root of the mass matrix.  For a general matrix, we would have to compute a spectral decomposition of M to do this.  However, if we use a lumped mass matrix, its square root is easily found $–$ just take the square root of all the diagonal elements!  Secondly, computing all the eigenvalues and eigenvectors of a general dynamical matrix H is an exceedingly time consuming process.  To make this a viable scheme, we normally extract only the lowest few eigenvalues and the corresponding eigenvectors (say 10 to 100), and discard the rest. 

 

 

 

8.2.11 Example 1D code with modal dynamics

 

As an example, a version of the 1D code with modal dynamics is provided in the MATLAB code FEM_1D_modal.m.  The codes can be downloaded from

https://github.com/albower/Applied_Mechanics_of_Solids

 

 

The code is set up to plot the displacements associated with the fundamental (lowest frequency) vibration mode, as well as to calculate the displacement of the end of the bar as a function of time.  The code will handle computations with a full mass matrix as well as a lumped matrix, so you can compare the results with the two cases. 

 

 

 

 

The figure above  shows the mode shape (displacements) for the fundamental vibration mode, and the predicted motion at the end of the bar.  The code will also animate the displacement in the bar (not shown here).

        

Here is a brief comparison of the Newmark and modal time integration schemes:

 

1. For linear problems the modal decomposition method usually beats direct time integration in both accuracy and speed.

 

2. If you are interested in looking in detail at transient wave propagation through the solid there is not much difference between the two methods $–$ in this case you need to retain a huge number of terms in the modal expansion for accurate results.

 

3. For most vibration problems, however, it is generally only necessary to retain a small number of modes in the expansion, in which case the modal decomposition method will be vastly preferable to direct time integration.  Moreover, knowledge of the vibration modes can itself be valuable information. You can use also use the modal approach to handle very complex forcing, including random vibrations, directly.

 

4. The main limitation of the modal decomposition approach is that it can only handle undamped, linear systems (one can introduce modal damping for vibration applications, but this does not model any real energy dissipation mechanism).  This is of course far too restrictive for any real application except the simplest possible problems in vibration analysis.

 

 

 

8.2.12 Example 2D and 3D FEM code to compute mode shapes and natural frequencies

 

As always, a sample code can be downloaded from

https://github.com/albower/Applied_Mechanics_of_Solids

 

1. An example program is provided in the file FEM_2Dor3D_modeshapes.m (MATLAB);

 

2. The code can be run with the file Linear_elastic_dynamic_beam.txt $–$ it computes and plots mode shapes for the beam as shown below (you need to edit the code to select which mode is displayed)