共振现象的数值模拟
▲图1
图1所示的单自由度体系，在简谐荷载 作用下的无阻尼受迫振动的运动方程为
方程(1)的解为
式(2)中的前三项都是频率为 的自由振动。其中第一、二两项是由初始条件引起的；第三项与初始条件无关，是伴随激振力的作用而产生的，称为伴生自由振动。第四项则是由激振力所引起并与其频率相同的振动，称为纯受迫振动。当有阻尼存在时前三项所代表的自由振动都将迅速衰减。因此，在实际问题中具有重要性的主要是纯受迫振动。由于它的振幅和频率都是恒定的，因而称为稳态受迫振动。由此可见，单自由度体系在简谐荷载作用下的稳态响应为
式中 即为将动力荷载幅值 作为静力荷载作用于体系时所引起的静位移，而
代表了动位移幅值与静位移之比，称为动力系数，它反映了惯性力的影响，其图像如图2所示
▲图2
当 时, ,即体系的振幅将趋于无穷大。实际结构由于阻尼的存在，振幅不可能趋于无穷大，但它仍将远大于静位移的值，这种现象称为共振。在工程设计中应尽量避免共振现象的发生。
下面用数值模拟来说明这一现象。
例1
▲图3
图3所示等截面两端固定梁，跨中作用有动荷载 。将梁划分为长度相等的四个单元,此时单元长度 ,并将各单元的质量 分为两半集中到单元两端的结点上,如图4所示。各单元局部坐标系的原点均设在单元的左端。
▲图4
采用先处理法求解,结点位移向量为
结构的刚度矩阵和质量矩阵分别为
其一阶频率为 ，且一阶振型是对称分布的。若一阶频率无限接近动荷载的频率，即 ，那么会发生共振现象。令 ，用Newmark方法解振动方程
其中
# Newmark方法模拟梁的共振现象
import numpy as np
import matplotlib.pyplot as plt
l = 1
l2 = l*l
l3 = l2*l
EI = 1
m = 1
t2 = EI/l3
M = 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.28
alpha = 0.5; beta = 0.25
c0 = 1/(beta * dt**2)
c1 = alpha /(beta * dt)
c2 = 1/(beta * dt)
c3 = 1/(2*beta) - 1
c4 = alpha /beta - 1
c5 = dt * (alpha /(2*beta ) - 1)
c6 = dt * (1- alpha )
c7 = alpha * dt
# 有效刚度矩阵
KK = K + c0 * M
invKK = np.linalg.inv(KK)
step = 260 #求解步数
data = np.zeros( step ) # 存储各时间步长下v3的位移
for i in range(step):
H = c0 * U_0 + c2 * V_0 + c3 * A_0
F_1 = np.sin(1.394 * dt*i) * F + M @ H #一阶频率等于动荷载的频率
U_1 = invKK @ F_1
data[i ] = U_1[2,0]
## t+dt时刻的速度和加速度
A_1 = c0 *(U_1 - U_0) - c2*V_0 -c3*A_0
V_1 = V_0 + c6* A_0 + c7* A_1
A_0 = A_1
U_0 = U_1
V_0 = V_1
time = np.zeros(step)
for i in range(step):
t = 0.28 * (i)
time[i] = t
fig, 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方法模拟框架的共振现象
import numpy as np
import matplotlib.pyplot as plt
M = 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.28
alpha = 0.5; beta = 0.25
c0 = 1/(beta * dt**2)
c1 = alpha /(beta * dt)
c2 = 1/(beta * dt)
c3 = 1/(2*beta) - 1
c4 = alpha /beta - 1
c5 = dt * (alpha /(2*beta ) - 1)
c6 = dt * (1- alpha )
c7 = alpha * dt
# 有效刚度矩阵
KK = K + c0 * M
invKK = np.linalg.inv(KK)
step = 450 #求解步数
data2 = np.zeros( step ) # 存储各时间步长下第二层的水平位移
data3 = np.zeros( step ) # 存储各时间步长下第3层的水平位移
for i in range(step):
H = c0 * U_0 + c2 * V_0 + c3 * A_0
F_1 = np.sin(0.386 * dt*i) * F + M @ H #一阶频率等于动荷载的频率
U_1 = invKK @ F_1
data2[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_0
V_1 = V_0 + c6* A_0 + c7* A_1
A_0 = A_1
U_0 = U_1
V_0 = V_1
time = np.zeros(step)
for i in range(step):
t = 0.28 * (i)
time[i] = t
fig, (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
