I

Inerter

V1

2023/03/08阅读：9主题：默认主题

「Python与地震工程」瑞利-里兹法与子空间迭代法

程序代码

# mdof.pyimport numpy as npimport scipy.linalg as splfrom numpy.linalg import solveimport sysnp.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 = Xif __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]]


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 %

Python

I
V1