「Python与地震工程」单自由度体系求解之频域法（傅里叶变换法）

Baron Jean Baptiste Joseph Fourier

原理

频域法（傅里叶变换法）求解体系响应的基本流程

Baron Jean Baptiste Joseph Fourier

图中， p(t)↔P(ω) 为激励的Fourier变换对；u(t)↔U(ω) 为响应的Fourier变换对； H(iω) 为频域响应函数，对于地震作用下的单自由度体系

$H\left( i\omega \right) =\frac{1}{\omega _{0}^{2}}\cdot \frac{1}{1-\left( \frac{\omega}{\omega _0} \right) ^2+i\left( 2\zeta \frac{\omega}{\omega _0} \right)}$

程序代码

# mdof.py
import numpy as np
from numpy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt


plt.style.use("ggplot")
plt.rcParams['font.sans-serif'] = ['Microsoft Yahei']
plt.rcParams['axes.unicode_minus'] = False


next2pow = lambda x: \
    int(round(2**np.ceil(np.log(x)/np.log(2.))))


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_freq(omg, zeta, ag, dt):

    n = len(ag)
    Nfft = next2pow(n)*2 # Fourier变换数据点数取为2的整数次幂并大于等于原始数据点数的2倍 
    af = fft(ag, Nfft) # 快速Fourier变换
    f = fftfreq(Nfft, d=dt) # 离散频率点
    Omg = 2.0*np.pi*f # 离散圆频率点
    H_u = -1.0/(omg**2-Omg**2+2.*zeta*omg*Omg*1.0j) # 频响函数
    u = ifft(af*H_u, Nfft).real # 快速Fourier逆变换并取实部
    v = ifft(af*Omg*H_u, Nfft).real # 快速Fourier逆变换并取实部

    u = u[:n]
    v = v[:n]

    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_freq(omg, zeta, ag, dt)
    draw_response("Seismic Response -- freq", t0, acc0, t, u) # 绘图