Newmark方法解运动方程

1. 算法分析

如果将位移按照泰勒公式展开，只保留一阶导数项，    时刻的速度    和位移    可由前一步    时刻的速度和位移表示：

给定初值    和    后，由    时刻的运动方程

可求出    时刻加速度    ，然后由    分别求    时刻的速度    和位移    ，这种方法称为欧拉法，它满足    时刻的运动方程，因此是一种显式格式。欧拉法由前一步的已知值可求下一步的值，故为一步法，可以自起步(self-starting)。但是欧拉法在位移表达式中只保留了    的一阶项，位移的截断误差是    ，因此是一阶精度格式，一般只能用于起步或者与其他方法配合使用。欧拉法可通过在    用导数的平均值代替t时刻的导数值来改进，即

将    带入    得到

   可理解为，在位移和速度表达式中，近似取平均加速度    ,因此该方法也叫平均加速度法。Newmark引入如下的速度和位移关系：

参数    取不同的值可以得到多种方法。如    时即为中心差分法；    时即为平均加速度法。

2.算法流程

Newmark方法的求解过程如下：

(一)初始计算

1. 形成刚度矩阵        ，质量矩阵        和阻尼矩阵        
2. 由初值        和        求解        
3. 由参数        ，时间步长        计算
4. 形成有效刚度矩阵        
5. 分解        

(二) 对每一时间步

1. 计算        时刻的有效载荷
2. 解        时刻位移        
3. 计算        时刻的速度和加速度.

3.算例

用Newmark方法解运动方程，时间步长    

其中

初始条件

   ,常数

有效刚度矩阵

在每个时间步中计算有效载荷

求解方程组

可得到各时刻位移。

4. 编程实现

# Newmark方法
# @表示矩阵乘法
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')
U_0 = np.mat('0; 0')    #初始位移
V_0 = np.mat('0; 0')    #初始速度
A_0 = np.mat('0; 10')   #初始加速度
# 参数
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)
for i in range(step):
    H = c0 * U_0 + c2 * V_0 + c3 * A_0
    Q_1 = Q + M @ H
    U_1 = invKK @ Q_1
    data[0:, [i] ] = U_1
## 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
#解线性方程组采用求逆的方法，计算规模大的，可以用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()