||
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()
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 10:21
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社