多级弹簧-质量系统瞬态分析（基于Newmark）

 该程序主要用于实现多自由度下动态系统的隐式分析，计算模型参见多级弹簧-质量系统的瞬态计算（基于中心差分法）。下面是具体的Python程序。

#coding:utf-8
"""
Created on Tue Nov 20 12:41:40 2018

@author: yujin.wang
"""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys
sys.setrecursionlimit(10000) # 10000 is an example, try with different values
################Force###############
def CenterDiff(t,endtime,m,k,fext,u,delta,du0):
	fint = k.dot(u)
	ddu = m.I.dot((fext - fint))
	du = du0 + ddu * delta
	u = u + du*delta
	t += delta
	t_list.append(t);u1_list.append(u[0,0]);u2_list.append(u[1,0])
	if t > endtime:
		return t_list,u1_list,u2_list
	else:
		return CenterDiff(t,endtime,m,k,fext,u,delta,du)
	
def NewmarkSolver(t,endT,deltaT,TOL,u,du,ddu,k,c,m,iternum):
	t += deltaT
	a1 = m/(beta*deltaT*deltaT) + gamma *c/(beta*deltaT)
	a2 = m/(beta*deltaT) + c*(gamma/beta-1)
	a3 = m*(1./(2*beta)-1) + deltaT*c*(gamma/(2*beta)-1)

	k_ = k + a1

	f_ = f + a1*u + a2*du + a3*ddu
	un = k_.I*f_
	Rn = f_ - k*un - a1*un
#	print Rn
	print '##############TIME: %f####################' %(t)
	print time.ctime()
	print 'Iter.:%d Rn:%f' %(iternum,abs(Rn).max())
#Newton-Raphson
	while abs(Rn).max() > TOL:
		iternum += 1
		k_= k_ + a1
		dun = k_.I*Rn
		un = un + dun
		Rn = f_ - k*un - a1*un
		print 'Iter.:%d Rn:%f' %(iternum,abs(Rn).max())
	dun = gamma*(un-u)/(beta*deltaT) + du*(1-gamma/beta) + deltaT*ddu*(1-gamma/(2*beta))
	ddun = (un-u)/(beta*deltaT*deltaT) - du/(beta*deltaT) - ddu*(1./(2*beta)-1)
	t_list.append(t)
	u1_list.append(un[0,0])
	u2_list.append(un[1,0])
	if t+deltaT > endT:
		return  t_list,u1_list,u2_list
	else:
		print t,un
		return NewmarkSolver(t,endT,deltaT,TOL,un,dun,ddun,k,c,m,iternum)
	
if __name__ == "__main__":
	iternum = 0
	gamma=0.5;
	beta=0.25;
	TOL =1e-3
	deltaT = 0.02
	m1=1.;m2=1.
	k1 = 100.;k2=100.;k3=0.001
	f1=0.;f2=1.
	m = np.matrix([[m1,0],[0,m2]])
	c = np.matrix([[0,0],[0,0]])
	k = np.matrix([[k1+k2,-k2],[-k2,k2+k3]])
	f = np.matrix([[f1],[f2]])

	u0 = np.matrix([[0],[0]])
	du0 = np.matrix([[0],[0]])
	u1_list = [];u2_list = []
	t_list = []
	un_list = []

	ddu0 = np.linalg.inv(m)*(f - c*du0 )
	t_list,u1_list,u2_list = NewmarkSolver(0,25.,deltaT,TOL,u0,du0,ddu0,k,c,m,iternum)
	plt.figure(1)
	plt.plot(t_list,u1_list,'red')
	plt.plot(t_list,u2_list,'blue')
	plt.legend(['X1','X2'])
	plt.grid()
	plt.show()

设置阻尼矩阵C=[[1,0],[0,1]]