使用 Python 进行平面桁架的模态分析
5.00/5 (6投票s)
Python 中结构动力响应的基本模态分析。
引言
结构在动态载荷下会振动。这些振动对于分析师和设计者来说至关重要，因为动态载荷通常会引起比静载荷大得多的结构响应。静力分析相对简单，有现成的解决方案。动力学分析需要一套不同的线性代数运算。为了获得结构上的动态载荷，需要进行模态分析。本文以一个二维桁架问题为例，描述了对结构进行模态分析的步骤，并提供了Python代码。
Python代码结构清晰，易于理解，并允许修改代码以实现更多附加功能。完整的Python代码已附在此处，供您参考，以帮助您更好地理解文章。
背景
虽然本文旨在提供完整的信息，但读者应具备结构力学的基础知识。另外，建议读者阅读一些关于结构力学和结构动力学的基本教材。对有限元分析一文的快速回顾：https://codeproject.org.cn/Articles/1207954/Finite-Element-Analysis-in-VB-Net 会非常有帮助。尽管如此，下文将介绍为动力学系统求解的基本方程。
对于如下所示具有分布式质量的实体
动能由下式给出
$ T = \frac{1}{2} \int_V{u^T u \rho DVD} $
其中向量u的元素是变形，由u = Nq给出
对于结构动力学问题，q向量是随时间变化的，而N是结构模型获得的形函数。那么单元的动能T为
$ T_e = \frac{1}{2} \dot{q}^T \left[\int_e{\rho N^T N dV }\right]\dot{q} $
方括号内的表达式是单元质量矩阵
$ m_e = \int_e{\rho N^T N dV} $
由于此质量矩阵与单元形函数一致，因此称为一致质量矩阵。将整个域的动能相加（对每个单元求和）得到
$ 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} $
势能由下式给出
$ \Pi = \frac{1}{2} Q^T K Q - Q^T F $
使用拉格朗日量L；
$ L = T - \Pi $
这最终导致了运动方程
$ M \ddot{Q} + K Q = F $
对于自由振动，力F = 0。因此，
$ M \ddot{Q} + K Q = 0 $
这最终导向了广义特征值问题
$ K U = \lambda M U $
其中，\(U\)是特征向量，代表与振动模态相对应的归一化变形；\(\lambda\)是对应的特征值，是圆频率\(\omega\)的平方。模态分析涉及求解此处给出的特征值问题。该解不是直接的，需要迭代过程。
Using the Code
附带的代码是用Visual Studio 2017编写的，但应该可以在任何Python 3.x解释器上运行。整个Visual Studio解决方案已在此处附上。代码需要numpy和matplotlib库。已尽力遵守单一职责原则。为每种任务类型创建了单独的类。这里解释了最重要的类。开始之前，您需要添加numpy和matplotlib环境。
主文件是PlaneTruss.py。该文件导入所需的各种模块和库。下面显示了PlaneTruss.py文件中调用每个例程以执行特征值求解的部分：您可以在附件中看到完整的文件。
# 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.")
此处执行的步骤如下：
- 读取输入。
- 生成单元刚度和质量矩阵。
- 组装单元刚度和质量矩阵，形成整体刚度和整体质量矩阵。
- 生成载荷向量。
- 施加边界条件。请注意，在此步骤中仅修改刚度矩阵。
- 检索整体矩阵。
- 进行静力分析。这主要是为了确认结构已正确建模，并且边界条件已适当地应用。
- 进行特征值分析，以获得特征值和特征向量，它们给出了模态形状。
每个实体由其自己的对象管理。例如，Element 类负责管理完全定义一个单元所需的变量。本例考虑二维桁架单元，因此二维桁架的刚度和质量矩阵开发如下：
$\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} $
其中，*A*是截面积，*E*表示桁架构件材料的弹性模量，*l*是单元/构件的长度。*c*和*s*表示构件的方向余弦和正弦。
一致质量矩阵由下式给出：
$\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} $
上面显示的矩阵由下面的Element 类生成
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
特征值分析使用反迭代法进行。反迭代法返回最低的特征值。为了获得其他特征值和特征向量，使用Gram-Schmidt过程生成正交试探向量。有关该主题的详细处理，在尝试编程这些方法之前，必须阅读更专业的文献。此处显示的程序灵感来自于《工程有限元导论》，Tirupathi R Chandrupatla，Ashok D Belegundu，第3版，Prentice Hall，2002。Eigen是最关键的类，下面将进行说明。
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
eigen类执行以下步骤：
- 假设/估计一个初始试探特征向量"u0"
- 通过将试探向量乘以质量矩阵得到右侧向量"v"
- 求解方程组K u = v。得到u
- 估计此步骤的特征值
- 归一化特征向量
- 检查容差，如果未满足容差标准，则迭代
通过执行这些步骤，得到最低的特征值。采用Gram-Schhmidt过程生成与前一个试探向量正交的后续试探向量。
最后，绘制特征向量。对于桁架示例，下面显示了两个作为模态形状绘制的特征向量。
输入文件格式是对结构的简单表示。描述输入文件格式的文件包含在附带的zip文件中，并在此处重新生成以供完整性考虑。
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
在该格式中，指示的字段和变量是：
- nn = 模型中的节点数
- nele = 模型中的单元数
- nbc（与nbc节点相同）= 具有边界条件（支座）的总节点数
- nloads（与nloaded节点相同）= 施加了载荷的总节点数
- x1, y1, x2, y2...是描述节点的x或y坐标。因此，在对结构进行初始描述后，需要指定节点坐标。
- 单元描述 follows 节点坐标。每个单元（构件）应由节点1（端点i处的节点）、节点2（端点j处的节点）、截面积、材料的弹性模量和材料密度（每单位体积的质量）描述。
- 支座紧随单元描述之后。每个支座描述包括施加支座的节点，接着是ux方向约束的逻辑，该逻辑指示ux方向（ux）的位移是否受到约束（1表示约束，0表示不约束），后跟规定的ux值。必须给出规定的ux。如果约束逻辑=0，则规定的值将被忽略。接下来是v方向（uy）约束的逻辑。指定1表示在y方向上的平移受到约束，0表示不约束。接着是规定的uy*的值。同样，与ux*一样，必须指定此值。如果约束uy=0，则规定的uy值将被忽略。
- 最后，应为每个加载节点指定载荷。载荷描述包括施加载荷的节点、x方向的载荷（Fx）和y方向的载荷（Fy）。
请注意，输入文件中不允许有空格或注释。Input类可以非常轻松地修改以接受注释和/或空行。另请注意，输入文件中的节点编号从1开始，而不是像代码中那样从0开始。
关注点
特征值解是一个迭代过程。每次迭代都需要一个正交试探向量，为此必须执行Gram-Schmidt过程。否则，解最终将收敛到最低的特征值和相应的特征向量。
提供的程序旨在解释并提供一个代码库，该代码库可以开始开发自己的结构动力学仿真程序的旅程。为便于展示和开发基本代码，将刚度和质量矩阵以方形格式存储。将其修改为考虑带状或天际线存储或其他更有效的方案并不困难。
历史
- 2018年9月3日：版本1.0.1。（图片已更新）