CHE.rar_highestn3k_matlab刚度退化_滞回曲线 matlab_累积滞回耗能_退化 matlab

clc; close all; clear; %%%%%%%%%通过实验数据直接生成骨架曲线 不包含弹塑性时程分析程序%%%%%%%%%%% data= xlsread('F:\MATLAB WorkSpace\WangTao\试验数据.xlsx');%读取数据 D=data(:,1);% unit: mm位移 F=data(:,2);% unit: KN力 outdata=[D F]; LineNum=size(outdata,1); EndCircle=1; %% %==============找各级循环分割点位置=================% for ii=3:LineNum if(outdata(ii,1)*outdata(ii-1)<=0)&&(outdata(ii-1,1)<0) %默认初始加载从正向开始 EndCirclePoints(EndCircle)=ii;%分割点的行号EndCirclePoints EndCircle=EndCircle+1; %中间分割点数为EndCircle-1 end end %==============检查最后是否存在不完整滞回环=========% if EndCirclePoints(EndCircle-1)<LineNum LoopNum=length(EndCirclePoints)+1; %LoopNum 仅为完整循环的个数。还有最后一个不完整循环 else LoopNum=length(EndCirclePoints); end %==================分割各个滞回环==================% for k=1:LoopNum if k==1 LoopCircles(k)={outdata(1:EndCirclePoints(k),:)}; elseif k<LoopNum LoopCircles(k)={outdata(EndCirclePoints(k-1):EndCirclePoints(k),:)}; else LoopCircles(k)={outdata(EndCirclePoints(k-1):LineNum,:)}; %后续步骤中，最后一段应区分是少于半周还是少于整周 end end %========提取各个滞回环荷载骨架曲线点，储存在矩阵 FramePointsPositive 和 FramePointsNegative 中 FramePointsPositive=zeros(LoopNum,2);FramePointsNegative=zeros(LoopNum,2); for k=1:LoopNum A=LoopCircles{k}; if k<LoopNum [ColMaxValue,LineMax]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax [ColMinValue,LineMin]=min(A); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax FramePointMax=A(LineMax(2),:); FramePointMin=A(LineMin(2),:); FramePointsPositive(k,:)=FramePointMax; %骨架点储存在FramePoints矩阵中 FramePointsNegative(k,:)=FramePointMin; %======区分最后一段是完整周，多于半周，还是少于半周=========% else if EndCirclePoints(EndCircle-1)==LineNum % 完整周 [ColMaxValue,LineMax]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax [ColMinValue,LineMin]=min(A); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax FramePointMax=A(LineMax(2),:); FramePointMin=A(LineMin(2),:); FramePointsPositive(k,:)=FramePointMax; %骨架点储存在FramePoints矩阵中 FramePointsNegative(k,:)=FramePointMin; elseif outdata(LineNum,1)>=0 % 仅有正半周，仅有一个正向最大值 [ColMaxValue,LineMax]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax FramePointMax=A(LineMax(2),:); FramePointsPositive(k,:)=FramePointMax; %骨架点储存在FramePoints矩阵中 else % 有负半周，既有最大值，又有最小值 [ColMaxValue,LineMax]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax [ColMinValue,LineMin]=min(A); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax FramePointMax=A(LineMax(2),:); FramePointMin=A(LineMin(2),:); FramePointsPositive(k,:)=FramePointMax; %骨架点储存在FramePoints矩阵中 FramePointsNegative(k,:)=FramePointMin; end end end %% %滞回耗能 for k=1:LoopNum S(k,1)=polyarea(LoopCircles{k}(:,1), LoopCircles{k}(:,2)); Es(k,1)=FramePointsPositive(k,1)*FramePointsPositive(k,2); kesai(k,1)=S(k,1)/Es(k,1)/(2*pi); end CEs=sum(S);%累计滞回耗能 %% %极限强度、屈服强度、破坏强度,延性系数的计算 [ultimate_load1,LineMax1]=max(FramePointsPositive(:,2));%正向极限强度 yeild_load=ultimate_load1*0.8; [error_yeild_load1,LineMin1]=min(abs(FramePointsPositive(:,2)-yeild_load)); yeild_load1=FramePointsPositive(LineMin1,:);%正向屈服强度及屈服位移 先位移后强度 damaged_load=ultimate_load1*0.85; [error_damaged_load1,LineMin2]=min(abs(FramePointsPositive(:,2)-damaged_load)); damaged_load1=FramePointsPositive(LineMin2,:);%正向破损强度 [ultimate_load2,LineMax2]=min(FramePointsNegative(:,2));%负向极限强度 yeild_load=ultimate_load2*0.8; [error_yeild_load2,LineMin3]=min(abs(FramePointsNegative(:,2)-yeild_load)); yeild_load2=FramePointsNegative(LineMin3,:);%负向屈服强度及屈服位移 damaged_load=ultimate_load2*0.85; [error_damaged_load2,LineMin4]=min(abs(FramePointsNegative(:,2)-damaged_load)); damaged_load2=FramePointsNegative(LineMin4,:);%负向破损强度 [ultimate_displacement1,LineMax1]=max(FramePointsPositive(:,1)); Ductility_coefficient1=damaged_load1(1)/yeild_load1(1);%正向延性系数 [ultimate_displacement2,LineMax2]=min(FramePointsNegative(:,1)); Ductility_coefficient2=damaged_load2(1)/yeild_load2(1);%负向延性系数 Du1=yeild_load1(1)*10;%正向单调加载下的极限位移 Du2=yeild_load2(1)*10;%负向单调加载下的极限位移 %% % % 摘取特征点 % jiazailvli2=[0,4.611429]; % for k=1:LoopNum % A=LoopCircles{k}; % num=length(A); % num1=1; % num2=1; % num3=1; % num4=1; % for ii=1:num % if A(ii,2)>=0 % AFP(num1,:)=A(ii,:); % num1=num1+1; % else % AFN(num2,:)=A(ii,:); % num2=num2+1; % end % if A(ii,1)>=0 % ADP(num3,:)=A(ii,:); % num3=num3+1; % else % ADN(num4,:)=A(ii,:); % num4=num4+1; % end % end % [ColMaxValue,LineMax]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax % [ColMinValue,LineMin]=min(A); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax % FramePointMax=A(LineMax(2),:); % FramePointMin=A(LineMin(2),:); % DATA(1,:)=FramePointMax; % DATA(2,:)=AFN(1,:); % DATA(3,:)=ADN(1,:); % DATA(4,:)=FramePointMin; % DATA(5,:)=AFN(end,:); % DATA(6,:)=ADN(end,:); % jiazailvli(:,:,k)=DATA; % jiazailvli2=[jiazailvli2;DATA]; %end % %强度退化率，刚度退化率 % for k=1:LoopNum/2 % lamda1(k,1)=FramePointsPositive(2*k,2)/FramePointsPositive(2*k-1,2); % lamda2(k,1)=FramePointsNegative(2*k,2)/FramePointsNegative(2*k-1,2); % K1(k,1)=(FramePointsPositive(2*k,2)+FramePointsPositive(2*k-1,2))/(FramePointsPositive(2*k,1)+FramePointsPositive(2*k-1,1)); % K2(k,1)=(FramePointsNegative(2*k,2)+FramePointsNegative(2*k-1,2))/(FramePointsNegative(2*k,1)+FramePointsNegative(2*k-1,1)); % end % %============比较最大位移，删除骨架曲线中同一荷载级的较小骨架点==========% % for k=LoopNum:-1:2 % A=LoopCircles{k}; % B=LoopCircles{k-1}; % [ColMaxValueA,LineMaxA]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax % [ColMinValueA,LineMinA]=min(A); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax % [ColMaxValueB,LineMaxB]=max(B); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax % [ColMinValueB,LineMinB]=min(B); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax % if abs(ColMaxValueA(1)-ColMaxValueB(1))<0.4 % 判断是否为同一荷载级，容差取4mm % if ColMaxValueA(2)<ColMaxValueB(2) % 比较荷载值，删除较小者 % FramePointsPositive(k,:)=[]; % else % FramePointsPositive(k-1,:)=[]; % end % end % if abs(ColMinValueA(1)-ColMinValueB(1))<0.4 % 判断是否为同一荷载级，容差取4mm % if ColMinValueA(2)<ColMinValueB(2) % 比较荷载值，删除较小者 % FramePointsNegative(k-1,:)=[]; % else % FramePointsNegative(k,:)=[]; % end % end % end % 得到加载履历 % MaxD=zeros(LoopNum,1);MinD=zeros(LoopNum,1); % for k=1:LoopNum % A=LoopCircles{k};NumA=length(A); % %======获得两个力的峰值，分别记为MaxForce和MinForce % [ColMaxValue,LineMax]=max(A); % 获得各列最大值ColMaxValue,及各列最大值所在行号LineMax % [ColMinValue,LineMin]=min(A); % 获得各列最小值ColMaxValue,及各列最小值所在行号LineMax % MaxD(k)=A(LineMax(1),1); % MinD(k)=A(LineMin(1),1); % end % DIS=zeros(2*LoopNum+1,1); % for k=2:2:LoopNum*2 % DIS(k)=MaxD(k/2); % DIS(k+1)=MinD(k/2); % end %% figure(1) plot(D,F); xlabel('位移(mm)') ylabel('力(KN)') title('滞回曲线'); grid on hold on figure(2) plot(D(6783:end),F(6783:end)); xlabel('位移(mm)') yl