关于python:Python与地震工程单自由度体系求解之Newmarkβ法

「Python与地震工程」单自由度体系求解之Newmark-β法

原理

Newmark-β法是地震工程畛域最经典的逐渐积分算法。

推导过程请查阅构造动力学或地震工程学教材，此处仅简略列出逐渐递推公式。

已知第 \( i \) 步响应，则第 \( i+1 \) 步位移响应可按下式计算：

$$
u_{i+1}=\frac{\hat{p}_{i+1}}{\hat{k}}
$$

其中

$$
\hat{k}=k+\frac{\gamma}{\beta \Delta t}c+\frac{1}{\beta \left( \Delta t \right) ^2}m
$$

$$
\hat{p}_{i+1}=p_{i+1}+\left[ \frac{1}{\beta \left( \Delta t \right) ^2}m+\frac{\gamma}{\beta \Delta t}c \right] u_i+\left[ \frac{1}{\beta \Delta t}m+\left( \frac{\gamma}{\beta}-1 \right) c \right] \dot{u}_i
\\
+\left[ \left( \frac{1}{2\beta}-1 \right) m+\Delta t\left( \frac{\gamma}{2\beta}-1 \right) c \right] \ddot{u}_i
$$

则第 i+1 步速度、加速度响应为

$$
\dot{u}_{i+1}=\frac{\gamma}{\beta \Delta t}\left( u_{i+1}-u_i \right) +\left( 1-\frac{\gamma}{\beta} \right) \dot{u}_i+\Delta t\left( 1-\frac{\gamma}{2\beta} \right) \ddot{u}_i
$$

$$
\ddot u_{i + 1} = \frac{{{p_{i + 1}} – c{{\dot u}_{i + 1}} – k{u_{i + 1}}}}{m}
$$

对于地震激励

$$p_{i+1} = -ma_{g,i+1}$$

程序代码

import numpy as np
import matplotlib.pyplot as plt

plt.style.use("ggplot")
plt.rcParams['font.sans-serif'] = ['Microsoft Yahei']
plt.rcParams['axes.unicode_minus'] = False

def draw_response(title, ta, a, t, u):
    plt.figure(title,(12,6))
    plt.subplot(2,1,1)
    plt.plot(ta,a,label=r"输出地震波加速度时程")
    plt.grid(True)
    plt.legend()
    plt.xlim(0,t[-1])
    plt.subplot(2,1,2)
    plt.plot(t,u,label=r"SDOF体系位移响应时程")
    plt.xlabel(r"工夫 (s)")
    plt.grid(True)
    plt.legend()
    plt.xlim(0,t[-1])
    plt.show()


def solve_sdof_eqwave_nmk(omg, zeta, ag, dt):

    n = len(ag)
    omg2 = omg*omg
    
    gma = 0.5
    bta = 0.25

    u0 = 0.0
    v0 = 0.0

    c1 = 1.0/bta/dt/dt
    c2 = 1.0/bta/dt
    c3 = gma/bta/dt
    c4 = 1.0 - gma/bta
    c5 = 1.0 - 0.5*gma/bta
    c6 = 0.5/bta - 1.0

    c = 2.0*zeta*omg
    c7 = c1 + c3*c
    c8 = c2 - c4*c
    c9 = c6 - dt*c5*c

    u = np.zeros(n)
    v = np.zeros(n)
    a = np.zeros(n)

    kh = omg2 + c7
    u[0] = u0
    v[0] = v0
    a[0] = -ag[0]-c*v[0]-omg2*u[0]

    for i in range(n-1):
        ph = -ag[i+1] + c7*u[i] + c8*v[i] + c9*a[i]
        u[i+1] = ph/kh
        v[i+1] = c3*(u[i+1]-u[i]) + c4*v[i] + dt*c5*a[i]
        a[i+1] = -ag[i+1]-c*v[i+1]-omg2*u[i+1]
    
    return u, v


if __name__ == '__main__':
    
    acc0 = np.loadtxt("EQ-S-3.txt") # 读取地震波
    dt = 0.005 # 工夫距离
    n = len(acc0)
    t0 = np.linspace(0.0,dt*(n-1),n)

    # 对时程数据进行补零以显示地震完结后一段时间内的自在振动衰减状况
    ne = round(n*1.2)
    t = np.linspace(0.0,dt*(ne-1),ne)
    ag = np.zeros(ne)
    ag[:n] = acc0

    omg = 2.0*np.pi # 自振圆频率
    zeta = 0.05 # 固有阻尼比

    u, v = solve_sdof_eqwave_nmk(omg, zeta, ag, dt)

    draw_response("Seismic Response -- nmk", t0, acc0, t, u) # 绘图

文中所用地震波下载：
EQ-S-3.txt

后果

转载本文请注明出处。科研成果中应用本文代码请援用作者课题组的钻研论文（请戳“浏览原文”）。

本文由mdnice多平台公布