基于Clough–Penzien谱模拟平稳加速度时程和非平稳地震序列.rar

import numpy as np import matplotlib.pyplot as plt def wg(t,w0,a,T): wg=w0-a*t/T return wg def xig(t,xi0,b,T): xig=xi0-b*t/T return xig def S0(t,amax,r,wg,xig): S0=2*amax**2/(r**2*np.pi*wg*(2*xig+1/(2*xig))) return S0 def A2(t,c,d): A2=(t/c*np.exp(1-t/c))**d return A2 def Sxg(t,w,wg,xig,wf,xif,A2,S0): Sxg=A2*\ (wg**4+4*xig**2*wg**2*w**2)/((w**2-wg**2)**2+4*xig**2*wg**2*w**2)*\ (w**4)/((w**2-wf**2)**2+4*xif**2*wf**2*w**2)\ *S0 return Sxg def XkYk(N): j=np.arange(N) theta=np.random.uniform(-np.pi,np.pi) k=j Xk=np.sqrt(2)*np.cos(k*theta+np.pi/4) Yk=np.sqrt(2)*np.cos(k*theta+np.pi/4) return Xk,Yk def XkYk_double(nsel): Theta1=np.random.uniform(0,2*np.pi,nsel) Theta2=np.random.uniform(0,2*np.pi,nsel) Xk=np.cos(nsel*Theta1)+np.sin(nsel*Theta1) Yk=np.cos(nsel*Theta2)+np.sin(nsel*Theta2) return Xk,Yk def Xg_nonstationary(t,deltaw,wk,Xk,Yk,Sxg): P=(np.cos(wk*t)).T*Xk+(np.sin(wk*t)).T*Yk print(P.shape) O=np.sqrt(2*Sxg*deltaw) print(O.shape) X_m=P.T*O print(X_m.shape) return X_m w0=11;xi0=0.7;a=8;b=0.15;c=7;d=2;T=30;r=2.65;nsel=1601;deltaw=0.2;amax=200 #单变量 _Xn=XkYk(nsel)[0];_Yn=XkYk(nsel)[0] kkk=np.random.permutation(nsel)#随机打乱一个数字序列 Xk=_Xn[kkk] Yk=_Yn[kkk] wg=5*np.pi;xig=0.60;wf=0.5*np.pi;xif=0.60;amax=200;r=2.8;wu=240;T=30 deltaw=0.15; tt=np.arange(0,T,0.01);k=np.arange(0,nsel);w=k*deltaw t,wk=np.meshgrid(tt,w) def Sug(S0,wg,xig,wf,xif,w): Sug=S0*\ (wg**4+4*xig**2*wg**2*w**2)/((w**2-wg**2)**2+4*xig**2*wg**2*w**2)*\ (w**4)/((w**2-wf**2)**2+4*xif**2*wf**2*w**2) return Sug def Xg_stationary(t,deltaw,wk,Xk,Yk,Sug): P=(np.cos(wk*t)).T*Xk+(np.sin(wk*t)).T*Yk print(P.shape) O=np.sqrt(2*Sug*deltaw) print(O.shape) X_m=O*P print(X_m.shape) return X_m S0=amax**2/(r**2*(np.pi*wg*(2*xig+1/(2*xig)))) Sug1=Sug(S0,wg,xig,wf,xif,w) #平稳时间序列的生成 X_sat=Xg_stationary(t,deltaw,wk,Xk,Yk,Sug1) result_sationary=np.sum(X_sat,axis=1) import pandas as pd Station_series=pd.DataFrame(result_sationary,index=tt,columns=['Stationary acceleration']) plt.figure(2) Station_series.plot(figsize=(6,2),color='k',lw=1,xlabel='Time',ylabel=r'$Accelerate (cm/s^2)$') plt.tight_layout() plt.show(block=True)