「Python与地震工程」复模态分析

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)

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.dot(self.M).dot(self.phi))
                C_ = np.diag(2.0*self.zeta*self.omega/M_n)
                self.C += self.M.dot(self.phi).dot(C_).dot(self.phi.T).dot(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 complex_mode_analyze(self, check_orthogonality=False):
        self.O = np.zeros((self.ndof, self.ndof))

        self.A = np.block( [
            [self.C, self.M], 
            [self.M, self.O] ] )
        self.B = np.block( [
            [self.K, self.O], 
            [self.O,-self.M] ] )
        
        lamda, Phi = spl.eig(-self.B, self.A)
        # eig 默认将复特征值按模从大到小排列，按结构分析惯例将其逆序排列
        self.lamda = lamda[::-1]
        self.Phi = np.flip(Phi, axis=1)
        
        # 计算自振圆频率及阻尼比
        tmp = self.lamda[::2]
        self.omega = abs(tmp)
        self.zeta = -tmp.real/self.omega
        
        # 验证正交性
        if check_orthogonality:
            print("正交性验证")
            print("ΦAΦ:")
            print(self.Phi.T @ self.A @ self.Phi)
            print("ΦBΦ:")
            print(self.Phi.T @ self.B @ self.Phi)


if __name__ == '__main__':

    ndof = 3
    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("复模态分析：")
    mdof.complex_mode_analyze(True)
    for i in range(ndof*2):
        print("第 %d 阶 复特征值："%(i+1), mdof.lamda[i])
    
    for i in range(ndof):
        print("第 %d 阶 自振圆频率：%f rad/s    阻尼比：%f"%(i+1, mdof.omega[i], mdof.zeta[i]))