I

Inerter

V1

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

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

Collapsed building (Generated via OpenAI)
Collapsed building (Generated via OpenAI)

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

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

分类:

数学

标签:

Python

作者介绍

I
Inerter
V1