《Python与地震工程》Newmark-β法求解单自由度系统的原理Newmark-β法是地震工程领域最经典的逐步积分算法。推导过程可参考结构动力学或地震工程教材。这里只简单列出一步一步的递归公式。给定第\(i\)步的响应，第\(i+1\)步的位移响应可以计算如下：$$u_{i+1}=\frac{\hat{p}_{i+1}}{\hat{k}}$$其中$$\hat{k}=k+\frac{\gamma}{\beta\Deltat}c+\frac{1}{\beta\left(\Deltat\right)^2}m$$$$\hat{p}_{i+1}=p_{i+1}+\left[\frac{1}{\beta\left(\Deltat\right)^2}m+\frac{\gamma}{\beta\Deltat}c\right]u_i+\left[\frac{1}{\beta\Deltat}m+\left(\frac{\gamma}{\beta}-1\right)c\right]\dot{u}_i\\+\left[\left(\frac{1}{2\beta}-1\right)m+\Deltat\left(\frac{\gamma}{2\beta}-1\right)c\right]\ddot{u}_i$$那么第i+1步的速度和加速度响应为$$\dot{u}_{i+1}=\frac{\gamma}{\beta\Deltat}\left(u_{i+1}-u_i\right)+\left(1-\frac{\gamma}{\beta}\right)\dot{u}_i+\Deltat\left(1-\frac{\gamma}{2\beta}\right)\ddot{u}_i$$$$\ddotu_{i+1}=\frac{{{p_{i+1}}-c{{\dotu}_{i+1}}-k{u_{i+1}}}}{m}$$对于地震激发$$p_{i+1}=-ma_{g,i+1}$$程序代码importnumpyasnpiimportmatplotlib.pyplotaspltplt.style.use("ggplot")plt.rcParams['font.sans-serif']=['微软雅黑']plt.rcParams['axes.unicode_minus']=Falseheadefdraw_response(title,ta,a,t,u):plt.figure(title,(12,6))plt.subplot(2,1,1)plt.plot(ta,a,label=r"Part.grid(True)plt.legend()plt.xlim(0,t[-1])plt.subplot(2,1,2)plt.plot(t,2)plt.grid(True)plt.legend()u,label=r"SDOFSpecificEnhancementLabel")plt.xlabel(r"Speed(s)")plt.grid(True)plt.legend()plt.xlim(0,t[-1])plt.show()defsolve_sdof_eqwave_nmk(omg,zeta,ag,dt):n=len(ag)omg2=omg*omggma=0.5bta=0.25u0=0.0v0=0.0c1=1.0/bta/dt/dtc2=1.0/bta/dtc3=gma/bta/dtc4=1.0-gma/btac5=1.0-0.5*gma/btac6=0.5/bta-1.0c=2.0*zeta*omgc7=c1+c3*cc8=c2-c4*cc9=c6-dt*c5*cu=np.zeros(n)v=np.zeros(n)a=np.zeros(n)kh=omg2+c7u[0]=u0v[0]=v0a[0]=-ag[0]-c*v[0]-omg2*u[0]foriinrange(n-1):ph=-ag[i+1]+c7*u[i]+c8*v[i]+c9*a[i]u[i+1]=ph/khv[i+1]=c3*(u[i+1]-u[i])+c4*v[i]+dt*c5*a[i]a[i+1]=-ag[i+1]-c*v[i+1]-omg2*u[i+1]返回u,vif__name__=='__main__':acc0=np.loadtxt("EQ-S-3.txt")#读取地震波dt=0.005#时间间隔n=len(acc0)t0=np.linspace(0.0,dt*(n-1),n)#时间历史数据填零以显示时间段之后地震自由振动的衰减ne=round(n*1.2)t=np.linspace(0.0,dt*(ne-1),ne)ag=np.zeros(ne)ag[:n]=acc0omg=2.0*np.pi#自然振动频率zeta=0.05#固有阻尼比u,v=solve_sdof_eqwave_nmk(omg,zeta,ag,dt)draw_response("SeismicResponse--nmk",t0,acc0,t,u)#绘图文本Seismicwave下载使用：EQ-S-3.txt本文成果转载请注明科研成果出处，转载请注明作者课题组研究论文（请点击“阅读原文”）。本文由mdnice多平台发布
