▲图1图1所示的单自由度体系，在简谐荷载作用下的无阻尼受迫振动的运动方程为方程(1)的解为式(2)中的前三项都是频率为的自由振动。其中第一、二两项是由初始条件引起的；第三项与初始条件无关，是伴随激振力的作用而产生的，称为伴生自由振动。第四项则是由激振力所引起并与其频率相同的振动，称为纯受迫振动。当有阻尼存在时前三项所代表的自由振动都将迅速衰减。因此，在实际问题中具有重要性的主要是纯受迫振动。由于它的振幅和频率都是恒定的，因而称为稳态受迫振动。由此可见，单自由度体系在简谐荷载作用下的稳态响应为式中即为将动力荷载幅值作为静力荷载作用于体系时所引起的静位移，而代表了动位移幅值与静位移之比，称为动力系数，它反映了惯性力的影响，其图像如图2所示▲图2当时,,即体系的振幅将趋于无穷大。实际结构由于阻尼的存在，振幅不可能趋于无穷大，但它仍将远大于静位移的值，这种现象称为共振。在工程设计中应尽量避免共振现象的发生。下面用数值模拟来说明这一现象。例1▲图3图3所示等截面两端固定梁，跨中作用有动荷载。将梁划分为长度相等的四个单元,此时单元长度,并将各单元的质量分为两半集中到单元两端的结点上,如图4所示。各单元局部坐标系的原点均设在单元的左端。▲图4采用先处理法求解,结点位移向量为结构的刚度矩阵和质量矩阵分别为其一阶频率为，且一阶振型是对称分布的。若一阶频率无限接近动荷载的频率，即，那么会发生共振现象。令，用Newmark方法解振动方程其中#Newmark方法模拟梁的共振现象importnumpyasnpimportmatplotlib.pyplotaspltl=1l2=l*ll3=l2*lEI=1m=1t2=EI/l3M=np.array([[m,0,0,0,0,0],[0,0,0,0,0,0],[0,0,m,0,0,0],[0,0,0,0,0,0],[0,0,0,0,m,0],[0,0,0,0,0,0]])K=t2*np.array([[24,0,-12,6*l,0,0],[0,8*l2,-6*l,2*l2,0,0],[-12,-6*l,24,0,-12,6*l],[6*l,2*l2,0,8*l2,-6*l,2*l2],[0,0,-12,-6*l,24,0],[0,0,6*l,2*l2,0,8*l2]])F=np.array([[0],[0],[2],[0],[0],[0]])U_0=np.array([[0],[0],[0],[0],[0],[0]])#初始位移V_0=np.array([[0],[0],[0],[0],[0],[0]])#初始速度A_0=np.array([[0],[0],[0],[0],[0],[0]])#初始加速度#参数dt=0.28alpha=0.5;beta=0.25c0=1/(beta*dt**2)c1=alpha/(beta*dt)c2=1/(beta*dt)c3=1/(2*beta)-1c4=alpha/beta-1c5=dt*(alpha/(2*beta)-1)c6=dt*(1-alpha)c7=alpha*dt#有效刚度矩阵KK=K+c0*MinvKK=np.linalg.inv(KK)step=260#求解步数data=np.zeros(step)#存储各时间步长下v3的位移foriinrange(step):H=c0*U_0+c2*V_0+c3*A_0F_1=np.sin(1.394*dt*i)*F+M@H#一阶频率等于动荷载的频率U_1=invKK@F_1data[i]=U_1[2,0]##t+dt时刻的速度和加速度A_1=c0*(U_1-U_0)-c2*V_0-c3*A_0V_1=V_0+c6*A_0+c7*A_1A_0=A_1U_0=U_1V_0=V_1time=np.zeros(step)foriinrange(step):t=0.28*(i)time[i]=tfig,ax=plt.subplots(1,1,figsize=(5,12))ax.plot(time,data,'r-')ax.set_xlabel('$t/s$',fontsize=14)ax.set_ylabel('$v3/m$',fontsize=14)fig.savefig('f459.png',dpi=400)plt.show()▲图5例2▲图6图6所示三层刚架，各横梁为无限刚性，刚架的质量全部集中在横梁上，分别为，各层间侧移刚度分别为，第二层横梁上作用有水平简谐荷载。各层横梁分别发生单位侧移时体系的刚度系数分别为：结构的刚度矩阵和质量矩阵分别为于是，可求得自振频率对应的主振型如下▲图7结构的受迫振动方程为其中若一阶频率无限接近动荷载的频率，即，那么会发生共振现象。令，用Newmark方法求出内第二层的水平位移。#Newmark方法模拟框架的共振现象importnumpyasnpimportmatplotlib.pyplotaspltM=np.array([[1.5,0,0],[0,1,0],[0,0,1.5]])K=np.array([[2,-1,0],[-1,2,-1],[0,-1,1]])F=np.array([[0],[2],[0]])U_0=np.array([[0],[0],[0]])#初始位移V_0=np.array([[0],[0],[0]])#初始速度A_0=np.array([[0],[0],[0]])#初始加速度#参数dt=0.28alpha=0.5;beta=0.25c0=1/(beta*dt**2)c1=alpha/(beta*dt)c2=1/(beta*dt)c3=1/(2*beta)-1c4=alpha/beta-1c5=dt*(alpha/(2*beta)-1)c6=dt*(1-alpha)c7=alpha*dt#有效刚度矩阵KK=K+c0*MinvKK=np.linalg.inv(KK)step=450#求解步数data2=np.zeros(step)#存储各时间步长下第二层的水平位移data3=np.zeros(step)#存储各时间步长下第3层的水平位移foriinrange(step):H=c0*U_0+c2*V_0+c3*A_0F_1=np.sin(0.386*dt*i)*F+M@H#一阶频率等于动荷载的频率U_1=invKK@F_1data2[i]=U_1[1,0]data3[i]=U_1[2,0]##t+dt时刻的速度和加速度A_1=c0*(U_1-U_0)-c2*V_0-c3*A_0V_1=V_0+c6*A_0+c7*A_1A_0=A_1U_0=U_1V_0=V_1time=np.zeros(step)foriinrange(step):t=0.28*(i)time[i]=tfig,(ax1,ax2)=plt.subplots(2,1,figsize=(6,12))ax1.plot(time,data2,'r-')ax1.set_xlabel('$t/s$',fontsize=14)ax1.set_ylabel('$u_2/m$',fontsize=14)ax2.plot(time,data3,'r-')ax2.set_xlabel('$t/s$',fontsize=14)ax2.set_ylabel('$u_3/m$',fontsize=14)fig.savefig('f439.png',dpi=400)plt.show()▲图8来源：数值分析与有限元编程
