「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多平台发布