「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 %

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