借助 Python 进行 Newmark_beta 法计算结构响应

原理见《结构动力学》

已做面向对象处理，可直接调用

注意：该脚本面向层数大于1的 ''糖葫芦串'' 简化结构模型

响应初始值默认为0

其他相关说明见代码

import numpy as np
from scipy import linalg

class newmark_beta:
    # 该方法将结构简化为框架结构
    # self.m_matrax 返回质量矩阵
    # self.k_matrix 返回刚度矩阵
    # self.w_n 返回各阶基频(Hz)
    # self.c_matrax 返回阻尼矩阵
    # self.t 返回时刻矩阵[i,j],i=1，j为对应时刻
    # self.d 返回位移矩阵[i,j],i为对应层数，j为对应时刻
    # self.v 返回速度矩阵[i,j],i为对应层数，j为对应时刻
    # self.a 回加速度矩阵[i,j],i为对应层数，j为对应时刻
    def __init__(self,m,k,zeta,nt,dt,force,gama_newmark = 0.5, beta_newmark = 0.25):
        # m 简化为悬臂梁后的质量列表，由下至上
        # k 简化为悬臂梁后的刚度列表，由下至上
        # zeta 阻尼比
        # nt 总步数
        # dt 子步长
        # force 外力矩阵，由下至上，i = 层数，j = 子步数
        # gama_newmark: newmark_beta法稳定系数，取值0.5
        # beta_newmark: newmark_beta法稳定系数，取值.25
        if len(m) != len(k):
            print('质量和刚度列表长度不同')
            quit()

        num = len(k)
        self.m_matrax = np.diag(m) #返回质量矩阵
        
        k_matrix = np.zeros((num,num))
        for i in range(num):
            if i == 0:
                k_matrix[i][i] = k[0] + k[1]
                k_matrix[i][i+1] = -1 * k[1]
            elif i == num-1:
                k_matrix[i][i] = k[-1]
                k_matrix[i][i-1] = -1 * k[-1]
            else:
                k_matrix[i][i-1] = -1 * k[i]
                k_matrix[i][i] =  k[i] + k[i+1]
                k_matrix[i][i+1] = -1 * k[i+1]
        self.k_matrix = k_matrix #返回刚度矩阵

        self.w_n = 1/np.sqrt(linalg.eigvals(self.m_matrax,k_matrix).real)/(2 * np.pi) #各阶基频(Hz)

        omiga = self.w_n[0] * 2 * np.pi #圆频率单位的基频  
        alpha = omiga * zeta #阻尼系数1
        beta = zeta / omiga #阻尼系数2
        self.c_matrax = alpha * self.m_matrax + beta * k_matrix #返回阻尼矩阵

        a0 = 1 / beta_newmark / dt / dt
        a1 = gama_newmark / beta_newmark / dt
        a2 = 1 / beta_newmark / dt
        a3 = 1 / 2 / beta_newmark -1
        a4 = gama_newmark / beta_newmark - 1
        a5 = dt / 2 * ( gama_newmark / beta_newmark - 2 )
        a6 = dt * ( 1 - gama_newmark )
        a7 = dt * gama_newmark

        self.t = np.zeros([1, nt+1])  #返回时刻矩阵
        self.d = np.zeros([num , nt+1])  #返回位移矩阵
        self.v = np.zeros([num , nt+1])  #返回速度矩阵
        self.a = np.zeros([num , nt+1])  #返回加速度矩阵
        ke = k_matrix + a0 * self.m_matrax + a1*self.c_matrax #等效刚度

        for j in range(1,nt+1):  #迭代计算
            self.t[0][j]  = j * dt
            fe = np.reshape(force[:,j-1],[num,1]) \
                + np.dot(self.m_matrax , np.reshape((a0*self.d[:,j-1] + a2*self.v[:,j-1] + a3*self.a[:,j-1]),[num,1])) \
                + np.dot(self.c_matrax , np.reshape((a1*self.d[:,j-1] + a4*self.v[:,j-1] + a5*self.a[:,j-1]),[num,1]))
            self.d[:,j] = np.reshape(np.dot(np.linalg.inv(ke) , fe),[1,num])
            self.a[:,j] = a0 * (self.d[:,j] - self.d[:,j-1]) - a2 * self.v[:,j-1] - a3*self.a[:,j-1]
            self.v[:,j] = self.v[:,j-1] + a6 * self.a[:,j-1] + a7*self.a[:,j]


```