Python for Structural Analysis(1)

---Response Spectrum

对于结构地震分析而言，地震波的选择非常重要。而地震波选择的主要依据是地震波的频谱特性，频谱特性在实际应用中多以傅里叶谱（幅值谱、相位谱、功率谱等）及反应谱等来描述。傅里叶谱可以简单地用Matlab中的fft()函数得到，但工程中使用更为广泛的反应谱的计算则不是很方便，（用SeismoSignal可以很容易得到反应谱，但目前的版本只能算周期0~4s的反应谱；据说Midas也可以算但没用过）。早就想自己DIY一个反应谱的小程序，但苦于编程水平太一般，搁置许久。近来因学习ABAQUS而认识了Python，感觉非常好用，第一个Python程序就先从反应谱开始吧。

其实原理很简单，就是调用了Scipy.signal模块的几个函数，功能就是把同一文件夹下的地震波文件（时间-加速度）读入生成反应谱并保存为文本。计算结果倒是准确，就是慢了点儿。我取了上海人工波和著名的El-centro波试算了一下，结果如下图所示：

下面把代码晒出来献献丑：

from glob import *
from scipy import *
from scipy.signal import *
from matplotlib.pyplot import *
for wf in glob('*.txt'):
    dr = 0.05
    eq = loadtxt(wf)
    plot(eq[:,0],eq[:,1])
    show()  #绘制地震波
    t0,t,dt = 0.01,6,0.05  #设定反应谱的周期范围及间隔
    T = arange(t0,t+dt,dt)
    s=zeros(len(T))
    w=2.*pi/T
    num = (1.)
    for i in range(0,len(w)):
        den = (1.,2.*dr*w[i],w[i]*w[i])
        SYS = (num, den)
        t1, y, x = lsim(SYS,eq[:,1],eq[:,0])
        s[i]=max(abs(y))
    v=s*w    #伪速度谱
    a=v*w    #伪加速度谱
    X = zeros((len(T),4))
    X[:,0]=T;X[:,1]=a;X[:,2]=v;X[:,3]=s
    savetxt(wn.replace('.txt','.dat'), X, fmt='%4.4f', delimiter=' ')
    plot(T,a)  #绘制伪加速度谱
    show()