I
Inerter
V1
2023/03/08阅读:9主题:默认主题
「Python与地震工程」瑞利-里兹法与子空间迭代法

「Python与地震工程」瑞利-里兹法与子空间迭代法
原理
瑞利-里兹法可以用较小的代价求解结构前几阶自振频率。通过子空间迭代可以提高求解精度。
瑞利-里兹法
可将计算瑞利商所需的变形向量写成若干向量线性组合的形式
代入瑞利商的表达式中可得
瑞利商是真实自振圆频率的上限。上述瑞利商的最小值将是对圆频率的最佳估计。
其中
上式将N阶特征值问题变为L阶特征值问题。
子空间迭代法
子空间迭代法可以提高瑞利-里兹法的求解精度。
包含所有特征值向量的结构广义特征值问题可表示为
其中, 是特征向量组成的矩阵
是特征值组成的对角矩阵
广义特征值问题等式两边同时右乘
令 ( 相当于 中各阶特征值向量除以相应的特征值,即两个矩阵中的列向量间保持着比例关系),则
子空间迭代法的计算过程
确定初始向量矩阵
其中 为所需计算的特征值数量, 。
解线性方程组计算 矩阵
以矩阵 中的向量为基底,应用瑞利-里兹法。即令 ,求解以下广义特征值问题
其中
求解得到特征向量矩阵
进而得到下一次迭代的特征向量矩阵
重复以上过程,直至两次迭代计算的最高阶特征值的误差小于容许值。
程序代码
# mdof.py
import numpy as np
import scipy.linalg as spl
from numpy.linalg import solve
import sys
np.set_printoptions(precision=4, threshold=sys.maxsize,
linewidth=150)
class MDOF():
def __init__(self, m, k, c, zeta=0.05) -> None:
self.m = m
self.k = k
self.c = c
self.zeta = zeta
self.ndof = len(m)
def build_matrix(self, damp_mat="D", i=1, j=2):
'''
建立质量、刚度、阻尼矩阵和荷载向量\\
参数含义:
damp_mat: 固有阻尼矩阵的建立方法。“D”为直接构造(默认),“R”为瑞利阻尼。
i,j: 确定瑞利阻尼系数的两阶振型序号。默认i=1, j=2。
'''
self.M = np.diag(self.m)
k1 = np.zeros_like(self.k)
k1[:-1] = self.k[1:]
self.K = np.diag(self.k+k1) + np.diag(-self.k[1:],1) + np.diag(-self.k[1:],-1)
c1 = np.zeros_like(self.c)
c1[:-1] = self.c[1:]
self.C = np.diag(self.c+c1) + np.diag(-self.c[1:],1) + np.diag(-self.c[1:],-1)
self.E = np.ones(self.ndof)
if self.zeta > 0:
self.mode_analyze()
if damp_mat == "D":
M_n = np.diag(self.phi.T@self.M@self.phi)
C_ = np.diag(2.0*self.zeta*self.omega/M_n)
self.C += self.M@self.phi@C_@self.phi.T@self.M
elif damp_mat == "R":
w_i, w_j = self.omega[i-1], self.omega[j-1]
a_0, a_1 = self.zeta*2.0*w_i*w_j/(w_i+w_j), self.zeta*2.0/(w_i+w_j)
self.C += a_0*self.M + a_1*self.K
else:
self.C += np.zeros_like(self.M)
def mode_analyze(self, check_orthogonality=False):
omega2, self.phi = spl.eigh(self.K, self.M)
self.omega = np.sqrt(omega2)
self.period = 2.0*np.pi/self.omega
# 验证正交性
if check_orthogonality:
print("正交性验证")
print("ΦMΦ:")
print(self.phi.T @ self.M @ self.phi)
print("ΦKΦ:")
print(self.phi.T @ self.K @ self.phi)
def calc_eigen_Rayleigh_Ritz(self, Phi):
p_i = self.M @ Phi
Phi = solve(self.K, p_i)
Phi = Phi/Phi[-1,:]
K_s = Phi.T @ self.K @ Phi
M_s = Phi.T @ self.M @ Phi
omega2, Z = spl.eigh(K_s, M_s)
self.omega = np.sqrt(omega2)
self.period = 2.0*np.pi/self.omega
self.phi = Phi @ Z
def calc_eigen_subspace_iteration(self, num_vecs, num_iters=20, tol=1e-6):
X = np.zeros((self.ndof, num_vecs))
X[:,0] = 1.0
X[1:num_vecs,1:num_vecs] = np.eye(num_vecs-1)
Y = solve(self.K, self.M @ X)
K_s = Y.T @ self.K @ Y
M_s = Y.T @ self.M @ Y
omega2, Z = spl.eigh(K_s, M_s)
X = Y @ Z
omega2_n = omega2[-1]
for i in range(num_iters-1):
Y = solve(self.K, self.M @ X)
K_s = Y.T @ self.K @ Y
M_s = Y.T @ self.M @ Y
omega2, Z = spl.eigh(K_s, M_s)
X = Y @ Z
error = abs(omega2_n - omega2[-1])/omega2[-1]
if error < tol:
break
omega2_n = omega2[-1]
self.omega = np.sqrt(omega2)
self.period = 2.0*np.pi/self.omega
self.phi = X
if __name__ == '__main__':
ndof = 100
m = np.ones(ndof)
k = 1600.0*m
c = np.zeros(ndof)
zeta = 0.05 # 固有阻尼比
mdof = MDOF(m, k, c, zeta)
mdof.build_matrix("D")
print("实模态分析:")
mdof.mode_analyze()
for i in range(ndof):
print("第 %d 阶 自振圆频率:%f rad/s"%(i+1, mdof.omega[i]))
print("瑞利-里兹法求解自振频率:")
num_vecs = 6
omega = mdof.omega[:num_vecs]
vecs = np.zeros((ndof, num_vecs))
for i in range(1,num_vecs+1):
for j in range(1,ndof+1):
vecs[j-1,i-1] = 1.0 - np.cos((i-1/2)*np.pi/ndof*j)
# vecs[j-1,i-1] = np.sin((i-1/2)*np.pi/ndof*j)
mdof.calc_eigen_Rayleigh_Ritz(vecs)
print("| 阶数 | 圆频率 | 误差 |")
print("|------|-----------|-----|")
for i in range(num_vecs):
error = (mdof.omega[i] - omega[i])/omega[i]
print("| %d | %.6f | %.6f %% |"%(i+1, mdof.omega[i], error*100.0))
print("子空间迭代法求解自振频率:")
num_vecs = 6
mdof.calc_eigen_subspace_iteration(num_vecs,10,1e-20)
print("| 阶数 | 圆频率 | 误差 |")
print("|------|-----------|-----|")
for i in range(num_vecs):
error = (mdof.omega[i] - omega[i])/omega[i]
print("| %d | %.6f | %.6f %% |"%(i+1, mdof.omega[i], error*100.0))
结果
graph TD;
A((m))-- k ---B((m))-. k .-C((m));
C-- k ---D((m))-- k ---E[[ground]]
以100自由度结构为例,求解前6阶圆频率。 前6阶圆频率的精确解为:
阶数 | 圆频率 rad/s |
---|---|
1 | 0.625186 |
2 | 1.875406 |
3 | 3.125167 |
4 | 4.374166 |
5 | 5.622095 |
6 | 6.868651 |
瑞利-里兹法
-
基底向量取为
其中
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.646584 | 3.422616 % |
2 | 1.945637 | 3.744823 % |
3 | 3.236251 | 3.554469 % |
4 | 4.586019 | 4.843285 % |
5 | 5.868646 | 4.385393 % |
6 | 8.493746 | 23.659588 % |
-
基底向量取为
其中
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.625187 | 0.000205 % |
2 | 1.875441 | 0.001869 % |
3 | 3.125335 | 0.005364 % |
4 | 4.374653 | 0.011135 % |
5 | 5.623238 | 0.020333 % |
6 | 6.871211 | 0.037265 % |
计算基底向量的表达式是均匀纯剪切梁的振型函数。串联多自由度体系可视为离散化的纯剪切梁,因此计算所得圆频率误差很小。
子空间迭代法
-
子空间向量6个,迭代1次
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.627964 | 0.444278 % |
2 | 4.121268 | 119.753399 % |
3 | 21.554322 | 589.701354 % |
4 | 41.944486 | 858.914012 % |
5 | 61.238154 | 989.240795 % |
6 | 75.030882 | 992.367031 % |
-
子空间向量6个,迭代3次:
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.625186 | 0.000000 % |
2 | 1.877667 | 0.120552 % |
3 | 3.523650 | 12.750759 % |
4 | 10.743097 | 145.603348 % |
5 | 33.503885 | 495.932382 % |
6 | 59.208954 | 762.017176 % |
-
子空间向量6个,迭代10次:
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.625186 | 0.000000 % |
2 | 1.875406 | 0.000000 % |
3 | 3.125167 | 0.000000 % |
4 | 4.374166 | 0.000017 % |
5 | 5.623881 | 0.031768 % |
6 | 7.028587 | 2.328482 % |
-
子空间向量6个,迭代20次:
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.625186 | 0.000000 % |
2 | 1.875406 | 0.000000 % |
3 | 3.125167 | 0.000000 % |
4 | 4.374166 | 0.000000 % |
5 | 5.622095 | 0.000000 % |
6 | 6.868840 | 0.002744 % |
-
若将子空间向量增加到12个,10次迭代即可使前6阶圆频率达到足够精度:
阶数 | 圆频率 | 误差 |
---|---|---|
1 | 0.625186 | 0.000000 % |
2 | 1.875406 | 0.000000 % |
3 | 3.125167 | 0.000000 % |
4 | 4.374166 | 0.000000 % |
5 | 5.622095 | 0.000000 % |
6 | 6.868654 | 0.000044 % |
转载本文请注明出处。科研成果中使用本文代码请引用作者课题组的研究论文(请戳“阅读原文”)。
作者介绍
I
Inerter
V1