「Python与地震工程」荷载相关的里兹向量

原理

振型叠加法可以认为是以结构各阶振型为基底来描述结构各自由度的动力响应这一多维向量的变化。从数学角度来说，这样的基底其实不是唯一的（类似于空间坐标系不唯一）。这里即将讨论的里兹向量就是这样的一组基底向量，或者称为一种特殊的振型。它具有正交性的特征，能替代结构振型进行振型叠加的运算。由于它考虑了荷载的影响，因此可以使用较少的振型得到较高精度的结果。但当结构收到的动力荷载加载位置很多时（比如将非线性阻尼器作为等效荷载考虑时）则需要大量的里兹向量，使用时需要注意。

荷载相关的里兹向量

给定结构质量矩阵 、刚度矩阵 和荷载向量

求解结构的静力位移并进行归一化处理即得到第一阶里兹向量（Ritz Vector）

将第一阶里兹向量对应的惯性力施加到结构上，求解得结构的静力位移

中包含有与 线性相关的部分，需要剔除掉，故假设

其中， 是第二阶里兹向量（未归一化）。

两边同乘

由于 已经过归一化处理，故 ；考虑里兹向量的正交性， 。代入上式得

则有

归一化处理后得第二阶里兹向量

当需要计算第n阶（n>2）里兹向量时，将第n-1阶里兹向量对应的惯性力施加到结构上，求解得结构的静力位移

中包含有与 线性相关的部分，需要剔除掉，故假设

其中， 是第二阶里兹向量（未归一化）。

不失一般性，两边同乘 （$i < p>

<>

由于 已经过归一化处理，故 ；考虑里兹向量的正交性， ， 。代入上式得

则有

令

归一化处理后得第 n 阶里兹向量

假设前 L 阶里兹向量已得到，令 ，则有 。即所得向量关于质量矩阵是正交的。

以 作为广义瑞利-里兹法的基底向量，解特征值问题

得到特征向量（归一化）

显然所得向量是关于刚度矩阵正交的。

最终得到荷载相关的、满足正交条件的里兹向量如下

程序代码

# mdof.py
import numpy as np
import scipy.linalg as spl
import matplotlib.pyplot as plt

plt.style.use("ggplot")
plt.rcParams['font.sans-serif'] = ['Microsoft Yahei']
plt.rcParams['axes.unicode_minus'] = False

import sys
np.set_printoptions(precision=4, threshold=sys.maxsize,
                    linewidth=150)


def Ritz_vector(M,K,p,L=0):
    
    N = np.size(M,0)
    L = L if L>0 else N
    X = np.zeros((N,L))

    x = spl.solve(K,p)
    x = x/np.sqrt(x.T @ M @ x)

    X[:,0] = x

    for i in range(L-1):
        x_p = X[:,i]
        y = spl.solve(K,M @ x_p)
        
        X_ = X[:,:i+1]
        x = y - X_ @ X_.T @ M @ y
        x = x/np.sqrt(x.T @ M @ x)
        X[:,i+1] = x

    omega2, Z = spl.eigh(X.T @ K @ X)

    F = X @ Z

    return omega2, F


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_Ritz_vector(self,load_vec,num_vecs=0,check_orthogonality=False):
        omega2, self.phi = Ritz_vector(self.M,self.K,load_vec,num_vecs)
        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)


if __name__ == '__main__':

    ndof = 5
    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("里兹向量：")
    load_vec = np.ones(ndof)
    # load_vec = np.zeros(ndof)
    # load_vec[-1] = 1
    # load_vec[-2] = -2
    for num_vecs in range(1, ndof+1):
        print("计算 %d 阶里兹向量："%num_vecs)
        mdof.calc_Ritz_vector(load_vec, num_vecs)
        for i in range(num_vecs):
            print("第 %d 阶 特征圆频率：%f rad/s"%(i+1, mdof.omega[i]))

结果

以5自由度体系为例。

graph TD
A((m))-- k ---B((m))-- k ---C((m))
C-- k ---D((m))-- k ---E((m))
E-- k ---F[[ground]]

里兹向量对应特征圆频率（均匀侧向力）

阶数	实模态分析	里兹向量（1阶）	里兹向量（2阶）	里兹向量（3阶）	里兹向量（4阶）	里兹向量（5阶）
1	11.385187	11.451967	11.385239	11.385187	11.385187	11.385187
2	33.233201	-	34.772861	33.253419	33.233230	33.233201
3	52.388859	-	-	55.618607	52.459430	52.388859
4	67.300283	-	-	-	69.393876	67.300283
5	76.759438	-	-	-	-	76.759438

非对角线元素均为接近0的极小值，验证了里兹向量的正交性。

里兹向量对应特征圆频率（顶点集中力）

阶数	实模态分析	里兹向量（1阶）	里兹向量（2阶）	里兹向量（3阶）	里兹向量（4阶）	里兹向量（5阶）
1	11.385187	12.060454	11.386225	11.385187	11.385187	11.385187
2	33.233201	-	37.159701	33.305298	33.233344	33.233201
3	52.388859	-	-	57.332441	52.526596	52.388859
4	67.300283	-	-	-	69.856534	67.300283
5	76.759438	-	-	-	-	76.759438

里兹向量对应特征圆频率（顶部两层施加集中力）

顶层质点施加集中力1，次顶层质点施加集中力-2

阶数	实模态分析	里兹向量（1阶）	里兹向量（2阶）	里兹向量（3阶）	里兹向量（4阶）	里兹向量（5阶）
1	11.385187	14.322297	11.386293	11.385189	11.385187	11.385187
2	33.233201	-	59.376652	37.653626	33.262472	33.233201
3	52.388859	-	-	64.124680	53.078318	52.388859
4	67.300283	-	-	-	71.196135	67.300283
5	76.759438	-	-	-	-	76.759438

里兹向量可以用较少阶特征向量覆盖较大的频率范围，但有可能“跳过”某些频率。

当待求解里兹向量数等于结构自由度数时，里兹向量与实模态振型是一样的。

转载本文请注明出处。「Python与地震工程」荷载相关的里兹向量