「Python與地震工程」單自由度體系求解2

頻域法求解體系響應的基本流程

圖中， 為激勵的Fourier變換對； 為響應的Fourier變換對； 頻域響應函數，對於地震作用下的單自由度體系

以下是Python代碼

import scipy as spfrom numpy.fft import fft, ifft, fftfreqimport matplotlib.pyplot as pltnext2pow = lambda x: int(round(2**sp.ceil(sp.log(x)/sp.log(2.))))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 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(omg0 = 1.0*2.0*sp.pi, zeta = 0.05): acc0 = sp.loadtxt("EQ-S-3.txt") # 讀取地震波 dt = 0.005 # 時間間隔 n = len(acc0) t0 = sp.linspace(0.0,dt*(n-1),n) Nfft = next2pow(n)*2 # Fourier變換數據點數取為2的整數次冪並大於等於原始數據點數的2倍 af = fft(acc0, Nfft) # 快速Fourier變換 f = fftfreq(Nfft, d=dt) # 離散頻率點 omg = 2.0*sp.pi*f # 離散圓頻率點 H = -1.0/(omg0**2-omg**2+2.*zeta*omg0*omg*1.0j) # 頻響函數 u = ifft(af*H, Nfft).real # 快速Fourier逆變換並取實部 # 顯示地震結束後一段時間內的自由振動衰減情況 ne = round(n*1.2) t = sp.linspace(0.0,dt*(ne-1),ne) u = u[:ne] draw_response("Seismic Response -- freq", t0, acc0, t, u) # 繪圖if __name__ == __main__: solve_sdof_eqwave_freq()

文中所用人工波（EQ-S-3.txt）可前往下面鏈接下載：

https://pan.baidu.com/s/1SubFZHUB3t90gjdATfD7Rw?

pan.baidu.com

推薦閱讀：