「Python与地震工程」多自由度体系基本自振频率估算——瑞利法

原理

当仅需要多自由度体系的第一阶自振周期时可以使用计算量较小的估算方法。其中最具代表性的就是瑞利法。

瑞利商

无阻尼体系发生自由振动时符合机械能守恒定律。

多自由度体系的动能与势能可按下式计算

自由振动时

根据机械能守恒定律

即

瑞利商

瑞利商是真实自振圆频率的上限。利用上式估算结构自振圆频率的方法称为瑞利法。

瑞利商的迭代修正

要找出计算瑞利商所需合适的变形向量并不容易。实用中可通过迭代修正的方式得到更合理的向量。

• 初选变形向量
• 将变形向量 相关的惯性力施加到结构上
• 求解静力问题 ，得到新的变形向量 ，并计算瑞利商

• 重复以上两个步骤若干次，得到变形向量 ，及瑞利商

• 当前后两次迭代计算的瑞利商误差小于容许值时可结束计算。

该方法实际上等价于矩阵逆迭代法。

程序代码

# mdof.py
import numpy as np
import scipy.linalg as spl

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_Rayleigh_quotient(self, vec, num_iters=10, tol=1e-6):
        phi = np.copy(vec)
        omega2 = (phi.T@self.K@phi)/(phi.T@self.M@phi)
        
        for i in range(num_iters):
            p = self.M @ phi
            phi = solve(self.K, p)
            
            omega2_ = omega2
            omega2 = (phi.T@self.K@phi)/(phi.T@self.M@phi)
            error = abs(omega2_ - omega2)/omega2
            if error < tol:
                break
            
        self.omega_R = np.sqrt(omega2)
        self.phi_R = phi



if __name__ == '__main__':

    ndof = 12
    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")
    mdof.mode_analyze()

    print("瑞利法求解基频：")
    vec = np.ones(ndof)
    # vec = np.random.rand(ndof)
    print("| 迭代次数 | 圆频率 | 误差 |")
    print("|------|-----------|-----|")
    for i in range(6):
        mdof.calc_Rayleigh_quotient(vec, i)
        error = (mdof.omega_R - mdof.omega[0])/mdof.omega[0]
        print("| %d | %.6f | %.6f %% |"%(i, mdof.omega_R, error*100.0))

结果

12自由度结构

精确解为 。

瑞利法初始向量取为全1向量：

迭代次数	圆频率	误差
0	11.547005	129.871593 %
1	5.055601	0.644202 %
2	5.023592	0.006972 %
3	5.023246	0.000086 %
4	5.023242	0.000001 %
5	5.023242	0.000000 %

瑞利法初始向量取为随机向量：

迭代次数	圆频率	误差
0	15.427922	207.130806 %
1	5.038904	0.311799 %
2	5.023398	0.003120 %
3	5.023244	0.000039 %
4	5.023242	0.000000 %
5	5.023242	0.000000 %

24自由度结构

精确解为 。

瑞利法初始向量取为全1向量：

迭代次数	圆频率	误差
0	8.164966	218.430731 %
1	2.580914	0.654708 %
2	2.564304	0.006924 %
3	2.564128	0.000084 %
4	2.564126	0.000001 %
5	2.564126	0.000000 %

瑞利法初始向量取为随机向量：

迭代次数	圆频率	误差
0	34.966516	1263.681553 %
1	2.580012	0.619549 %
2	2.564302	0.006867 %
3	2.564128	0.000085 %
4	2.564126	0.000001 %
5	2.564126	0.000000 %

转载本文请注明出处。科研成果中使用本文代码请引用作者课题组的研究论文（请戳“阅读原文”）。