python结构_科学网—Python与结构分析(1)---反应谱 - 潘超的博文

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()