NewMark法求解动力学响应

问题描述：两个质量块m1 m2由三个弹簧连接k1 k2 k3

二自由度系统，m1为1kg m2为2kg k1为1000N/m k2为1000N/m k3为1000N/m

刚度阵为[[k1+k2,-k2],[-K2,K2+K3]]
质量阵用集中质量阵表示。
阻尼采用比例阻尼。
newmark法有两个参数β与γ(有的地方称呼不一样)。
该方法将连续时间进行离散，从初始条件x0 v0 a0开始，往后迭代求解下一个时间点的x v a。从而求出这一系列时间点处的响应。
调节两个参数可能会导致结果不一样，具体情况具体考虑。

===================================================

import numpy as np
import numpy.linalg as linalgm = np.array([[1,0],[0,2]])
c = 10*m   #修改阻尼观察结果
k = np.array([[2000,-1000],[-1000,2000]])
f0 = np.array([[0,0]]).T
x0 = np.array([[0,0]]).T
v0 = np.array([[0,0]]).T
m_inv = linalg.inv(m)# a0 = m'*(f0-c*v0-k*x0)
a0 = np.dot(m_inv,(f0-np.dot(c,v0)-np.dot(k,x0)))tm = 2 #总时间
dt = 0.01 #时间步长
nt = int(tm/dt) #步数
T = np.arange(0,tm,dt) #生成时间序列  共nt个元素 所以x v a 也共nt个#构建盛放计算结果的矩阵
empty = np.zeros([2,nt-1]) #定义一个空的矩阵 来放后续计算的x v a
x = np.hstack((x0,empty))
v = np.hstack((v0,empty))
a = np.hstack((a0,empty))β = 0.25
γ = 0.25α0 = 1/(γ*dt*dt)
α1 = β/(γ*dt)
α2 = 1/(γ*dt)
α3 = 1/(2*γ)-1
α4 = β/γ-1
α5 = dt/2*(β/γ-2)
α6 = dt-β*dt
α7 = β*dtf = np.array([[0,10]]).T #在第二个质量块上施加10Nfor i in np.arange(1,nt,1):   #i从1到nt-1"""因为python自己数组特性的缘故,我在此使用了reshape函数"""xi_1 = x[:,i-1].reshape(2,1)vi_1 = v[:,i-1].reshape(2,1)ai_1 = a[:,i-1].reshape(2,1)k_ = k + α0*m + α1*c f_ = f + np.dot(m,(α0*xi_1+α2*vi_1+α3*ai_1)) + np.dot(c,(α1*xi_1+α4*vi_1+α5*ai_1))xi = np.dot(linalg.inv(k_),f_) #位移=等效荷载/等效刚度x[:,i] = xi.reshape(2,)ai = α0*(xi-xi_1)-α2*vi_1-α3*ai_1 a[:,i] = ai.reshape(2,)vi = vi_1+α6*ai_1+α7*aiv[:,i] = vi.reshape(2,)#=====================================================
# 画图
import matplotlib.pyplot as pltplt.figure(1)
plt.plot(T,x[0,:],'r--',label='x1')
plt.plot(T,x[1,:],'b-',label='x2')
plt.legend()plt.figure(2)
plt.plot(T,v[0,:],'r--',label='v1')
plt.plot(T,v[1,:],'b-',label='v2')
plt.legend()plt.figure(3)
plt.plot(T,a[0,:],'r--',label='a1')
plt.plot(T,a[1,:],'b-',label='a2')
plt.legend()

NewMark法求解动力学响应相关推荐

【源码】采用空间有限元法和时间域Newmark法求解梁的振动问题
求解欧拉-伯努利梁的振动问题(包括静载自由简支梁). Solve the vibration of Euler-Bernoulli beam (including calmped-free and s ...
newmark法 matlab,newmark法和wilson法求解单自由度体系加速度反应谱
fid=fopen('77(0.01-4172).txt','r'); %以下为newmark法求解单自由度体系反应谱 [Acceleration,count]=fscanf(fid,'%g'); d ...
matlab 纽马克激励,用Newmark方法计算系统的动力学响应的matlab程序
请大家帮忙看看这个程序有什么问题?用Newmark方法计算系统的动力学响应,结果大的惊人. function[Q,V,AA]=newmarkb E=2.1e11;P=7850;D1=0.405;d1= ...
卷积法求解系统的零状态响应_连续LTI系统零状态响应求解方法的分析
连续 LTI 系统零状态响应求解方法的分析张淑敏 [摘要] [摘要]零状态响应是电子技术相关课程中的一个重要概念,本文将通过对时域分析法和(复)频域分析法求解连续 LTI 系统零状态响应的分析 ...
newmark法 matlab,newmark法程序newmark法程序.doc
newmark法程序newmark法程序用matlab编程法一.法原理 Newmark-?法是一种逐步积分的方法,避免了任何叠加的应用,能很好的适应非线性的反应分析. Newmark-?法假定: ...
Newmark法处理非线性项
Newmark法处理非线性项大家在做考虑非线性因素影响的耦合系统动力响应求解时可能会用到多种方法相结合的思路.有些甚至需要先将非线性项线性化处理.看过很多文章提到由于Newmark法是一步步迭代得到的 ...
牛顿-拉普森法求解线性方程组原理及matlab程序
牛顿-拉普森法求解线性方程组原理及matlab程序牛顿-拉普森法原理 Nowton-Raphson方法matlab程序? 牛顿-拉普森法原理在多变量微积分和矩阵理论的交叉点是求解非线性代数方程 ...
UA PHYS515 电磁理论II 静电场问题5 用Green函数法求解interior Dirichlet问题的例子
UA PHYS515 电磁理论II 静电场问题5 用Green函数法求解interior Dirichlet问题的例子例2 均匀金属空心外壳厚度可忽略的接地球球心位于原点,半径为aaa,用球坐标(r ...
UA PHYS515 电磁理论II 静电场问题4 用Green函数法求解Dirichlet问题
UA PHYS515 电磁理论II 静电场问题4 用Green函数法求解Dirichlet问题上一讲我们讨论过Dirichlet问题的积分解: Φ(r⃗)=∫Vρ(r⃗′)G(r⃗,r⃗′)dx′d ...

NewMark法求解动力学响应

NewMark法求解动力学响应相关推荐

最新文章

热门文章