import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
# ---------------------中心差分法----------------------- #
def Central_difference_method_displacenment( acceleration_g: list,time:list,step):
length_t = len(time)
displacement = [0.000] * length_t # 定义一个位移向量
speed = [0] * length_t # 定义一个速度向量
acceleration = [0] * length_t # 定义一个加速度向量
for i in range(2, length_t): # 求解位移，单位为mm
displacement[i] = (-acceleration_g[i] + (2 / step ** 2 - 4 * (math.pi) ** 2) * displacement[i - 1] + (
0.2 * math.pi / 2 / step - 1 / step ** 2) * displacement[i - 2]) / (
1 / step ** 2 + 0.2 * math.pi / 2 / step)
for i in range(length_t):
displacement[i]=displacement[i] *1000
displacement1 = displacement[:]
for i in range(length_t):
displacement1[i]=displacement1[i] /1000
for i in range(1,len(displacement)-1): # 求解速度，单位为m/s^2
speed[i]=(displacement1[i+1]-displacement1[i-1])/step/2
for i in range(1, len(displacement) - 1): # 求解加速度，单位为g
acceleration[i] = (displacement1[i + 1] -2*displacement1[i]+ displacement1[i - 1]) / step**2
acceleration[i]=acceleration[i]/9.8
return displacement,speed,acceleration
# ---------------------Newmark方法----------------------- #
def Newmark(acceleration_g:list,time:list,step): # Newmark方法，β取1/6
Newmark_length=len(time)*10
Newmark_displacement = [0]*Newmark_length # 定义一个位移向量
Nemark_speed = [0] * Newmark_length # 定义一个速度向量
Nemark_acceleration = [0]*Newmark_length # 求解位移，单位为mm
p1_m=[0]*(Newmark_length+1)
acceleration_g_linear=[0]*Newmark_length
for i in range(len(acceleration_g)-1):
deta=acceleration_g[i+1]-acceleration_g[i]
for j in range(i*10,i*10+10):
if j ==i:
acceleration_g_linear[j]=acceleration_g[i]
else:
acceleration_g_linear[j] = acceleration_g_linear[j-1]+0.1*deta
a1_m = 6 / (step **2) + 3 / step * 0.2 * math.pi #系数a1与质量m的比值
a2_m = 6 / step + 2 * 0.2 * math.pi #系数a2与质量m的比值
a3_m = 2 + step / 2 * 0.2 * math.pi #系数a3与质量m的比值
k1_m = 4 * math.pi ** 2 + a1_m
for i in range(1,Newmark_length-1):
p1_m[i]=-acceleration_g_linear[i]+a1_m*Newmark_displacement[i-1]+a2_m*Nemark_speed[i-1]+a3_m*Nemark_acceleration[i-1]
Newmark_displacement[i]=p1_m[i]/k1_m # 求解位移
Nemark_speed[i]=3/step*(Newmark_displacement[i]-Newmark_displacement[i-1])-2*Nemark_speed[i-1]-step/2
*Nemark_acceleration[i-1] # 求解速度
Nemark_acceleration[i]=6/step**2*(Newmark_displacement[i]-Newmark_displacement[i-1])-6/step*Nemark_speed[i-1]
-2*Nemark_acceleration[i-1] # 求解加速度
for i in range(len(Newmark_displacement)):
Newmark_displacement[i]=Newmark_displacement[i]*1000
for i in range(len(Nemark_acceleration)):
Nemark_acceleration[i]=Nemark_acceleration[i]/9.8
return Newmark_displacement, Nemark_speed, Nemark_acceleration
def Newmark_(acceleration_g:list,time:list,step): # Newmark方法，β取1/6
Newmark_length=len(time)
Newmark_displacement = [0]*Newmark_length # 定义一个位移向量
Nemark_speed = [0] * Newmark_length # 定义一个速度向量
Nemark_acceleration = [0]*Newmark_length # 求解位移，单位为mm
p1_m=[0]*(Newmark_length+1)
a1_m = 6 / (step **2) + 3 / step * 0.2 * math.pi #系数a1与质量m的比值
a2_m = 6 / step + 2 * 0.2 * math.pi #系数a2与质量m的比值
a3_m = 2 + step / 2 * 0.2 * math.pi #系数a3与质量m的比值
k1_m = 4 * math.pi ** 2 + a1_m
for i in range(1,Newmark_length-1):
p1_m[i]=-acceleration_g[i]+a1_m*Newmark_displacement[i-1]+a2_m*Nemark_speed[i-1]+a3_m*Nemark_acceleration[i-1]
Newmark_displacement[i]=p1_m[i]/k1_m # 求解位移
Nemark_speed[i]=3/step*(Newmark_displacement[i]-Newmark_displacement[i-1])-2*Nemark_speed[i-1]-step/2
*Nemark_acceleration[i-1] # 求解速度
Nemark_acceleration[i]=6/step**2*(Newmark_displacement[i]-Newmark_displacement[i-1])-6/step*Nemark_speed[i-1]
-2*Nemark_acceleration[i-1] # 求解加速度
for i in range(len(Newmark_displacement)):
Newmark_displacement[i]=Newmark_displacement[i]*1000
for i in range(len(Nemark_acceleration)):
Nemark_acceleration[i]=Nemark_acceleration[i]/9.8
return Newmark_displacement, Nemark_speed, Nemark_acceleration
# ---------------------调幅----------------------- #
def amplitude_modulation(acceleration_g,acceleration):
for i in range(len(acceleration_g)):
acceleration[i]=abs(acceleration[i]+acceleration_g[i]/9.8)
max_acc=2.5/max(acceleration)
for i in range(len(acceleration_g)):
acceleration_g[i]=acceleration_g[i]*max_acc
return acceleration_g
if __name__=="__main__":
# ---------------------------H1方向---------------------------------- #
data_H1= pd.read_excel('RSN-59.xlsx', 'H1')
acceleration_H1_g = list(data_H1['a'])
for i in range(len(acceleration_H1_g)):
acceleration_H1_g[i] = acceleration_H1_g[i] * 9.8 # m/s^2
t_H1 = data_H1['t']
t_H1_Newmark=np.arange(0,14.75,0.001)
displacement_H1, speed_H1, acceleration_H1 = Newmark_(acceleration_H1_g, t_H1, 0.01)
acceleration_g_H1=amplitude_modulation(acceleration_H1_g,acceleration_H1)
Cdisplacement_H1, Cspeed_H1, Cacceleration_H1 = Central_difference_method_displacenment(acceleration_g_H1,t_H1,0.01)
Newmark_displacement_H1, Newmark_speed_H1, Newmark_acceleration_H1 = Newmark(acceleration_g_H1, t_H1, 0.001)
plt.figure(num=1, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_H1, Cdisplacement_H1,'*',label='Central difference method')
plt.plot(t_H1_Newmark, Newmark_displacement_H1,label='Newmark')
plt.title('Displacement time history curve of Horizontal-1', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('displacement(mm)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right',fontsize=15)
plt.figure(num=2, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_H1, Cspeed_H1,'*',label='Central difference method')
plt.plot(t_H1_Newmark, Newmark_speed_H1,label='Newmark')
plt.title('velocity time history curve of Horizontal-1', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('velocity(m/s)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.figure(num=3, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_H1, Cacceleration_H1,'*',label='Central difference method')
plt.plot(t_H1_Newmark, Newmark_acceleration_H1,label='Newmark')
plt.title('Acceleration time history curve of Horizontal-1', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('acceleration(g)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.show()
# ---------------------------H2方向---------------------------------- #
data_H2 = pd.read_excel('RSN-59.xlsx', 'H2')
acceleration_H2_g = list(data_H2['a'])
for i in range(len(acceleration_H2_g)):
acceleration_H2_g[i] = acceleration_H2_g[i] * 9.8
t_H2 = data_H2['t']
t_H2_Newmark = np.arange(0, 14.75, 0.001)
displacement_H2, speed_H2, acceleration_H2 = Central_difference_method_displacenment(acceleration_H2_g, t_H2, 0.01)
acceleration_g_H2= amplitude_modulation(acceleration_H2_g, acceleration_H2)
Cdisplacement_H2, Cspeed_H2, Cacceleration_H2 = Central_difference_method_displacenment(acceleration_g_H2,t_H2,0.01)
Newmark_displacement_H2, Newmark_speed_H2, Newmark_acceleration_H2 = Newmark(acceleration_g_H2, t_H2, 0.001)
plt.figure(num=4, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_H2, Cdisplacement_H2, '*',label='Central difference method')
plt.plot(t_H2_Newmark, Newmark_displacement_H2,label='Newmark')
plt.title('Displacement time history curve of Horizontal-2', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('displacement(mm)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.figure(num=5, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_H2, Cspeed_H2, '*',label='Central difference method')
plt.plot(t_H2_Newmark, Newmark_speed_H2,label='Newmark')
plt.title('velocity time history curve of Horizontal-2', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('velocity(m/s)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.figure(num=6, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_H2, Cacceleration_H2, '*',label='Central difference method')
plt.plot(t_H2_Newmark, Newmark_acceleration_H2,label='Newmark')
plt.title('Acceleration time history curve of Horizontal-2', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('acceleration(g)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.show()
# ---------------------------V方向---------------------------------- #
data_V = pd.read_excel('RSN-59.xlsx', 'V')
acceleration_V_g = list(data_V['a'])
for i in range(len(acceleration_V_g)):
acceleration_V_g[i] = acceleration_V_g[i] * 9.8
t_V = data_V['t']
t_V_Newmark = np.arange(0, 14.75, 0.001)
displacement_V, speed_V, acceleration_V = Central_difference_method_displacenment(acceleration_V_g, t_V, 0.01)
acceleration_g_V = amplitude_modulation(acceleration_V_g, acceleration_V)
Cdisplacement_V, Cspeed_V, Cacceleration_V = Central_difference_method_displacenment(acceleration_g_V, t_V,0.01)
Newmark_displacement_V, Newmark_speed_V, Newmark_acceleration_V = Newmark(acceleration_g_V, t_V, 0.001)
plt.figure(num=7, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_V, Cdisplacement_V, '*',label='Central difference method') # 绘制位移时程曲线
plt.plot(t_V_Newmark, Newmark_displacement_V,label='Newmark')
plt.title('Displacement time history curve of Vertical', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('displacement(mm)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.figure(num=8, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_V, Cspeed_V, '*',label='Central difference method')
plt.plot(t_V_Newmark, Newmark_speed_V,label='Newmark')
plt.title('velocity time history curve of Vertical', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('velocity(m/s)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.figure(num=9, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
plt.plot(t_V, Cacceleration_V, '*',label='Central difference method')
plt.plot(t_V_Newmark, Newmark_acceleration_V,label='Newmark')
plt.title('Acceleration time history curve of Vertical', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('acceleration(g)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.legend(loc='upper right', fontsize=15)
plt.show()
plt.figure(num=10, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
for i in range(len(acceleration_g_H1)):
acceleration_g_H1[i]=acceleration_g_H1[i]/9.8
plt.plot(t_V, acceleration_g_H1)
acceleration_g_H1_min = np.argmin(acceleration_g_H1)
acceleration_g_H1_max = np.argmax(acceleration_g_H1)
plt.title(' PEER Ground Motion Database-H1(PGA=-1.84g)', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('acceleration(g)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.show()
plt.figure(num=11, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
for i in range(len(acceleration_g_H2)):
acceleration_g_H2[i]=acceleration_g_H2[i]/9.8
plt.plot(t_V, acceleration_g_H2)
plt.title(' PEER Ground Motion Database-H2(PGA=1.80g)', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('acceleration(g)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.show()
plt.figure(num=12, figsize=(15, 10), dpi=200, edgecolor='black', linewidth=0.2, frameon=False)
for i in range(len(acceleration_g_V)):
acceleration_g_V[i]=acceleration_g_V[i]/9.8
plt.plot(t_V, acceleration_g_V)
plt.title('PEER Ground Motion Database-V(PGA=1.56g)', fontsize=25)
plt.xticks(fontsize=15), plt.yticks(fontsize=15)
plt.ylabel('acceleration(g)', fontsize=20)
plt.xlabel('time(s)', fontsize=20)
plt.show()
