Newmark方法解运动方程
1. 算法分析如果将位移按照泰勒公式展开，只保留一阶导数项， 时刻的速度 和位移 可由前一步 时刻的速度和位移表示： 给定初值 和 后，由 时刻的运动方程 可求出 时刻加速度 ，然后由 分别求 时刻的速度 和位移 ，这种方法称为欧拉法，它满足 时刻的运动方程，因此是一种显式格式。欧拉法由前一步的已知值可求下一步的值，故为一步法，可以自起步(self-starting)。但是欧拉法在位移表达式中只保留了 的一阶项，位移的截断误差是 ，因此是一阶精度格式，一般只能用于起步或者与其他方法配合使用。欧拉法可通过在 用导数的平均值代替t时刻的导数值来改进，即 将 带入 得到 可理解为，在位移和速度表达式中，近似取平均加速度 ,因此该方法也叫平均加速度法。Newmark引入如下的速度和位移关系： 参数 取不同的值可以得到多种方法。如 时即为中心差分法； 时即为平均加速度法。2.算法流程Newmark方法的求解过程如下：(一)初始计算形成刚度矩阵 ，质量矩阵 和阻尼矩阵 由初值 和 求解 由参数 ，时间步长 计算形成有效刚度矩阵 分解 (二) 对每一时间步计算 时刻的有效载荷解 时刻位移 计算 时刻的速度和加速度. 3.算例用Newmark方法解运动方程，时间步长 其中 初始条件 ,常数有效刚度矩阵 在每个时间步中计算有效载荷 求解方程组 可得到各时刻位移。4. 编程实现# Newmark方法# @表示矩阵乘法import numpy as npimport matplotlib.pyplot as pltstep = 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.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 * 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] = tlabels =['Δ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()来源：数值分析与有限元编程
