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

原理见《结构动力学》

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

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

响应初始值默认为0

其他相关说明见代码

import numpy as np
from scipy import linalgclass 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法稳定系数，取值.25if 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 #阻尼系数1beta = zeta / omiga #阻尼系数2self.c_matrax = alpha * self.m_matrax + beta * k_matrix #返回阻尼矩阵a0 = 1 / beta_newmark / dt / dta1 = gama_newmark / beta_newmark / dta2 = 1 / beta_newmark / dta3 = 1 / 2 / beta_newmark -1a4 = gama_newmark / beta_newmark - 1a5 = dt / 2 * ( gama_newmark / beta_newmark - 2 )a6 = dt * ( 1 - gama_newmark )a7 = dt * gama_newmarkself.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 * dtfe = 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]```

借助 Python 进行 Newmark_beta 法计算结构响应相关推荐

matlab求刚度,matlab直接刚度法计算结构频率
matlab直接刚度法计算结构频率 syms E I K r l x cja1=-E*I*K*r^3*(cos(x)*sinh(x)+sin(x)*cosh(x))jc1=-E*I*K*r^2*sin ...
matlab用jc法计算可靠度,用MATLAB实现JC法计算结构可靠度程序
软件设计开发本栏目责任编辑 :谢媛媛 ComputerKnowledgeand Technology 电脑知识与技术第 5 卷第 29 期 (2009 年 10 月) 用 MATLAB 实现 JC 法 ...
python莱布尼茨法计算π_python圆周率计算（带进度条）
3.波尔文四次迭代式这个公式由乔纳森·波尔文和彼得·波尔文于1985年发表的. bailey-borwein-plouffe算法这个公式简称BBP公式,由David Bailey, Peter B ...
python莱布尼茨法计算π_酷叮猫少儿编程讲堂——Python 用莱布尼茨等式求π
原标题:酷叮猫少儿编程讲堂--Python 用莱布尼茨等式求π Python 用莱布尼茨等式求π 2018-08-01 德国大数学家莱布尼茨Leibniz在研究圆周率π的过程中发现一个数学公式是这样的 ...
【无标题】python利用公式法计算圆周率
# 计算圆周率(公式法) i = 1 j = 1 s = 0 print('******圆周率公式法******') for i in range(1, 100):print(f'循环{i}次,第一次 ...
Python使用筛选法计算小于给定数字的所有素数
代码思路:首先列出指定范围内所有候选数字,然后从前往后依次选择一个数字去除以后面所有数字,能够被整除的肯定不是素数,把这些数字过滤掉,然后重复这个过程,直到选择的除数大于最大数字的平方根为止.代码主要 ...
python求线段长度_python微元法计算函数曲线长度的方法
计算曲线长度,根据线积分公式: ,令积分函数f(x,y,z) 为1,即计算曲线的长度,将其微元化: 其中根据此时便可在python编程实现,给出4个例子,代码中已有详细注释,不再赘述 ''' 计算曲 ...
python计算长方形周长代码_python微元法计算函数曲线长度的方法
计算曲线长度,根据线积分公式: ,令积分函数f(x,y,z) 为1,即计算曲线的长度,将其微元化: 其中根据此时便可在python编程实现,给出4个例子,代码中已有详细注释,不再赘述 ''' 计算曲 ...
Python与Ansys apdl有限元系列二：矩阵位移法计算桁架结构
矩阵位移法计算三种桁架(静定,一次,二次超静定) 1 为了便于计算 , 桁架弹性模量与截面积乘积 E A = 1 , 所有杆件均为二力杆 \qquad为 ...

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

借助 Python 进行 Newmark_beta 法计算结构响应相关推荐

最新文章

热门文章