「Python与地震工程」单自由度体系求解之通用微分方程数值解法（基于scipy.integrate.odeint）

「Python与地震工程」单自由度体系求解之通用微分方程数值解法（基于scipy.integrate.odeint）

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

原理

单自由度体系分析是结构动力学和地震工程学习的启蒙知识，也是复杂结构动力分析的基石，具有重要的作用和研究价值。本文将探讨如何用Python的科学计算模块来实现单自由度体系的动力响应分析。

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

$$
mddot{x}left( t right) + cdot{x}left( t right) + kxleft( t right) = Pleft( t right)
$$

或

$$
ddot{x}left( t right) + 2 zeta omega_0 dot{x}left( t right) + omega_0^2 xleft( t right) = pleft( 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}}\ {pleft( t right) – 2zeta 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()

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

结果

转载本文请注明出处。

「Python与地震工程」单自由度体系求解之通用微分方程数值解法

科研成果中使用本文代码请引用作者课题组的研究论文：

作者课题组发表的研究论文列表（持续更新中……）

服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net