bouc_wen.zip_Bouc wen_bouc_bouc-wen模型_wen_wen model

import numpy as np import matplotlib.pylab as plt import xlrd, xlwt import copy class Parameters(): def __init__(self, alpha, k0, Fy, A, n, beta, gamma, Sk, Sp, Sc, Sb, Sm, Sl, mu_c): self.alpha = alpha self.k0 = k0 self.Fy = Fy self.A = A self.n = n self.beta = beta self.gamma = gamma self.Sk = Sk self.Sp = Sp self.Sc = Sc self.Sb = Sb self.Sm = Sm self.Sl = Sl self.mu_c = mu_c self.mu_y = Fy/k0 def Eta( mu, mu_max, mu_min, Sk ): if mu >= 0 : e = 1.0 + Sk * ( mu_max + mu ) / 2 if mu < 0 : e = 1.0 + Sk * ( abs(mu_min) + abs(mu) ) /2 return e def Mu(x, mu_y, mu_c): mu = x / mu_y for i in range(len(mu)): if abs(mu[i]) >= mu_c: mu[i] = np.sign(mu[1])*mu_c return mu def BW(self,x,dt): mu = Mu(x, self.mu_y, self.mu_c) d_mu = np.zeros(mu.shape) for i in range(len(mu)-1): d_mu[ i ] = (mu[ i + 1 ] - mu[ i ])/dt z = np.zeros(d_mu.shape) dz = np.zeros(d_mu.shape) mu_max = np.max(mu) mu_min = np.min(mu) eta = Eta (mu[0], mu_max, mu_min, self.Sk) dz[0] = d_mu[0] * ( self.A - abs(z[0])**self.n * ( self.beta * np.sign(d_mu[0] * z[0] ) + self.gamma ) ) \ / (eta * ( 1 +abs(mu[0])**self.Sm * self.Sc / self.Sb * np.sqrt( 2/np.pi ) * np.exp(-( z[0] / self.Sb )**2 )) * \ ( self.A - abs(z[0])**self.n * (self.beta * np.sign(d_mu[0] * z[0]) - self.gamma ) )) z[1] = z[0] + dz[0] * dt eta = Eta(mu[1], mu_max, mu_min, self.Sk) dz[1] = d_mu[1] * ( self.A - abs(z[1])**self.n * ( self.beta * np.sign(d_mu[1] * z[1] ) + self.gamma ) ) \ / (eta * ( 1 +abs(mu[1])**self.Sm * self.Sc / self.Sb * np.sqrt( 2/np.pi ) * np.exp( -( z[1] / self.Sb )**2 )) * \ (self.A - abs(z[1])**self.n * (self.beta * np.sign(d_mu[1] * z[1]) - self.gamma ) )) for i in range(2,len(mu)): eta = Eta(mu[i], mu_max, mu_min, self.Sk) z[i] = z[i-1] + dz[i-1] * dt + ( dz[i-1] - dz[i-2] )/ 2 * dt dz[i] = d_mu[i] * ( self.A - abs(z[i])**self.n * ( self.beta * np.sign(d_mu[i] * z[i] ) + self.gamma ) ) \ / (eta * ( 1 +abs(mu[i])**self.Sm * self.Sc / self.Sb * np.sqrt( 2/np.pi ) * np.exp( -( z[i] / self.Sb )**2 )) * \ ( self.A - abs(z[i])**self.n * (self.beta * np.sign(d_mu[i] * z[i]) - self.gamma ) )) F = self.alpha*self.k0*x+(1-self.alpha) * self.Fy * z return x, F def save_data(F, x): n = len(F) print(n) book = xlwt.Workbook() sheet = book.add_sheet('F-x', cell_overwrite_ok=True) for i in range(n-1): sheet.write(i, 0, x[i+1]) sheet.write(i, 1, F[i+1]) for i in range(n-1): sheet.write(i, 2, x[i]) sheet.write(i, 3, F[i]) for i in range(n-1): dx = (x[i+1] - x[i])/0.02 sheet.write(i, 4, dx) book.save('test/hysteretic_curve.xls') def f_r(x): dx_1 = np.zeros(len(x)) for i in range(len(x) - 1): dx_1[i] = (x[i + 1] - x[i]) / 0.02 r = np.zeros(x.shape) for i in range(1, len(x)): r[i] = r[i-1] + 1.01e6*(x[i] - x[i-1]) + 0.02*1.39e6*dx_1[i] - 0.02*3.18e3*np.abs(dx_1[i])*np.abs(r[i-1]-1.01e6*x[i-1])**1.09 * \ np.sign(r[i-1]-1.01e6*x[i-1]) - 0.02*41.49*dx_1[i]*np.abs(r[i-1] - 1.01e6*x[i-1])**1.09 return r def main(): y = np.zeros(1000) for i in range(1000): y[i] = 2000 / 2000000 * np.sin(np.pi / 250.0 * i) # x1 = np.zeros(500) # for i in range(500): # x1[i] = i / 500 # # x2 = copy.deepcopy(x1[::-1]) # x3 = copy.deepcopy(-x1) # x4 = copy.deepcopy(-x2) # # x = np.zeros(2000) # x[:500] = x1 # x[500:1000] = x2 # x[1000:1500] = x3 # x[1500:2000] = x4 # # y = np.zeros(20000) # for i in range(5): # y[i * 2000:(i + 1) * 2000] = (i*3.2)/ 800 * x # for i in range(5, 10): # y[i * 2000:(i + 1) * 2000] = 14.4 / 800 * x parameters = Parameters( alpha=0.009, k0=4.09e6, Fy=11170, A=1, n=2, beta=0.95, gamma=0.05, Sk=0.3, Sp=0.2, Sc=0.2, Sb=0.1, Sm=1.0, Sl=0.2, mu_c = 16.4 ) # _, F = BW(parameters, y, 0.02) F = f_r(y) plt.rcParams['savefig.dpi'] = 600 # 图片像素 plt.rcParams['figure.dpi'] = 600 # 分辨率 plt.plot(np.arange(len(F))*0.02, y) plt.xlabel('Time(s)') plt.ylabel('Curvature(1/m)') plt.show() M = F C = y plt.plot(C, M) plt.xlabel('Curvature(1/m)') plt.ylabel('Moment(kN.m)') plt.show() save_data(F, y) if __name__ == "__main__": main()