Newmark方法解运动方程

fem178

发布于 2023-01-30 15:10:58

4750

发布于 2023-01-30 15:10:58

文章被收录于专栏：数值分析与有限元编程数值分析与有限元编程

1. 算法分析

如果将位移按照泰勒公式展开，只保留一阶导数项，

t+\Delta t

时刻的速度

\mathbf {\dot a}_{t+\Delta t}

和位移

\mathbf {a}_{t+\Delta t}

可由前一步

时刻的速度和位移表示：

\mathbf {a}_{t+\Delta t} = \mathbf {a}_t + \mathbf {\dot a}_t{\Delta t}\cdots(1)

\mathbf {\dot a}_{t+\Delta t} = \mathbf {\dot a}_t + \mathbf {\ddot a}_t{\Delta t} \cdots(2)

给定初值

\mathbf {a}_0

和

\mathbf {\dot a}_0

后，由

时刻的运动方程

\mathbf M \mathbf {\ddot a}_t + \mathbf C \mathbf {\dot a}_t + \mathbf K \mathbf a_t = \mathbf Q_t\cdots(3)

可求出

时刻加速度

\mathbf {\ddot a}_t

，然后由

(1)(2)

分别求

t+\Delta t

时刻的速度

\mathbf {\dot a}_{t+\Delta t}

和位移

\mathbf {a}_{t+\Delta t}

，这种方法称为欧拉法，它满足

时刻的运动方程，因此是一种显式格式。欧拉法由前一步的已知值可求下一步的值，故为一步法，可以自起步(self-starting)。但是欧拉法在位移表达式中只保留了

{\Delta t}

的一阶项，位移的截断误差是

O({\Delta t}^2)

，因此是一阶精度格式，一般只能用于起步或者与其他方法配合使用。欧拉法可通过在

(1)(2)

用导数的平均值代替t时刻的导数值来改进，即

\mathbf {a}_{t+\Delta t} = \mathbf {a}_t + \frac{1}{2}(\mathbf {\dot a}_t+\mathbf {\dot a}_{t+\Delta t}){\Delta t}\cdots(4)

\mathbf {\dot a}_{t+\Delta t} = \mathbf {\dot a}_t + \frac{1}{2}(\mathbf {\ddot a}_t+\mathbf {\ddot a}_{t+\Delta t}){\Delta t}\cdots(5)

将

(5)

带入

(4)

得到

\mathbf {a}_{t+\Delta t} = \mathbf {a}_t + \mathbf {\dot a}_t{\Delta t} + \frac{1}{4}(\mathbf {\ddot a}_t+\mathbf {\ddot a}_{t+\Delta t}){\Delta t}^2 \cdots(6)

(5)(6)

可理解为，在位移和速度表达式中，近似取平均加速度

\frac{1}{2}(\mathbf {\ddot a}_t+\mathbf {\ddot a}_{t+\Delta t})

,因此该方法也叫平均加速度法。 Newmark引入如下的速度和位移关系：

\mathbf {\dot a}_{t+\Delta t} = \mathbf {\dot a}_t + (1-\alpha)\mathbf {\ddot a}_t{\Delta t}+ \alpha \mathbf {\ddot a}_{t+\Delta t}{\Delta t}\cdots(7)

\mathbf {a}_{t+\Delta t} = \mathbf {a}_t + \mathbf {\dot a}_t{\Delta t} + (\frac{1}{2} - \beta)\mathbf {\ddot a}_t {\Delta t}^2 + \beta \mathbf {\ddot a}_{t+\Delta t}{\Delta t}^2 \cdots(8)

参数

\alpha, \beta

取不同的值可以得到多种方法。如

\alpha = \frac{1}{2}，\beta =0

时即为中心差分法；

\alpha = \frac{1}{2}，\beta =\frac{1}{4}

时即为平均加速度法。

2.算法流程

Newmark方法的求解过程如下：

(一)初始计算

形成刚度矩阵

\mathbf K

，质量矩阵

\mathbf M

和阻尼矩阵

\mathbf C

由初值

\mathbf a_0

和

\mathbf {\dot a}_0

求解

\mathbf {\ddot a}_0

由参数

\alpha, \beta

，时间步长

\Delta t

计算

c_0=\frac{1}{\beta {\Delta t}^2},c_1=\frac{\alpha}{\beta \Delta t},c_2= \frac{1}{\beta \Delta t},c_3=\frac{1}{2\beta}-1,c_4=\frac{\alpha}{\beta}-1, c_5=\Delta t(\frac{\alpha}{2\beta}-1),c_6=\Delta t(1-\alpha),c_7=\Delta t \alpha

形成有效刚度矩阵

\overline {\mathbf K}= K + c_0\mathbf M+c_1\mathbf C

分解

\overline {\mathbf K}= \mathbf {LD}\mathbf L^T

(二) 对每一时间步

计算

t+\Delta t

时刻的有效载荷

\overline {\mathbf Q}_{t+\Delta t} = \mathbf Q_{t+\Delta t} +\mathbf M(c_0 \mathbf {a}_t+ c_2\mathbf {\dot a}_t+ c_3\mathbf {\ddot a}_t) + \mathbf C(c_1 \mathbf {a}_t+ c_4\mathbf {\dot a}_t+ c_5\mathbf {\ddot a}_t)

t+\Delta t

时刻位移

(\mathbf L\mathbf D\mathbf L^T )\mathbf a_{t+\Delta t} = \overline {\mathbf Q}_{t+\Delta t}

计算

t+\Delta t

时刻的速度和加速度.

\mathbf {\ddot a}_{t+\Delta t} = c_0(\mathbf {a}_{t+\Delta t}-\mathbf {a}_{t}) - c_2\mathbf {\dot a}_{t} - c_3\mathbf {\ddot a}_{t}

\mathbf {\dot a}_{t+\Delta t} = \mathbf {\dot a}_{t}+ c_6 \mathbf {\ddot a}_{t} + c_7 \mathbf {\ddot a}_{t+\Delta t}

3.算例

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

\Delta t =0.28

\mathbf M \mathbf {\ddot a}(t) + \mathbf K \mathbf a(t) = \mathbf Q

其中

\mathbf M= \begin{bmatrix} 2 & 0 \\ 0 & 1 \\ \end{bmatrix} , \mathbf K= \begin{bmatrix} 6 & -2 \\ -2 & 4 \\ \end{bmatrix} , \mathbf Q= \begin{bmatrix} 0 \\ 10 \end{bmatrix}

初始条件

\mathbf {a}(0) =\begin{bmatrix} 0 \\ 0 \end{bmatrix} ,\mathbf {\dot a}(0) = \begin{bmatrix} 0 \\ 10 \end{bmatrix}

\alpha=0.5, \beta=0.25

,常数

c_0=51.0，c_1=7.14，c_2=14.3，c_3=1，c_4=1，c_5=0，c_6=0.14，c_7=0.14

有效刚度矩阵

\overline {\mathbf K} = \begin{bmatrix} 108 & -2\\ -2 & 55 \\ \end{bmatrix}

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

\overline {\mathbf Q}_{t+\Delta t} = \begin{bmatrix} 0 \\ 10 \\ \end{bmatrix} + \begin{bmatrix} 2 & 0\\ 0 & 1 \\ \end{bmatrix} (51\mathbf {a}_t + 14.3 \mathbf {\dot a}_t + 1.0 \mathbf {\ddot a}_t)

求解方程组

\overline {\mathbf K}\mathbf a_{t+\Delta t} = \overline {\mathbf Q}_{t+\Delta t}

可得到各时刻位移。

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()

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2023-01-03，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自数值分析与有限元编程微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度