Modal Analysis of Plane Truss using Python

Bhairav Thakkar

5.00/5 (6 votes)

Sep 3, 2018

CPOL

8 min read

19953

508

Basic modal analysis for dynamic response of structures in Python

Download source code - 9.8 KB

Introduction

Structures vibrate under dynamic loads. These vibrations are of vital significance to the analyst and designer as dynamic loads often induce much higher structural response than static loads. Static analysis is comparatively simpler and solutions are available. Dynamic analysis requires a different set of linear algebraic operations. To obtain dynamic loads on a structure, modal analysis is required to be carried out. This article describes the steps to be carried out for performing modal analysis on structures by taking a 2D truss problem and showing Python code alongside.

The Python code has been structured for ease of understanding and allows modifying the code for more for implementing additional features. The complete Python code has been attached here for ready reference and to help follow the article better.

Background

While the article intends to be complete, the reader is expected to be equipped with fundamental understanding of structural mechanics. The reader is also urged to read through some fundamental texts on structural mechanics and structural dynamics. A quick review of the article on finite element analysis: https://www.codeproject.com/Articles/1207954/Finite-Element-Analysis-in-VB-Net would certainly help. Nonetheless, the basic equations that are solved for a dynamic system are presented below.

For a solid body with distributed mass as shown below:

Solid Body With Distributed Mass

The kinetic energy is given by:

$ T = \frac{1}{2} \int_V{u^T u \rho DVD} $

where elements of vector u are the deformations and are given by u = Nq

For structural dynamics problems, the q vector is time dependent and N are the spatial shape functions as obtained for the structural model. The kinetic energy T for an element is then obtained as:

$ T_e = \frac{1}{2} \dot{q}^T \left[\int_e{\rho N^T N dV }\right]\dot{q} $

The expression in square bracket is the element mass matrix:

$ m_e = \int_e{\rho N^T N dV} $

Since this mass matrix is consistent with the element shape functions, it is known as consistent mass matrix. Summing up the kinetic energy over the whole domain (summing up for each element) gives:

$ T = \sum_e{T_e} = \sum_e {\frac{1}{2} \dot{q}^T m_e \dot{q} } = \frac{1}{2}\dot{Q}^T M \dot{Q} $

The potential energy is given by:

$ \Pi = \frac{1}{2} Q^T K Q - Q^T F $

Using the Lagrangean L;

$ L = T - \Pi $

This eventually leads to the equation of motion:

$ M \ddot{Q} + K Q = F $

For free vibrations, force F = 0. Thus,

$ M \ddot{Q} + K Q = 0 $

This eventually leads us to the generalized eigenvalue problem

$ K U = \lambda M U $

Where, $U$ is the eigenvector representing the normalized deformations corresponding to the vibrating mode and $\lambda$ is corresponding eigenvalue which is the square of circular frequency $\omega$. The modal analysis involves solution of the eigenvalue problem stated here. The solution is not direct and requires iterative procedures.

Using the Code

The accompanying code has been written in Visual Studio 2017, but should be able to run on just about any Python 3.x interpreter. The entire Visual Studio solution has been attached here. The code requires numpy and matplotlib libraries. Efforts have been made to adhere to the Single Responsibility Principle. Separate classes have been made for each type of task. The most important classes are explained here. Before you start, you would need to add numpy and matplotlib environments.

The main file is PlaneTruss.py. The file imports various modules and libraries required. The part of the PlaneTruss.py file which calls each routine to perform the eigenvalue solution is shown below: You can see the complete file from the attachment.

# Read the input
inp = Input.Input()
inp.Read()

nnodes= inp.nnodes
nodes = inp.nodes

nele = inp.nele
elements = inp.elements
nbcnodes = inp.nbcnodes
supports = inp.supports
loads = inp.loads
filename = inp.filename

# Solve the truss

ndof = nnodes * 2
assembly = Assembly.Assembly()
assembly. GenerateKgAndMg(ndof,elements)

# generate load vector
assembly.GenerateLoadVector(loads)

# apply boundary conditions
assembly.ApplyBoundaryConditions(supports)

# Retrieve global stiffness matrix, mass matrix and load vector from the assembly object
kg = assembly.kg
mg = assembly.mg
r = assembly.r

# perform linear solution to get deformations for applied load
print("Performing linear solution... ", end="")
d = np.linalg.solve(kg,r)
print ("done.")

# lets reshape the displacements so that u and v components can be seen separately
disp = d.reshape(nnodes,2)

# write the output to the outputfile
output = Output.Output()
output.WriteOutputFile(inp.outfilename, nodes, elements, disp)

# write the output in matlab format so that we can plot it
mfile = filename + ".m"
dispscale = 400
output.WriteModelOnMatlabFile(mfile,nodes,elements,disp)

# compute eigenvalues and eigen vectors
e = Eigen.Eigen()
print("Performing eigenvalue solution... ", end="")
evec = e.Solve(kg,mg,1e-6)
#evec = e.RescaleEigenVectors(evec)

eval = e.eigenvalues
neval = len(eval)
print ("done.")

The steps performed here are as follows:

Read the input.
Generate elemental stiffness and mass matrices.
Assemble the elemental stiffness and mass matrices to form global stiffness and global mass matrices.
Generate the load vector.
Apply boundary conditions. Please note that only the stiffness matrix is modified in this step.
Retrieve the global matrices.
Perform static analysis. This is primarily to confirm that the structure has been correctly modeled and the boundary conditions have been applied appropriately.
Carry out eigen value analysis to get eigen values and eigen vectors which give the mode shapes.

Each entity is managed by its own object. For example, the Element class is responsible for managing state variables required to completely define an element. The present example considers 2D Truss element and hence stiffness and mass matrices for 2D truss are developed as:

$\begin{aligned} K_e = \frac{AE}{l} \begin{bmatrix} c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2 \end{bmatrix} \end{aligned} $

Where, A is the area of cross section, E denotes the modulus of elasticity of truss member material and l is the length of element / member. c and s indicate the direction cosine and sine of the member.

The consistent mass matrix is given by:

$\begin{aligned} M_e = \frac{1}{6} \rho A l \begin{bmatrix} 2 & 0 & 1 & 0 \\ 0 & 2 & 0 & 1 \\ 1 & 0 & 2 & 0 \\ 0 & 1 & 0 & 2 \end{bmatrix} \end{aligned} $

The above shown matrices are generated by the Element class shown below:

import math
import numpy as np
class Element(object):
    """2D Truss element. Contains various definitions required for the element and also computes
    element matrices"""

    def __init__(self, x1, y1, x2, y2, a, e, n1, n2, rho):
        self.x1 = x1
        self.y1 = y1
        self.x2 = x2
        self.y2 = y2
        self.a = a
        self.e = e
        self.n1 = n1
        self.n2 = n2
        self.l = math.sqrt((x2-x1)**2 + (y2-y1)**2)
        self.rho = rho

    def GetKe(self):
        nDofs = 4;

        c = (self.x2-self.x1)/self.l
        s = (self.y2-self.y1)/self.l
       
        cc = c * c
        ss = s * s
        cs = c * s

        s = (nDofs,nDofs)       
        ke = np.zeros(s)
        ke[0,0] = cc
        ke[0,1] = cs
        ke[1,1] = ss
        ke[1,0] = cs

        for r in range(2,4):
            for c in range(2,4):
                ke[r,c] = ke[r-2,c-2]

        for r in range(2,4):
            for c in range(0,2):
                ke[r,c] = -ke[r-2,c]

        for r in range(0,2):
            for c in range(2,4):
                ke[r,c] = -ke[r,c-2]
        
        return ke

    def GetMe(self, rho):
        nDofs = 4;

        s = (nDofs,nDofs)       
        me = np.zeros(s)
        me[0,0] = 2
        me[1,1] = 2
        me[2,2] = 2
        me[3,3] = 2

        me[0,2] = 1
        me[1,3] = 1

        me[2,0] = 1
        me[3,1] = 1

        me = me * (1/6 * rho * self.a * self.l)

        #transform the matrix to global coordinates
        s = (nDofs, nDofs)
        t = np.zeros(s)
        c = (self.x2-self.x1)/self.l
        s = (self.y2-self.y1)/self.l

        t[0,0] = c
        t[0,1] = s
        t[1,0] = -s
        t[1,1] = c

        t[2,2] = c
        t[2,3] = s
        t[3,2] = -s
        t[3,3] = c

        m = np.transpose(t) @ me @ t

        return m

Eivenvalue analysis is carried out using inverse iteration method. Inverse iteration method returns the lowest eigenvalue. To obtain other eigenvalues and eigenvectors, the Gram-Schmidt process is used to generate orthogonal trial vectors. For elaborate treatment of the subject, more specialized literature must be read before attempting to program the methods. The programs shown here have been inspired from Introduction to FInite Elements in Engineering, Tirupathi R Chandrupatla, Ashok D Belegundu, 3ed, Prentice Hall, 2002. Eigen is the most crucial class and is described below.

import numpy as np
import math

class Eigen(object):
    """Computes eigen values and eigen vectors for a given
    stiffness marix and mass matrix. Note that boundary conditions
    are required to have been already applied. The stiffness and
    mass matrices are expected to be in square format. Stiffness
    and mass matrices must be numpy matrices"""

    def __init__(self, *args, **kwargs):
        return super().__init__(*args, **kwargs)

    def Rescale(self, x):
        max = x.max()
        min = x.min()
        s = x.shape
        n = s[0]
        amax = max
        if abs(min) > max: amax = abs(min)  
        for i in range(n):
            x[i] = x[i] / max
        return x

    def RescaleEigenVectors(self, evec):
        dims = evec.shape
        ndofs = dims[0]
        for i in range(ndofs):
            evec[:,i] = self.Rescale(evec[:,i])
        return evec

    def GetOrthogonalVector(self, ndofs, trial, mg, ev,evec):
        const = 0
        s = [ndofs]
        sumcu= np.zeros(s)
        for e in range(ev  ):
            U = evec[:,e]
            const += trial @ mg @ U
            cu = [x * const for x in U]
            sumcu += cu
        trial = trial - sumcu
        return trial

    # Computes eigen values and eigen vectors using inverse iteration process
    def Solve(self,kg, mg, tolerance = 0.0000000000000001 ):
        dims = kg.shape
        ndofs = dims[0]

        # initialize the eigen vector holder and eigen values hlder
        s = (ndofs,ndofs)       
        evec = np.zeros(s)
        s = (ndofs)       
        eval = np.zeros(ndofs)

        # Step 0:
        # -------
        # prepare the initial guess of the lowest eigen vector u0
        trial = np.ones(ndofs)
        #uk_1 = np.zeros(ndofs)
        eigenvalue0 = 0
        # start the loop
        for ev in range(ndofs):
            print("Computing eigenvalue and eigen vector " + str(ev) + "... " , end="")

            converged = False            
            uk_1 = trial
            # set iteration index - k
            k = 0
            while converged == False:
                # Step 1: lets start with the procedure
                k += 1

                # if there are more eigenvalues and eigenvectors to be computed,
                # we apply Gram-Schmidt rpocess to compute an orthogonal trial
                # vector for eigenvector to obtain the next eigenvalue / eigenvector

                if ev > 0: # and ev < ndofs:
                    uk_1 = self.GetOrthogonalVector(ndofs,uk_1,mg,ev,evec)

                # Step 2: determine the right hand side
                vk_1 = mg @ uk_1 # multiply matrix m with vector u(k-1)

                # Step 3: compute u by solving equations
                uhatk = np.linalg.solve(kg,vk_1)

                # Step 4: denote vk = mu
                vhatk = mg @ uhatk

                # Step 5: estimate eigenvalue
                uhatkt = np.transpose(uhatk)
                eigenvalue = (uhatkt @ vk_1)/(uhatkt @ vhatk)

                # Step 6: normalize eigenvector
                denominator = math.sqrt(uhatkt @ vhatk)
                uk = uhatk/denominator

                # Step 7: check for tolerance
                tol = abs((eigenvalue - eigenvalue0) / eigenvalue)
                if tol <= tolerance:
                    converged = True
                    evec[:,ev] = uk
                    eval[ev] = eigenvalue
                    print("Eigenvalue = " + str(eigenvalue) + " ... done.")                        
                else:
                    eigenvalue0 = eigenvalue
                    uk_1 = uk

                    if k > 1000:
                        evec[:,ev] = uk
                        eval[ev] = eigenvalue
                        print ("could not converge. Tolerance = " + str(tol))
                        break

        self.eigenvalues = eval
        return evec

The class eigen performs the following steps:

Assume/estimate an initial trial eigenvector "u0"
Determine the right hand side vector "v" by premultiplying the trial vector by mass matrix
Solve the system of equations K u = v. Obtain u
Estimate eigenvalue at this step
Normalize the eigen vector
Check for tolerance and iterate if the tolerance criteria is not satisfied

By performing these steps, the lowest eigenvalue is obtained. Gram-Schhmidt process is adopted to generate subsequent trial vectors which are orthogonal to the previous trial vector.

Finally, the eigenvectors are plotted. For a truss example, two eigenvectors plotted as a mode shape are shown below:

Example: Vibration Modes 2 and 3

Example: Vibration Modes 4 and 5

The input file format is a simple representation of the structure. A file describing the input file format is included in the zip file attached and reproduced here for completeness:

nn, nele, nbc nodes, nloaded nodes
x1, y1
x2, y2
x3, y3
.
.
.
xnn, ynn
n11, n12, a1, e1, rho1
n21, n22, a2, e2, rho2
n31, n32, a3, e3, rho3
n41, n42, a4, e4, rho4
.
.
.
nnele1, nnele2, anele, enele, rhonele
node_bcnode1, ux restrained (1/0), uy restrained (1/0), prescibed ux, prescribed uy
node_bcnode2, ux restrained (1/0), uy restrained (1/0), prescibed ux, prescribed uy
node_bcnode3, ux restrained (1/0), uy restrained (1/0), prescibed ux, prescribed uy
.
.
.
node_nbc nodes, ux restrained (1/0), uy restrained (1/0), prescibed ux, prescribed uy
node_loadednode1, fx, fy
node_loadednode2, fx, fy
.
.
.
node_nloadednodes, fx, fy

In the format, the indicated fields and variables are:

nn = Number of nodes in the model
nele = Number of Elements in the model
nbc (same as nbc nodes)= Total number of nodes having boundary conditions (supports)
nloads (same as nloaded nodes) = Total number of nodes on which loads have been applied
x1, y1, x2, y2... are x or y coordinates describing the nodes. Thus, after initial description of the structure, the nodal coordinates are required to be specified.
Element description follows the node coordinates. Each element (member) shall be described by node 1 (node at end i), node 2 (node at end j), area of cross section, modulus of elasticity of the material and density of the material (mass per unit volume)
Supports immediately follow the element description. Each support description includes the node on which the support is applied, next is a logic of restraint in ux direction which indicates whether displacement in u direction (ux) is restrained or not (1 to indicate restrained, 0 to indicate unrestrained), followed by the value of prescribed ux. Prescribed ux must be given. If logic of restraint = 0, the prescribed value will be ignored. Next will be logic of restraint in v direction (uy). Specify 1 to indicate restrained, 0 to indicate unrestrained translation in y direction. This shall be followed by the value of prescribed uy. Again, just like ux, this value must be specified. If restraint uy = 0, prescribed uy value will be ignored.
Finally, the loads shall be specified for each loaded node. The load description includes the node on which the load is applied, load in x direction (Fx), and load in y direction (Fy)

Please note that no blank space or comments in the input file are permitted. The Input class can be very easily modified to accept comments and/or blank lines. Also note that node numbering in the input file starts from 1 and not 0 as in the code.

Points of Interest

Eigenvalue solution is an iterative process. Each iteration requires an orthogonal trial vector, for which, Gram-Schmidt process has to be carried out. This must be executed in each iteration, otherwise the solution would eventually converge to the lowest eigenvalue and corresponding eigenvector.

The program presented is intended to explain and offer a code base which can begin the journey to developing your own programs for structural dynamics simulations. The stiffness and mass matrices are stored in square format for simplicity in presentation and development of the basic code. It should not be difficult to modify the program to consider banded or skyline storage or other more efficient schemes.

History

3^rd September, 2018: Version 1.0.1. (with figures updated)