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:
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.
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
ndof = nnodes * 2
assembly = Assembly.Assembly()
assembly. GenerateKgAndMg(ndof,elements)
assembly.GenerateLoadVector(loads)
assembly.ApplyBoundaryConditions(supports)
kg = assembly.kg
mg = assembly.mg
r = assembly.r
print("Performing linear solution... ", end="")
d = np.linalg.solve(kg,r)
print ("done.")
disp = d.reshape(nnodes,2)
output = Output.Output()
output.WriteOutputFile(inp.outfilename, nodes, elements, disp)
mfile = filename + ".m"
dispscale = 400
output.WriteModelOnMatlabFile(mfile,nodes,elements,disp)
e = Eigen.Eigen()
print("Performing eigenvalue solution... ", end="")
evec = e.Solve(kg,mg,1e-6)
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)
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)
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)
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
def Solve(self,kg, mg, tolerance = 0.0000000000000001 ):
dims = kg.shape
ndofs = dims[0]
s = (ndofs,ndofs)
evec = np.zeros(s)
s = (ndofs)
eval = np.zeros(ndofs)
trial = np.ones(ndofs)
eigenvalue0 = 0
for ev in range(ndofs):
print("Computing eigenvalue and eigen vector " + str(ev) + "... " , end="")
converged = False
uk_1 = trial
k = 0
while converged == False:
k += 1
if ev > 0:
uk_1 = self.GetOrthogonalVector(ndofs,uk_1,mg,ev,evec)
vk_1 = mg @ uk_1
uhatk = np.linalg.solve(kg,vk_1)
vhatk = mg @ uhatk
uhatkt = np.transpose(uhatk)
eigenvalue = (uhatkt @ vk_1)/(uhatkt @ vhatk)
denominator = math.sqrt(uhatkt @ vhatk)
uk = uhatk/denominator
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:
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
- 3rd September, 2018: Version 1.0.1. (with figures updated)
