I

Inerter

V1

2023/03/25阅读：22主题：默认主题

「Python与地震工程」多自由度体系的阻尼矩阵

程序代码

import numpy as npimport scipy.linalg as splimport matplotlib.pyplot as pltplt.style.use("ggplot")plt.rcParams['font.sans-serif'] = ['Microsoft Yahei']plt.rcParams['axes.unicode_minus'] = Falseimport sysnp.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)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") # "R"    print("质量矩阵：")    print(mdof.M)    print("刚度矩阵：")    print(mdof.K)    print("阻尼矩阵：")    print(mdof.C)

结果

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


「Python与地震工程」多自由度体系的阻尼矩阵

Python

I
V1