同样可得向后差分公式:
以上两式相减和相加分别得到:
以上两式忽略高阶小量,可得到 时刻速度和加速度表达式:
为了求解 时刻的位移,将 代入 时刻动力学方程
得到
其中
若已经求得
在向后差分公式
其中
中心差分法解动力学方程的算法可归纳为
(一)初始计算
形成刚度矩阵 ,质量矩阵 和阻尼矩阵 由初值 和 求解 和 由时间步长 计算 计算 计算有效质量矩阵 对 进行分解
(二) 对每一时间步
计算 时刻的有效载荷 解 时刻位移 如果需要,按照 计算 时刻的速度和加速度
其中
初始条件
将初始条件代入方程,解得
每个时间步长计算
以及
# 中心差分法
# @表示矩阵乘法
import numpy as np
import matplotlib.pyplot as plt
step = 10#求解步数
data = np.zeros( (2,step) ) # 存储各时间步长的数据
K = np.mat('6,-2;-2,4')
M = np.mat('2,0;0,1')
Q = np.mat('0;10')
a_0 = np.zeros((2,1))
velocity_0 = np.zeros((2,1))
# 求初始加速度
acceleration_0 = np.linalg.inv(M) *(Q - K @ a_0)
dt = 0.28
c0 = 1/(dt**2)
c1 = 0.5 * dt
c2 = 2 * c0
c3 = 1 / c2
a_dt = a_0 - dt * velocity_0 + c3 * acceleration_0
# 有效质量矩阵
MM = c0 * M
invMM = np.linalg.inv(MM)
TMP1 = K - c2*M
A_t_sub_dt = a_dt
A_t = a_0
A_t_plus_dt = np.zeros((2,1))
for i in range(step):
TMP2 = M @ A_t_sub_dt
QQ = Q - TMP1 @ A_t - c0 * TMP2
A_t_plus_dt = invMM @ QQ
data[0:,[i] ] = A_t_plus_dt
A_t_sub_dt = A_t
A_t = A_t_plus_dt
#解线性方程组采用求逆的方法,计算规模大的,可以用LDLT分解.
#
time = np.zeros((step))
for i in range(step):
t = 0.28 * (i+1)
time[i] = t
labels =['Δt','2Δt','3Δt','4Δt','5Δt','6Δt','7Δt','8Δt','9Δt','10Δt']
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5,12) )
ax1.plot(time, data[0, :], 'r-')
ax1.set_xticks(time, labels)
ax1.set_ylabel('$a_1$',fontsize = 14)
ax2.plot(time, data[1, :], 'b-')
ax2.set_xticks(time, labels)
ax2.set_ylabel('$a_2$',fontsize = 14)
fig.savefig('./f429.png', dpi = 400)
plt.show()