• liu448698584
    了解作者
  • matlab
    开发工具
  • 1KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 10 积分
    下载积分
  • 5
    下载次数
  • 2019-03-19 20:14
    上传日期
粘弹性阻尼器等效标准固体力学模型计算分析程序,通过MATLAB编程进行参数识别
Xu_standard_solid_model.rar
  • Xu_standard_solid_model.m
    4.5KB
内容介绍
%% Equivalent standard solid model under sine Wave % Huang Xing-Huai 2018-08-08 clc clear close all hidden % (1) Prepare constants and nodes (this is the longest part of the script): % A = 1; B = 1; C = 1; % coefficients of the Bagley- Torvik equation w=2; %激励频率,excitation frequency T=-10; % temperature of the environmental % dt = 0.002; % step of discretization dt = 0.01/w/2; % step of discretization N = 8001; % number of nodes tmax=(N-1)*dt; % global i gama t X2 t = 0:dt:tmax; % nodes t=t';% T_node = T_node'; % nodes %%%%%%%%%% 9050A viscoelastic material:initial value %%%%%%%% p1=0.0048; %p1 不要调小,非常容易发散。 p1 is sensitive to failure q0=2.7405e6; q1=3.825e6; c=1.34; d=0.58; d=c/2; %the results makes sense only when c=c/2 T0=94; %%%%%%%%%%% ZN22 viscoelastic material:initial value %%%%%%%%% p1=0.0089; %beta(1) p1 不要调小,非常容易发散。 p1 is sensitive to failure q0=2.7853e6; %beta(2) q1=3.3825e5; %beta(3) c=1.18; %beta(4) % c=2; d=0.34; d=c/2; %beta(5) T0=119; %beta(6) alphaT=10^(-12*(T-T0)/(525+(T-T0))); % data fitting by Huang G1=(q0+p1*q1*alphaT^c*w^c)/(1+p1^2*alphaT^c*w^c); G2=(q1-p1*q0)*alphaT^d*w^d*(q0+p1*q1*alphaT^c*w^c)/(q0+p1*q1*alphaT^(2*d)*w^(2*d))/(1+p1^2*alphaT^c*w^c); yeta=G2/G1; yeta=(q1-p1*q0)*alphaT^d*w^d/(q0+p1*q1*alphaT^(2*d)*w^(2*d)); %%% material 9050A experimental test value x=zeros(10,2); x(:,1)=[1;1;2;0.1;0.5;1;5;1;2;1]; %frequency x(:,2)=[-20;-10;-10;0;0;0;0;10;10;20]; % temperature % x1=[1;1;2;0.1;0.5;1;5;1;2;1]; % x2=[-20;-10;-10;0;0;0;0;10;10;20]; y1=[17;5.8;10;2.5;3.3;3.8;7.7;3.0;3.4;2.7]; % G1 Unit: Mpa y2=[1.38;1.39;1.40;0.40;0.90;1.10;1.39;0.71;0.92;0.40]; % ita %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % x=x';y=y';z=z'; GG = @(beta) (beta(2)+beta(1)*beta(3)*10.^(-12*beta(4)*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^beta(4))./(1+beta(1)^2*10.^(-12*beta(4)*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^beta(4)); iita = @(beta) (beta(3)-beta(1)*beta(2))*10.^(-12*beta(4)/2*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^(beta(4)/2)./(beta(2)+beta(1)*beta(3)*10.^(-12*2*beta(4)/2*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^(2*beta(4)/2)); % % % GG0=abs(GG-y1); iita0=abs(iita-y2); % GG0 = @(beta) abs((beta(2)+beta(1)*beta(3)*10.^(-12*beta(4)*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^beta(4))./(1+beta(1)^2*10.^(-12*beta(4)*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^beta(4))-y1); % iita0 = @(beta) abs((beta(3)-beta(1)*beta(2))*10.^(-12*beta(4)/2*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^(beta(4)/2)./(beta(2)+beta(1)*beta(3)*10.^(-12*2*beta(4)/2*(x(:,2)-beta(5))./(525+(x(:,2)-beta(5)))).*x(:,1).^(2*beta(4)/2))-y2); % G_ita2={GG0 iita0}; betalow=[0.08;2;0.03;1;1];%待定系数的预估值 p1, q0, q1, c, T0 beta0up=[0.8;20;3;3;100];%待定系数的预估值: p1, q0, q1, c, T0 beta_bond=[betalow'; beta0up']; alfa2=0.6;beta2=1-alfa2; G_ita=@(beta)(alfa2*sum(((y1-GG(beta))./y1).^2)+beta2*sum(((y2-iita(beta))./y2).^2)); options = optimoptions('ga','PopulationSize',200,'PlotFcn',@gaplotbestf ); [xx,fval,exitflag]= ga(G_ita,5,[],[],[],[],betalow,beta0up,[],options); beta=xx'; % abs((y1-GG(beta))./y1) % abs((y2-iita(beta))./y2) % G_ita2=@(beta)[sum(((y2-iita(beta))./y2).^2);sum(((y1-GG(beta))./y1).^2);]; % options = optimoptions('gamultiobj','PopulationSize',200,'PlotFcn',@gaplotpareto); % [xx22,fval22,exitflag] = gamultiobj(G_ita2,5,[],[],[],[],betalow,beta0up,options); % beta=xx22'; % abs((y1-GG(beta))./y1) % abs((y2-iita(beta))./y2) %%%% the following results are obtained from GA toolbox method, sigle %%%% The p0 of first beta is too small to cause error of Laplace method. % beta=[0.0800024437116451;2.53888368177681;1.79999998694332;1.58103952011942;29.1182879014171]; % beta=[0.0237758950591116;2.84361561185834;0.594148913821854;1.56471129777263;45.6387614669100]; beta=roundn(beta,-3); % beta=x'; figure plot(1:10,y1,'ko',1:10,GG(beta),'b-') legend('Data','Best fit') xlabel('t') ylabel('exp(-tx)') figure plot(1:10,y2,'ko',1:10,iita(beta),'b-') legend('Data','Best fit') xlabel('t') ylabel('exp(-tx)') abs((y1-GG(beta))./y1) abs((y2-iita(beta))./y2) % pause %%%%%%%%%% best value of 9050A viscoelastic material %%%%%%%% p1=beta(1); %p1 不要调小,非常容易发散。 p1 is sensitive to failure q0=beta(2); q1=beta(3); c=beta(4); d=c/2; %the results makes sense only when c=c/2 T0=beta(5); gama0=1; gama=gama0*sin(w*t); figure plot(t,gama) % pause close all hidden
评论
    相关推荐