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