基于简化雨流计数法处理随机载荷谱
简介
机械结构在服役条件下通常承受振动载荷的作用，振动载荷可以视为一种随机载荷。在对机械结构进行疲劳可靠性分析时，直接使用随机载荷非常不便，因此往往需要将随机载荷通过特定的计数方法转化为若干恒幅载荷，并编制为由若干恒幅载荷谱块构成的变幅载荷谱。恒幅载荷作用下的疲劳寿命估算，可以借助S-N曲线解决。而通过结合Miner线性累积损伤法则，则可应用于变幅载荷下的疲劳寿命预测，最终实现估算实际服役条件下机械结构的疲劳寿命。
简化雨流计数法
将不规则的、随机的载荷-时间历程，转化为一系列循环的方法，即为“循环计数法”。在ASTM E1049-85 Standard Practices for Cycle Counting in FatigueAnalysis中，规定了若干种用于疲劳分析的循环计数方法，其中应用最为广泛的方法为雨流计数法（Rain-flow Counting Method）。相比于其他计数方法，雨流计数方法具有更明确的物理意义，这是因为通过雨流计数法获得的恒幅载荷，恰好对应了在该随机载荷作用下的应力-应变响应中的封闭滞回环，如图1所示。
图1 雨流计数法的物理意义
在图1（a）中，阴影部分的循环为通过雨流计数法计数得到的循环，在图1（b）中，阴影部分为与之对应的封闭滞回环。由于封闭滞回环代表材料产生的疲劳损伤，因此通过雨流计数法获得的恒幅载荷所产生的疲劳损伤与变幅载荷作用下产生的疲劳损伤应该是等效的。
雨流计数法又可细分为三点和四点雨流计数法，这里本文仅讨论三点雨流计数法，也称为简化雨流计数法。简化雨流计数法适用于以典型载荷谱段为基础的重复历程，也就是说，机械构件承受的随机载荷谱可以看为某一典型载荷谱段的不断重复，因此可以只针对该典型载荷谱段进行雨流计数法即可。此外，简化雨流计数法要求将典型载荷谱段的最大峰值或谷值放置在载荷谱段的起始位置，否则将会出现计数错误。
之所以得名雨流计数法，这是因为该方法的计数原理可类比雨滴沿多层屋顶流动的过程，其计数方法如下所示：
（1）从随机载荷谱中选取典型谱段，并且使得典型谱段的最大峰值或谷值位于典型谱段的起始位置，如图2所示。需要注意的是，在从随机载荷谱中选取典型谱段之前，还需要首先去除随机载荷谱中的非峰谷值，使得随机载荷谱完全由峰值和谷值构成。
图2 雨流计数典型段的选取
（2）将典型谱段旋转90°放置，如图3所示。将载荷历程看作多层屋顶，设想有雨滴沿最大峰或谷处开始往下流。若无屋顶阻挡，则雨滴反向，继续流至端点。图中雨滴从A处开始，沿AB流动，至B点后落至CD屋面，继续流至D处；因再无屋顶阻挡，雨滴反向沿DE流动至E处，下落至屋面JA'，至A'处流动结束。所示的雨流路径为ABDEA'。
（3）记下雨滴流过的最大峰、谷值，作为一个循环。图中第一次流经的路径，给出的循环为ADA'。循环的幅值和均值可由图中读出。如ADA'循环的载荷幅值为ΔS=5-(-4)=9，平均载荷Sm=[5+(-4)]/2=0.5。
（4）从载荷历程中删除雨滴流过的部分，对各剩余历程段，重复上述雨流计数。直至再无剩余历程为止。如图所示，第二次雨流得到BCB'和EHE'循环；第三次雨流得到FGF'和IJI'循环；计数完毕。
图3 雨流计数过程
简化雨流计数法计算规则
上述图形方式虽然较为直观，但不便于程序编制。对于存在大量循环数的随机谱段，可采用下述步骤编制为简化雨流计数程序：
（1）读取由最大峰或谷处起止的典型段，按载荷谱顺序输入各峰、谷值，直至数据完毕。
（2）读入下一峰、谷值。若数据完毕，则停止。
（3）若数据点数少于3，则返回步骤（2）；若数据点数大于等于3，则由最后读入的3个峰、谷值，计算幅值X和Y。这三点中，第一点与第二点之差的绝对值为Y；第二点与第三点之差的绝对值为X。
（4）比较X和Y的大小。若X<Y，则返回（2）；若X≥Y，则进行步骤（5）。
（5）将幅值Y记作一个循环，删除与Y相应的峰、谷值，返回步骤（3）。
计算示例
这里，本文通过一个简单示例来展示简化雨流计数法的计数规则。考虑如下表所示的一段谱段，谱段的峰、谷点和对应的峰、谷值列于表中。图4给出了与表格对应的用于雨流计数的载荷谱，为了满足简单雨流计数法的要求，载荷谱的最大峰值位于起始点，并且起始点和终止点的数值相同，保证该载荷谱是可重复的。
图4 用于雨流计数的原始载荷谱
首先读入A，B，C三点，计算变程X和Y，有：
由于X<Y，因此返回步骤（2），需要读取下一峰谷值，即读取D。
此时B，C，D为最后读入的三个峰谷值，由此重新计算X和Y，有：
此时满足X≥Y，进入步骤（5），即将此时的变程Y（B→C）记为一循环，并删除B，C点，返回步骤（3）。
删除B，C点后的载荷谱如下所示。
此时剩余的峰谷点为A和D，数量小于3，因此需读取E，并根据A，D，E计算变程X和Y，即有：
由于满足X<Y，因此进入步骤（2），读取峰谷点F，并根据D，E，F计算变程X和Y，即有：
由于满足X<Y，返回步骤（2），读取峰谷点G，并根据E，F，G计算变程X和Y，即有：
由于满足X<Y，返回步骤（2），读取峰谷点H，并根据F，G，H计算变程X和Y，即有：
由于X≥Y，进入步骤（5），将变程Y（G→F）记为一循环，并删除G，F点，返回步骤（3）。
删除G，F点之后的载荷谱如下所示。
此时剩余的峰谷点为A，D，E，H，由于数据点数大于3，因此基于D，E，H点计算变程X和Y，即有：
由于X<Y，因此返回步骤（2），读取峰谷点I，并根据E，H，I点计算变程X和Y，即有：
由于X<Y，返回步骤（2），读取峰谷点J，并根据H，I，J点计算变程X和Y，即有：
由于X<Y，返回步骤（2），读取峰谷点A'，并根据I，J，A'计算变程X和Y，即有：
由于X≥Y，进入步骤（5），将变程Y（I→J）记为一循环，并删除I和J点，返回步骤（3）。
删除I和J点之后的载荷谱如下所示。
此时剩余的峰谷点为A，D，E，H和A'，数量大于3，因此根据E，H，A'点计算变程X和Y，即有：
由于X≥Y，因此将变程Y（E→H）记为一循环，并删除E和H点，返回步骤（3）。
删除E，H点之后的载荷谱如下所示。
此时剩余的峰谷点为A，D和A'，数据点数大于等于3，根据A，D，A'计算变程X和Y，即有：
此时满足X≥Y，进入步骤（5），将变程Y（A→D）记为一循环，并删除A和D，返回步骤（3）。
此时剩余峰谷点的数量不足3，返回步骤（3），需读取下一峰谷点，但此时数据已读取完毕，计数过程结束。
最终简化雨流计数法的计数结果如下所示。
程序实现
为了实现通过简化雨流计数法处理随机载荷谱，本文基于上述计数原理编制了相应的MATLAB程序。该程序由Rainflow_PreArrange和Rainflow_Count两个文件构成。Rainflow_PreArrange文件用于读取用户指定的载荷序列（以列数据的形式指定），去除数据中的非峰谷值，并调整载荷序列的顺序，使得最大峰值位于载荷序列的起始位置。例如，图5给出了包含非峰谷值的原始载荷谱数据，通过Rainflow_PreArrange程序去除非峰谷值和调整载荷顺序之后的载荷谱如图6和图7所示。
图5 包含非峰谷值的原始载荷谱数据
图6 去除非峰谷值之后的载荷谱数据
图7 调整载荷顺序之后的载荷谱数据（适用于简化雨流计数法）
在获得适用于简化雨流计数法的载荷谱数据之后，即可通过Rainflow_Count程序对其进行简化雨流计数。例如，对于图7中所示的载荷谱数据，Rainflow_Count程序给出的简化雨流计数法结果如下所示：
Range Mean Start value Target value
4.0000 1.0000 -1.0000 3.0000
3.0000 -0.5000 -2.0000 1.0000
7.0000 0.5000 4.0000 -3.0000
9.0000 0.5000 5.0000 -4.0000 事实上，MATLAB内置的信号处理工具箱也内置了rainflow函数用于对载荷谱进行雨流计数。对于上述载荷谱数据（存放于变量Spectrum_Data中），通过下述命令对其进行雨流计数：
[c,hist,edges,rmm,idx] = rainflow(Spectrum_Data);
T = array2table(c,'VariableNames',{'Count','Range','Mean','Start','End'})得到对应的输出结果如下所示：
T =
5×5 table
Count Range Mean Start End
_____ _____ ____ _____ ___
1 4 1 2 3
1 3 -0.5 6 7
1 7 0.5 5 8
0.5 9 0.5 1 4
0.5 9 0.5 4 9 注意，上述结果中的Start和End代表载荷数据的索引值，而非载荷值，Count中的0.5代表半个循环。通过对比由本文程序给出的计数结果可以发现，二者是完全一致的，这证明本文编制的简化雨流计数程序是正确的。
附录
Rainflow_PreArrange程序
%程序用于对输入的载荷谱进行预排序 以满足简化雨流计数法的要求(ASTM E1049-85)
%输入载荷谱数据的要求:(1)数据必须只包含峰谷值 (2)首尾数据的数值必须一致(构成可重复的谱块)
clear,clc
%------------------------------载荷谱数据导入-------------------------------
%导入原始载荷谱数据
Spectrum_Data_Origin=importdata("Spectrum_Data1.txt");
%绘制原始载荷谱数据
plot(Spectrum_Data_Origin,"-o");
title("Original spectrum data");xlabel("Inversal");ylabel("Load");
%--------------------------------------------------------------------------
%-----------------------------载荷谱数据预检查1------------------------------
%统计原始数据点的数量
Origin_Data_Num=size(Spectrum_Data_Origin,1);
%检查数据维度是否符合要求(一维列向量数据)
if size(Spectrum_Data_Origin,2)~=1
error("Error:The input load spectrum data must be in "+...
"the form of a one-dimensional column vector..");
end
%检查首尾数据的数值是否一致(必须构成可重复的谱块)
if Spectrum_Data_Origin(1)~=Spectrum_Data_Origin(end)
error("Error: The initial and final data points differ,"+...
"which obstructs the creation of a repeatable load spectrum block.");
end
%--------------------------------------------------------------------------
%-----------------------------载荷谱数据预检查2------------------------------
%检查数据是否只由峰谷值构成
%查找数据中的所有峰值
[~,locs_pks]=findpeaks(Spectrum_Data_Origin);
%查找数据中的所有谷值(通过取负值来实现)
[~,locs_vls]=findpeaks(-1*Spectrum_Data_Origin);
%合并和排序峰值和谷值的位置(不含首尾位置)
CombinedLocs=sort([locs_pks;locs_vls]);
%合并和排序峰值和谷值的位置(增加首尾位置)
CombinedLocs=[1;CombinedLocs;Origin_Data_Num];
%检查峰值和谷值的位置是否覆盖原始数据的全部位置
if isequal(CombinedLocs,(1:Origin_Data_Num)')
else
warning("Warning: The input load spectrum data includes "+...
"other points besides peaks and valleys.");
end
%重新组合剔除非峰谷值之后的数据
Spectrum_Data_Filtered=Spectrum_Data_Origin(CombinedLocs);
%--------------------------------------------------------------------------
%-------------------------重新绘制预检查之后的数据---------------------------
figure;
plot(Spectrum_Data_Filtered,"-o");
title("Filtered spectrum data");xlabel("Inversal");ylabel("Load");
%--------------------------------------------------------------------------
%------------------------------载荷谱数据重排序------------------------------
%移动载荷谱中最大峰值至起始点 使其满足简化雨流计数法的要求
%统计剔除非峰谷值之后的载荷谱数据点数量
Filtered_Data_Num=size(Spectrum_Data_Filtered,1);
%判断剔除非峰谷值之后的载荷谱数据点数量是否满足要求(不能小于3个点)
if Filtered_Data_Num<3
error("Error:The number of valid data points in "+...
"the input load spectrum is insufficient.");
end
%重新确定剔除非峰谷值后的载荷谱中的峰谷值及对应索引
[pks,locs_pks]=findpeaks(Spectrum_Data_Filtered);
[vls,locs_vls]=findpeaks(-1*Spectrum_Data_Filtered);
vls=-1*vls;
%判断起始点和终止点为峰值还是谷值(上述判断未包含起始和终止点)
if Spectrum_Data_Filtered(1)<Spectrum_Data_Filtered(2)
%起始点和终止点为谷值
vls=[Spectrum_Data_Filtered(1);vls;Spectrum_Data_Filtered(1)];
locs_vls=[1;locs_vls;Filtered_Data_Num];
else
%起始点和终止点为峰值
pks=[Spectrum_Data_Filtered(1);pks;Spectrum_Data_Filtered(1)];
locs_pks=[1;locs_pks;Filtered_Data_Num];
end
%判断最大峰值所在索引
[~,index]=max(pks);
locs_max_pks=locs_pks(index);
%移动载荷谱使其满足简化雨流计数法的要求
%初始化存放重排序之后载荷谱的列向量数组
Spectrum_Data_PreArranged=nan(Filtered_Data_Num,1);
%载入数据至重排序之后的载荷谱
%移动最大峰值之后的数据
Data_Moved=Spectrum_Data_Filtered(locs_max_pks:end);
Data_Moved_Num=size(Data_Moved,1);
Spectrum_Data_PreArranged(1:Data_Moved_Num)=Data_Moved;
%移动最大峰值之前的数据
if locs_max_pks==1
%最大峰值出现在起始点 无需进一步处理
elseif locs_max_pks==2
%最大峰值出现在第二个点 只需在数据末尾补充最大峰值点
Data_Moved=Spectrum_Data_PreArranged(1);
Spectrum_Data_PreArranged(Data_Moved_Num+1)=Data_Moved;
else
%最大峰值未出现在前两个点 将前面的数据补至数据末尾
Data_Moved=Spectrum_Data_Filtered(2:locs_max_pks);
Spectrum_Data_PreArranged(Data_Moved_Num+1:end)=Data_Moved;
end
%--------------------------------------------------------------------------
%----------------------------绘制重排序之后的载荷谱--------------------------
figure;
plot(Spectrum_Data_PreArranged,"-o");
title("Rearranged spectrum data");xlabel("reversal");ylabel("Load");
%输出载荷谱数据
fprintf("Processed load spectrum data:\n");
for i=1:Filtered_Data_Num
fprintf("%-12.4f\n",Spectrum_Data_PreArranged(i));
end
%--------------------------------------------------------------------------
Rainflow_Count程序
%程序用于基于简化雨流计数法(ASTM E1049-85)对载荷谱进行计数
%输出载荷谱必须满足采用简化雨流计数法进行计数的要求(可重复的载荷谱)
clear,clc
%------------------------------载荷谱数据导入-------------------------------
%导入载荷谱数据
Spectrum_Data=importdata("Spectrum_Data2.txt");
%绘制载荷谱数据
plot(Spectrum_Data,"-o");
title("Spectrum data");xlabel("Inversal");ylabel("Load");
%--------------------------------------------------------------------------
%--------------------------------变量初始化---------------------------------
%统计数据点的总数
Spectrum_Data_Num=size(Spectrum_Data,1);
%初始化计数得到的循环总数
Rainflow_Counted_Num=0;
%初始化存放计数后的循环的起始和终止索引的数组(数组维度小于等于数据点总数)
Rainflow_Start_Index=nan(Spectrum_Data_Num,1); %存放起始索引
Rainflow_Target_Index=nan(Spectrum_Data_Num,1); %存放终止索引
%初始化存放计数后的循环的幅值和均值的数组(数组维度小于等于数据点总数)
Rainflow_Amplitude=nan(Spectrum_Data_Num,1); %存放循环幅值
Rainflow_Mean=nan(Spectrum_Data_Num,1); %存放循环均值
%初始化存放计数后的循环的起始和终止值的数组(数组维度小于等于数据点总数)
Rainflow_Start_Val=nan(Spectrum_Data_Num,1); %存放循环起始值
Rainflow_Target_Val=nan(Spectrum_Data_Num,1); %存放循环终止值
%初始化存放临时读取的峰谷点的数组
Temp_PV_Point=nan(Spectrum_Data_Num,1);
%初始化存放临时读取的峰谷点索引值的数组
Temp_PV_Point_Index=nan(Spectrum_Data_Num,1);
%--------------------------------------------------------------------------
%------------------------采用简化雨流计数法进行计数--------------------------
%判断当前载荷谱中的数据点数量是否满足计数要求(数量需大于等于3)
if Spectrum_Data_Num<3
error("Error:Insufficient number of data points "+...
"in the input load spectrum.");
end
%读取载荷谱中的三个初始峰谷点的数值
Temp_PV_Point(1:3)=Spectrum_Data(1:3);
%读取载荷谱中三个初始峰谷点的索引值
Temp_PV_Point_Index(1:3)=[1;2;3];
%记录此时读取的峰谷点的总数
Temp_PV_Point_Num=3;
%设定是否读取下一个峰谷点的标志符
ReadNext=0;
%循环读取其余峰谷点用于计数
while(1)
%判断是否需要读取下一峰谷值(步骤2)
if ReadNext==0
%无需读取下一峰谷值 跳过当前步骤
else
%判断所有峰谷点是否均已读取完毕
%获取当前已读取的峰谷值的最大索引
PV_Max_Index=Temp_PV_Point_Index(Temp_PV_Point_Num);
%判断是否能够读取下一峰谷点
if PV_Max_Index>=Spectrum_Data_Num
break;
end
%读取下一峰谷点
%读取下一峰谷点的数值至数组
Temp_PV_Point(Temp_PV_Point_Num+1)=Spectrum_Data(PV_Max_Index+1);
%读取下一峰谷点的索引至数组
Temp_PV_Point_Index(Temp_PV_Point_Num+1)=PV_Max_Index+1;
%更新存放的峰谷点数量
Temp_PV_Point_Num=Temp_PV_Point_Num+1;
end
%判断峰谷点数量是否满足要求
if Temp_PV_Point_Num<3
continue;
end
%由读取到的最后三个峰谷点计算变程X和Y(步骤3)
X=abs(Temp_PV_Point(Temp_PV_Point_Num)-...
Temp_PV_Point(Temp_PV_Point_Num-1));
Y=abs(Temp_PV_Point(Temp_PV_Point_Num-1)-...
Temp_PV_Point(Temp_PV_Point_Num-2));
%比较变程X和Y的大小(步骤4)
if X<Y
%返回步骤(2) 读取下一峰谷值
ReadNext=1;
else
%记变程Y为一循环 删除与Y对应的峰谷值(步骤5)
%更新雨流计数得到的循环总数
Rainflow_Counted_Num=Rainflow_Counted_Num+1;
%记录变程Y的起始和终止索引
Rainflow_Start_Index(Rainflow_Counted_Num)=...
Temp_PV_Point_Index(Temp_PV_Point_Num-2);
Rainflow_Target_Index(Rainflow_Counted_Num)=...
Temp_PV_Point_Index(Temp_PV_Point_Num-1);
%记录变程Y的起始和终止值
Rainflow_Start_Val(Rainflow_Counted_Num)=...
Temp_PV_Point(Temp_PV_Point_Num-2);
Rainflow_Target_Val(Rainflow_Counted_Num)=...
Temp_PV_Point(Temp_PV_Point_Num-1);
%计算变程Y的幅值和均值
Rainflow_Amplitude(Rainflow_Counted_Num)=...
abs(Rainflow_Start_Val(Rainflow_Counted_Num)-...
Rainflow_Target_Val(Rainflow_Counted_Num));
Rainflow_Mean(Rainflow_Counted_Num)=...
0.5*(Rainflow_Start_Val(Rainflow_Counted_Num)+...
Rainflow_Target_Val(Rainflow_Counted_Num));
%删除与变程Y对应的峰谷值
%复 制数据至临时数组
temp_Array=Temp_PV_Point(1:Temp_PV_Point_Num);
%重置存放峰谷值的数组
Temp_PV_Point(1:Temp_PV_Point_Num)=nan;
%重新存放剔除变程Y的峰谷值至数组
j=1;
for i=1:Temp_PV_Point_Num
if (i==Temp_PV_Point_Num-1)||(i==Temp_PV_Point_Num-2)
%不存入变程Y对应的峰谷值
else
Temp_PV_Point(j)=temp_Array(i);
j=j+1;
end
end
%删除与变程Y对应的索引
% 复 制 数据至临时数组
temp_Array=Temp_PV_Point_Index(1:Temp_PV_Point_Num);
%重置存放峰谷值索引的数组
Temp_PV_Point_Index(1:Temp_PV_Point_Num)=nan;
%重新存放剔除变程Y的峰谷值索引值至数组
j=1;
for i=1:Temp_PV_Point_Num
if (i==Temp_PV_Point_Num-1)||(i==Temp_PV_Point_Num-2)
%不存入变程Y对应的峰谷值
else
Temp_PV_Point_Index(j)=temp_Array(i);
j=j+1;
end
end
%更新数组中存放的峰谷点数量
Temp_PV_Point_Num=Temp_PV_Point_Num-2;
%更新是否读取下一峰谷值的标志符
if Temp_PV_Point_Num<3
%剩余峰谷点数量小于3需读入下一峰谷值
ReadNext=1;
else
ReadNext=0;
end
end
end
%--------------------------------------------------------------------------
%-------------------------输出简化雨流计数法计数结果--------------------------
%输出标题行
fprintf("%-12s %-12s %-12s %-12s\n","Range","Mean","Start value","Target value");
%输出数据
for i=1:Rainflow_Counted_Num
fprintf("%-12.4f %-12.4f %-12.4f %-12.4f\n",Rainflow_Amplitude(i),...
Rainflow_Mean(i),Rainflow_Start_Val(i),Rainflow_Target_Val(i));
end
%--------------------------------------------------------------------------
