精细积分法的算例

对于一个单自由度的体系，质量  ，刚度  ，阻尼比为  。初始状态下 体系的处于静止状态，当  时刻作用一个大小为  的简谐荷载，作用时间为  ，求体系的位移时程，并将结果与采用  解析积分法作比较。

 积分的基本原理是将任意载荷视为一系列脉冲激励的叠加，系统的响应可以表示成对这些脉冲激励的响应函数的叠加。在单自由度系统中，  积分的形式为：

其中，  是系统在时间  的位移响应，  是外载荷函数，   是质量，   是无阻尼自然频率，  为有阻尼时结构的频率，  是阻尼比‌.

  
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import numpy as np


defcompute_t(n, H, eta):
    ## 计算矩阵T
    N = 20
    tau = eta/(2**N)
    H_tau1 = H*tau
    H_tau2 = H @ H_tau1 *tau
    H_tau3 = H @ H_tau2 *tau
    H_tau4 = H @ H_tau3 *tau
    T_a = H_tau1 + H_tau2/2.0 + H_tau3/6.0 + H_tau4/24.0
    for i inrange(N):
        T_a = 2.0*T_a + T_a @ T_a

    T = np.eye(n) + T_a  

    return T 

defprecise_integration_solver( H,F,ug,step,dt):
    # 精细积分法
    invH = np.linalg.inv(H)
    T = compute_t(2, H, dt)

    x  = np.zeros((step, 2))
    xi = np.zeros(2)

    for i inrange(step-1):
        r0 = -F*(-ug[i])
        r1 = -F*(-ug[i+1]+ug[i]) / dt

        xi = T.dot(xi+invH.dot(r0+invH.dot(r1))) - invH.dot(r0+invH.dot(r1)+dt*r1)
        x[i,:] = xi

    u = x[:,0]
    return u


# 杜哈梅积分函数
deffunc(tau):   
    return np.sin(tau) * np.exp( -zeta*omg*(t-tau) ) * np.sin(omg*(t- tau ) )

if __name__ == "__main__":
    plt.style.use("ggplot")
    plt.rcParams['font.sans-serif'] = ['MiSans']    # 这个字体能够支持中英文字符
    plt.rcParams['axes.unicode_minus'] = False      # 正常显示负号

    omg = 1.0          # 自振频率
    zeta = 0.05          # 固有阻尼比

      
    dt = 0.005          # 时间间隔
    t0 = np.arange(0.0, 40.0, dt)
    step = len(t0)
    dpm = np.zeros( step )

    ug = np.sin(t0)   # 简谐荷载
    

    
    H = np.array( [ [0.0,1.0], [-omg*omg,  -2.0*zeta*omg] ] )
    F = np.array([0.0, 1.0])

    u = precise_integration_solver( H,F,ug,step,dt)

    for i inrange(step):
        t = t0[i]
        res = integrate.quad(func,0,t)[0] # 调用scipy计算杜哈梅积分

        dpm[i] = 1 / omg * res

    fig, (ax0,ax1) = plt.subplots(nrows=2, ncols=1, sharex=True, )
    
    ax0.plot(t0, ug)
    ax0.set( title= '简谐荷载')
    ax0.grid()

    ax1.plot(t0, u, label='精细积分法')
    ax1.plot(t0, dpm, linestyle='--', color='#1661ab',label='杜哈梅积分法')
    ax1.set( title= '位移时程')
    ax1.grid()
    ax1.legend(frameon=True,         # 启用边框
        framealpha=0.8,             # 设置透明度
        fancybox=False,             # 边框圆角与否
        shadow=True,                # 是否启用阴影
        facecolor='lightgrey',      # 背景颜色
        edgecolor='dimgrey',       # 边框颜色
        loc= 'upper left') 

    fig.savefig("11.png",dpi=300)
    plt.show()