动力方程的精细积分法

非齐次动力学方程常用的数值积分方法有中心差分法、Wilson类方法、Houbolt法、Newmark类方法等，虽然应用十分广泛，但是在长时间的积分过程中会表现出能量耗散和相位误差。钟万勰院士提出的精细积分法，由于其具有高精度和强稳定性的特点，为非齐次动力方程的时程积分提供了一种新的计算途径。

微分方程的降阶

结构动力学方程的一般形式可表达为

式中，  和  分别表示系统的质量矩阵、阻尼矩阵和刚度矩阵，  表示施加于系统的外力。当质量矩阵正定时，引入变量  或

令新的状态变量为  ，式(1)可以重新写为

式中  为常数矩阵,  为非线性广义外力项.

在一个积分步长  内，记  ，(2)两边同时乘以  ，非齐次动力方程(2)的解为

指数矩阵的计算

(4)中有  需要计算。根据指数的运算

由于  本来是不大的时间区段,则  将是非常小的一个时间区段(通常取  )。采用  级数展开(4项)近似计算矩阵指数

其中,  为单位矩阵。(6)代入(5)

由(6)可知，  是一个小量的矩阵，当它与单位阵  相加时，在计算机舍入操作时将很可能被丢失掉。考虑递推关系

其中

(9)的N次矩阵乘法迭代相当于

  for(k =0; k <N; k++) {
    Ta = 2Ta + TaTa
}

 为常数矩阵,求解过程中只计算一次。当以上循环语句全部计算结束后，然后再计算:

即可求出指数矩阵  的值。参考《计算机科学计算》第七章。

非齐次方程特解的计算

由于  为常数矩阵，对于(1)的齐次方程的通解形式为：

由上述方法精细计算得矩阵  后，时程积分就变为:

而对于非齐次方程(2)，若非齐次项  在时间步  内为线性，即方程为:

当  时，  ，  可表示为:

凑微分的时候，

算例

如图所示的结构，设横梁为无限刚性,体系的质量全部集中在各横梁上,各层间侧移刚度为1，自振频率分别为  。用精细积分法计算位移时程。

  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(20):
        T_a = 2.0*T_a + T_a @ T_a

    T = np.eye(n) + T_a  

    return T 


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

    ug = np.loadtxt("EQ-S-3.txt") # 读取地震波    
    dt = 0.005# 时间间隔
    n = len(ug)
    t0 = np.linspace(0.0, dt*(n-1), n)

    eta = 0.05

    zeta = 0.05         # 固有阻尼比
    omg1 = 0.386
    omg2 = 1.036

    K = np.array([ [2,-1,0],[-1,2,-1],[0,-1, 1] ])
    M = np.array([ [1.5, 0, 0],[0, 1, 0],[0, 0, 1.5] ])

    # 一阶和二阶频率以及瑞利阻尼构造阻尼矩阵
    a0 = 2*zeta*omg1*omg2 / (omg1+omg2 )
    a1 = 2*zeta/ (omg1+omg2 )
    C  = a0*M + a1*K

    F = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0])

    invM = np.linalg.inv(M)

    H0 =  np.dot(C, invM)
    H1 =  np.dot(invM, C)
    H2 = 0.25 * np.dot(C, H1) -K

    H3 =  np.hstack( (-0.5*H1, invM) )
    H4 =  np.hstack( (H2, -0.5*H0) )

    H = np.vstack( (H3, H4) )
    invH = np.linalg.inv(H)

    T = compute_t(6, H, eta)

    x  = np.zeros((n, 6))
    xn = np.zeros(6)

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

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

    u = x[:,0]

    fig, (ax0,ax1) = plt.subplots(nrows=2, ncols=1, sharex=True,
                                    figsize=(16, 6) )
    
    ax0.plot(t0, ug)
    ax0.set( title= '地震波加速度时程')
    ax0.grid()

    ax1.plot(t0, u)
    ax1.set( title= '位移时程')
    ax1.grid()

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


if __name__ == "__main__":
    main( )