2017年11月伊朗克尔曼沙省7.3级地震（AFP Photo/Atta Kenare/File）

「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 @ 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)

结果

以3自由度体系为例。

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

对于阻尼矩阵是否实现了预设的阻尼比可通过复模态分析验证，参见复模态分析一文中的探讨。

柯西阻尼矩阵的建立留作读者练习。

转载本文请注明出处。

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

科研成果中使用本文代码请引用作者课题组的研究论文：

作者课题组发表的研究论文列表（持续更新中……）