Python基于中心差分及Newmark法计算地震波反映谱

发表于 更新于 Waline： 阅读次数：

Python基于中心差分及Newmark法计算地震波反映谱

地震工程作业代码

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
import os
import math
import numpy as np
import matplotlib.pyplot as plt
import scipy

class spectrum:
    sw=np.array([]) #地震波a
    v=np.array([]) #地震波速度
    su=np.array([]) #地震波位移
    tf=np.array([]) #调幅后地震波a
    
    u=np.array([])  
    u0=np.array([])
    u00=np.array([])
    
    saT=np.array([])
    sa=np.array([]) #谱加速度
    sv=np.array([]) #谱速度
    sd=np.array([]) #谱位移

    PGA=0
    PGD=0
    PGV=0
    t=0
    
    def __init__(self,path,t):#路径，步长
        file_4=open(path,"r")
        str4="" #文件前四行
        for p in range(0,4):
            str4=str4+file_4.readline()
        data=file_4.readlines()# 读取后面全部数据
        file_4.close() #获取前四行
        lit=[]
        for d in data:
            tmpd=d.split()
            for num_tmpd in range(0,len(tmpd)):
                lit.append(float(tmpd[num_tmpd]))
        self.sw=np.array(lit)
        self.t=t
        return
    
    def newmark_beta(self,T,ksi=0.05,gamma=0.5,beta=0.25):#步长，周期，地震波
        if T==0:
            T=0.00000001
        n=len(self.sw)
        m=1.0
        c=2*m*ksi*(2*np.pi/T)
        #分别为加速度，速度，位移
        u00=np.zeros(n)
        u0=np.zeros(n)
        u=np.zeros(n)
        k=(2*np.pi/T)**2*m
        u00[0]=-self.sw[0]-c*u0[0]-k*u[0]
        a1=m/(beta*(self.t**2))+gamma*c/(beta*self.t)
        a2=m/(beta*self.t)+(gamma/beta-1)*c
        a3=(1/(2*beta)-1)*m+self.t*(gamma/(2*beta)-1)*c
        k_hat=k+a1 #k^/hat

        for i in range(1,n):
            p_hat=-self.sw[i]+a1*u[i-1]+a2*u0[i-1]+a3*u00[i-1]
            u[i]=p_hat/k_hat
            u0[i]=gamma/(beta*self.t)*(u[i]-u[i-1])+(1-gamma/beta)*u0[i-1]+self.t*(1-gamma/(2*beta))*u00[i-1]
            u00[i]=1/(beta*(self.t**2))*(u[i]-u[i-1])-u0[i-1]/(beta*self.t)-(1/(2*beta)-1)*u00[i-1]
        self.u=u
        self.u0=u0
        self.u00=u00+self.sw
        return u,u0,u00+self.sw
    
    def central_difference(self,T,ksi=0.05,gamma=0.5,beta=0.25):
        m=1.0
        n=len(self.sw)
        u00=np.zeros(n)
        u0=np.zeros(n)
        u=np.zeros(n+1)
        c=2*m*ksi*(2*np.pi/T)
        k=(2*np.pi/T)**2*m
        u00[0]=-self.sw[0]-c*u0[0]-k*u[0]
        u_1=u[0]-self.t*u0[0]+self.t**2*u00[0]/2
        k_hat=m/(self.t**2)+c/(2*self.t)
        a=m/(self.t**2)-c/(2*self.t)
        b=k-2*m/(self.t**2)
        for i in range(n):
            if i==1:
                p_hat=-self.sw[i]-a*u_1-b*u[i]
            else:
                p_hat=-self.sw[i]-a*u[i-1]-b*u[i]
            u[i+1]=p_hat/k_hat
            u0[i]=(u[i+1]-u[i-1])/(2*self.t)
            u00[i]=(u[i+1]-2*u[i]+u[i-1])/(self.t**2)
        self.u=u[:n]
        self.u0=u0
        self.u00=u00+self.sw
        return
            
    def get_PGA(self):
        self.PGA=max(abs(self.sw))
        return max(abs(self.sw))
    
    def get_PGV(self):
        self.PGV=max(abs(self.v))
        return max(abs(self.v))
    
    def get_PGD(self):
        self.PGD=max(abs(self.su))
        return max(abs(self.su))
    
    def get_sa(self,begin,end,step):
        T=np.arange(begin,end,step)
        sa=np.array([])
        sv=np.array([])
        su=np.array([])
        for i in T:
            u,v,a=self.newmark_beta(i)
            sa=np.append(sa,max(abs(a)))
            sv=np.append(sv,max(abs(v)))
            su=np.append(su,max(abs(u)))
        self.saT=T
        self.sa=sa
        self.sd=su
        self.sv=sv
        return
    
    def get_v(self):#步长，地震波
        v=[]
        for i in range(len(self.sw)):
            v.append(np.trapz(self.sw[:i+1],dx=self.t))
        self.v=np.array(v)
        return v
    
    def get_u(self):
        v=self.get_v()
        u=[]
        for i in range(len(self.sw)):
            u.append(np.trapz(v[:i+1],dx=self.t))
        self.su=np.array(u)
        return u
    
    def tiaofu_sa(self,T,target):# 地震波步长 地震波 要调幅到的周期和值'
        u,v,a=self.newmark_beta(T)
        sa=max(abs(a))
        self.tf=self.sw
        self.sw=self.sw*(target/sa)
        self.get_v()
        self.get_u()
        return target/sa
    
    def fypa(self,T,x,path):
        y=np.loadtxt(path)
        plt.rcParams['font.family'] = ['sans-serif']
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.figure(figsize=(10,5))
        plt.title("标准加速度反应谱")
        self.get_PGA()
        plt.plot(T,x/self.PGA,label="python")
        plt.plot(T,y/self.PGA,label="SPECTR")
        plt.xlabel("T(s)")
        plt.ylabel("$\\beta$")
        plt.legend()
        
    def fypv(self,T,x,path):
        y=np.loadtxt(path)/1000
        plt.rcParams['font.family'] = ['sans-serif']
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.figure(figsize=(10,5))
        plt.title("标准速度反应谱")
        self.get_PGV()
        plt.plot(T,x/self.PGV,label="python")
        plt.plot(T,y/self.PGV,label="SPECTR")
        plt.xlabel("T(s)")
        plt.ylabel("$\\beta$")
        plt.legend()
        
    def fypu(self,T,x,path):
        y=np.loadtxt(path)/1000
        plt.rcParams['font.family'] = ['sans-serif']
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.figure(figsize=(10,5))
        plt.title("标准位移反应谱")
        self.get_PGD()
        plt.plot(T,x/self.PGD,label="python")
        plt.plot(T,y/self.PGD,label="SPECTR")
        plt.xlabel("T(s)")
        plt.ylabel("$\\beta$")
        plt.legend()
        
    def sw_a(self):
        plt.rcParams['font.family'] = ['sans-serif']
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.rcParams['axes.unicode_minus']=False # 用来正常显示负号
        plt.figure(figsize=(10,5))
        plt.title("地震波加速度时程")
        t=np.linspace(0,len(self.sw)*self.t,len(self.sw))
        plt.plot(t,self.sw)
        plt.xlabel("t(s)")
        plt.ylabel("a(g)")
        #plt.legend()

使用实例

1
2
3
4
5
6
# 创建对象，并读取PEER地震波初始化，第一个参数为路径，第二个参数是地震波的步长，没有自动读取
n=spectrum(r"C:\ATC63\ATC63\远场\960\RSN960_NORTHR_LOS-UP.AT2",0.01)
# 使用Sa调幅，第一个参数为周期，第二个参数为要调幅的目标值 
tt=n.tiaofu_sa(2,0.5)
# 与其他软件输出的反应谱进行对比，其实前两个参数完全冗余，不想改了 
n.fypa(n.saT,n.sa,"esa.csv")

说明

  1. 写成了一个对象，整理成文档发现小问题好多呀，不想改了怎么办
  2. 中心差分法和Newmark-beta本身没有大问题，因为我是做标准反应谱，即除以PGA、PGV和PGD的谱，所以并没有乘以9.8依旧数据没问题，如果不是要标准谱的话，请自行乘以g值
  3. 还是有时间我自己好好写写吧

ToDo

  1. g值写进去
  2. 增加PGA调幅
  3. 自动读取PEER地震波的时程步长
  4. 删除类内输出反应谱函数的冗余参数
  5. 增加其他常见波参数，并直接输出excel方便查看所有地震波常见信息