|
原理见《结构动力学》
已做面向对象处理，可直接调用
响应初始值默认为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]
```
?
|