骨架曲线的提取
代码网上很多，但是却没有教程，所以代码拿来一跑就出错，为了以后方便复习及后面的师弟师妹好使用，在此做个手把手教程。在此感谢133办公室曾师兄，没有他的指导代码仍是跑不动的。
先上代码，点击matlab中的新建脚本，将此代码复制进去修改即可。
clear all
clc;
hyster=importdata('D:\matlab\bin\xjx.mat');%导入数据
LineNum=size(hyster,1);
EndCircle=1;
%--------------------找各级循环分割点位置----------------------负->正or负->零
for i=3:LineNum
if (hyster(i,1)*hyster(i-1,1)<=0)&&(hyster(i-1,1)<0);%--默认初始加载从正半周开始
EndCirclePoints(EndCircle)=i; %分割点的行号EndCirclePoints
EndCircle=EndCircle+1; %中间分割点个数为EndCircle-1
end;
end;
%--------------------检查最后是否存在不完整滞回环----------------------
if EndCirclePoints(EndCircle-1)<LineNum
LoopNum=length(EndCirclePoints)+1; %LoopNum仅为完整循环的个数，应还有最后一个不完整循环
else
LoopNum=length(EndCirclePoints);
end;
hyster=hyster(1:EndCirclePoints(EndCircle-1),:);
deflection=hyster(:,1);%位移
Force=hyster(:,2);%力
% figure()
% plot(deflection,Force,'r') %%绘图
% hold on
Vmp=max(Force(:,1));%%最大荷载(正向）
Vmn=min(Force(:,1));%%最大荷载(负向）
dmp=deflection(find(hyster(:,2)==Vmp),1);%%最大位移（正向）
dmn=max(deflection(find(hyster(:,2)==Vmn),1));%%最大位移（负向）
% plot(dmp,Vmp,'*',dmn,Vmn,'*')
% hold on
Vcrack_p=0.85*Vmp;%%峰值85%
Vcrack_n=0.85*Vmn;%%峰值85%
% plot([0,40],[Vmp,Vmp],'b',[0,40],[Vp,Vp],'k',[-40,0],[Vmn,Vmn],'b',[-40,0],[Vn,Vn],'k')
Vp=0.7*Vmp;%%峰值70%
Vn=0.7*Vmn;%%峰值70%
%% 正向
%% 找出相交临近点
count1=0
for i=2:length(deflection(:,1))
if Force(i-1)<=Vp && Force(i)>=Vp || Force(i-1)>=Vp && Force(i)<=Vp
count1=count1+1
record_section_range_1(count1,1)=i
end
end
%% 采用线性插值（其实可以直接从图上图）
for i=1:length(record_section_range_1)
point=record_section_range_1(i)
x1=hyster(point-1,1);
y1=hyster(point-1,2);
x2=hyster(point,1);
y2=hyster(point,2);
a=(y2-y1)/(x2-x1);
b=y1-a*x1;
d_section_p(i,:)=(Vp-b)/a;
end
%% 确定割线
d_section_p_min=min(d_section_p);
a_p=Vp/d_section_p_min;
%% 确定屈服位移
dy_p=Vmp/a_p;
count2=0
for i=2:length(deflection(:,1))
if deflection(i-1)<=dy_p && deflection(i)>=dy_p
count2=count2+1
record_section_range_2(count2,1)=i
end
end
%% 采用线性插值（其实可以直接从图上图）
for i=1:length(record_section_range_2)
point=record_section_range_2(i)
x1=hyster(point-1,1);
y1=hyster(point-1,2);
x2=hyster(point,1);
y2=hyster(point,2);
a=(y2-y1)/(x2-x1);
b=y1-a*x1;
V_section_p(i,:)=a*dy_p+b;
end
%% 屈服位移
Vy_p=max(V_section_p);
%% 确定破坏位移
%% 找出相交临近点
count5=0
for i=2:length(deflection(:,1))
if Force(i-1)>=Vcrack_p && Force(i)<=Vcrack_p || Force(i-1)<=Vcrack_p && Force(i)>=Vcrack_p
count5=count5+1
record_section_range_5(count5,1)=i
end
end
%% 采用线性插值（其实可以直接从图上图）
for i=1:length(record_section_range_5)
point=record_section_range_5(i)
x1=hyster(point-1,1);
y1=hyster(point-1,2);
x2=hyster(point,1);
y2=hyster(point,2);
a=(y2-y1)/(x2-x1);
b=y1-a*x1;
d_section_p_crack(i,:)=(Vcrack_p-b)/a;
end
dcrack_p=max(d_section_p_crack);
%% 负向屈服强度和位移
%% 找出相交临近点
count3=0
for i=2:length(deflection(:,1))-1
if Force(i+1)<=Vn && Force(i)>=Vn || Force(i+1)>=Vn && Force(i)<=Vn
count3=count3+1
record_section_range_3(count3,1)=i
end
end
%% 采用线性插值（其实可以直接从图上图）
for i=1:length(record_section_range_3)
point=record_section_range_3(i)
x1=hyster(point+1,1);
y1=hyster(point+1,2);
x2=hyster(point,1);
y2=hyster(point,2);
a=(y2-y1)/(x2-x1);
b=y1-a*x1;
d_section_n(i,:)=(Vn-b)/a;
end
%% 确定割线
d_section_n_max=max(d_section_n);
a_n=Vn/d_section_n_max;
%% 确定屈服位移
dy_n=Vmn/a_n;
count4=0;
for i=2:length(deflection(:,1))-1
if deflection(i+1)<=dy_n && deflection(i)>=dy_n
count4=count4+1
record_section_range_4(count4,1)=i
end
end
%% 采用线性插值（其实可以直接从图上图）
for i=1:length(record_section_range_4)
point=record_section_range_4(i)
x1=hyster(point+1,1);
y1=hyster(point+1,2);
x2=hyster(point,1);
y2=hyster(point,2);
a=(y2-y1)/(x2-x1);
b=y1-a*x1;
V_section_n(i,:)=a*dy_n+b;
end
%% 屈服位移
Vy_n=min(V_section_n);
%% 破裂
%% 找出相交临近点
count6=0
for i=2:length(deflection(:,1))-1
if Force(i+1)>=Vcrack_n && Force(i)<=Vcrack_n || Force(i+1)<=Vcrack_n && Force(i)>=Vcrack_n
count6=count6+1
record_section_range_6(count6,1)=i
end
end
%% 采用线性插值（其实可以直接从图上图）
for i=1:length(record_section_range_6)
point=record_section_range_6(i)
x1=hyster(point+1,1);
y1=hyster(point+1,2);
x2=hyster(point,1);
y2=hyster(point,2);
a=(y2-y1)/(x2-x1);
b=y1-a*x1;
d_section_n_crack(i,:)=(Vcrack_n-b)/a;
end
%% 确定割线
dcrack_n=min(d_section_n_crack);
figure()
plot(deflection,Force,'r') %%绘图
hold on
plot(dmp,Vmp,'*',dy_p,Vy_p,'*', dcrack_p,Vcrack_p,'*',dmn,Vmn,'*',dy_n,Vy_n,'*',dcrack_n,Vcrack_n,'*');
%% 汇总
M=[Vmp dmp Vy_p dy_p Vcrack_p dcrack_p;Vmn dmn Vy_n dy_n Vcrack_n dcrack_n] %%第一行为正向，第二行为负向；第一列为最大荷载值；第二列为最大位移；第三列为屈服荷载；第四列为屈服位移;第五列为破坏荷载；第六列为破坏位移
% save('极限和屈服数据.txt','M','ascii'),
fid=fopen('极限和屈服数据.txt','w');%写入文件路径
[m,n]=size(M);
for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%gn',M(i,j));
else
fprintf(fid,'%gt',M(i,j));
end
end
end
fclose(fid);
以上便是提取骨架曲线以及最大位移、最大承载力、屈服位移、屈服承载力的代码
要使用这段代码，应有以下步骤：
1.得把数据整理好，一列是位移，单位为mm，一列是力，单位为kN
2.将该数据导入matlab。在matlab中点击主页再点击新建变量
将准备好的两列数据保存进去，就像这样
3.记住这个文件的名字和路径，再来代码中修改文件导入的地方。
将第三行代码的路径改为刚保存数据的路径即可
4.点击运行
此时将会弹出滞回曲线的图形，由于软件提取的骨架曲线是会漏点的，故只有手动提取。
点击最大的点，位移和力的大小都会出现，记录在表格，骨架曲线最后单独画即可。
耗能能力的提取
耗能能力就是滞回曲线每一圈所包络的面积，这个如果自己处理才是无从下手，这个才是我们代码省时省力的地方！
上代码
% 导入文件并设定基本参数
Data=importdata('D:\matlab\bin\xjx.mat')
datanum = length(Data);
RervP(1) = 1;k=2;
% 识别每圈滞回曲线与x轴的交点
for j = 2 :datanum
if Data(j,1)*Data(j-1,1)<=0
RervP(k) = j;
k = k +1;
end
end
% 对判断点进行修正（对于每圈数据少于5个的滞回圈予以删除）
m=1;
for k = 2 : length(RervP)
if RervP(k) - RervP(k-1) > 5
RervPM(m) = RervP(k-1);
m = m+1;
end
end
RervPM(m) = RervP(k);
% 对滞回曲线进行分解，将每一圈的数据存储到数组 Loop{} 中
LoopID = RervPM(1:2:end);
LoopNum = length(LoopID);
Loop = cell(LoopNum,1);
for j = 1 : LoopNum
if j == LoopNum
Loop{j} = Data(LoopID(j):end,:);
else
Loop{j} = Data(LoopID(j):LoopID(j+1),:);
end
end
% 计算滞回曲线叠加耗能(存入EnerSum（）数组中。
Dispmax = zeros(LoopNum+1,1); EnerSum = zeros(LoopNum+1,1);
for j = 1: LoopNum
x = Loop{j}(:,1);
y = Loop{j}(:,2);
if length(x) == 1
EnerSum(j+1) = [];
Dispmax(j+1) = [];
else
EnerSum(j+1) = trapz(x,y) + EnerSum(j);
Dispmax(j+1) = max(x);
end
end
使用方法与提取骨架曲线的类似，主要是运行后去哪里找数据
运行完毕后EnerSum（）数组中便是每一圈滞回的耗能能力。
将里面数据拿到表格做处理即可。