关于python:Python与地震工程单自由度体系求解之通用微分方程数值解法基于scipyintegrateodeint

0次阅读

共计 2299 个字符，预计需要花费 6 分钟才能阅读完成。

背景图：2008 年 5 月汶川地震

单自由度体系剖析是构造动力学和地震工程学习的启蒙常识，也是简单构造能源剖析的基石，具备重要的作用和钻研价值。本文将探讨如何用 Python 的科学计算模块来实现单自由度体系的能源响应分析。

经典的弹性有阻尼单自由度体系在任意激励 P(t) 作用下的静止方程可写为

$$
m\ddot{x}\left(t \right) + c\dot{x}\left(t \right) + kx\left(t \right) = P\left(t \right)
$$

或

$$
\ddot{x}\left(t \right) + 2 \zeta \omega_0 \dot{x}\left(t \right) + \omega_0^2 x\left(t \right) = p\left(t \right)
$$

求解此静止方程的经典办法有杜哈梅积分法、Newmark- β 法、核心差分法等。借助 Python 求解此方程最快捷（编程者的工作量小，并非最疾速）的办法就是将方程降阶为状态空间模式调用 scipy.integrate.odeint() 函数。须要留神的是，对于地震激励，须要把离散的地震波“包装”成一个连续函数能力被 odeint 辨认。

上述方程的状态空间模式为

$$
\frac{{\rm{d}}}{{{\rm{d}}t}}\left({\begin{array}{c} {{y_1}}\\ {{y_2}} \end{array}} \right) = \left({\begin{array}{c} {{y_2}}\\ {p\left( t \right) – 2\zeta \omega_0 {y_2} – {\omega_0^2}{y_1}} \end{array}} \right)
$$

采纳 odeint() 进行单自由度体系共振简谐激励和地震激励作用下响应求解的残缺代码为

import scipy as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def pd(tn,*p_args):
    """等工夫距离离散信号 a 对应的间断化函数（线性插值办法实现）"""
    dt = p_args[0] # 离散信号工夫距离
    a  = p_args[1] # 离散信号数据

    ind = int(sp.floor(tn/dt))
    if (ind+1)>=len(a):
        return 0.0
    else:
        al = a[ind]; ar = a[ind + 1];
        k = (ar - al) / dt
        return al+k*(tn-ind*dt)

def solve_sdof_resonance(omg = 1.0*2.0*sp.pi, zeta = 0.02):

    y0 = sp.asarray([0.0,0.0]) # 初始条件
    p = lambda t: sp.sin(omg*t) # 共振简谐激励

    f_sdof = lambda y,t: sp.asarray([y[1], -p(t)-2.0*zeta*omg*y[1]-omg*omg*y[0] ])

    t = sp.linspace(0, 30, 1001)
    y = odeint(f_sdof, y0, t)

    # 绘图:
    plt.figure("Resonance Response",(12,4))
    plt.plot(t,y[:,0])
    plt.grid(True)
    plt.show()


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_odeint(omg0 = 1.0*2.0*sp.pi, zeta = 0.05):

    y0 = sp.asarray([0.0,0.0]) # 初始条件

    acc0 = sp.loadtxt("EQ-S-3.txt") # 读取地震波
    dt = 0.005 # 工夫距离
    n = len(acc0)
    t0 = sp.linspace(0.0, dt*(n-1), n)
    p = lambda t: pd(t, dt, acc0) # 线性间断化地震激励

    f_sdof = lambda y,t: sp.asarray([y[1], -p(t)-2.0*zeta*omg0*y[1]-omg0*omg0*y[0] ])

    # 显示地震完结后一段时间内的自在振动衰减状况
    ne = round(n*1.2)
    t = sp.linspace(0.0,dt*(ne-1),ne)
    y = odeint(f_sdof, y0, t)

    draw_response("Seismic Response -- scipy.integrate.odeint", t0, acc0, t, y[:,0]) # 绘图


if __name__ == '__main__':
    solve_sdof_resonance()
    solve_sdof_eqwave()