一种基于逆向雨流计数法的载荷谱重建方法

在机械结构的疲劳可靠性研究中，除了考虑如S-N曲线等材料参数的统计特性，结构承受载荷的统计特性也是需要考虑的内容。对于处于服役条件下的机械结构，其往往承受随机振动载荷的作用。为评估承受随机载荷作用下机械结构的疲劳可靠性，一种较为可行的方法为基于随机振动载荷的统计特性，生成若干组不同的随机载荷，然后通过蒙特卡洛法完成随机载荷作用下机械结构的疲劳分析，实现对机械结构的疲劳可靠性评估。上述方法依赖于建立一种随机载荷谱的重建方法。

考虑到雨流计数法在随机载荷的处理中具有明确的物理意义，本文旨在建立一种基于逆向雨流计数法的载荷谱重建方法，其基本思路为对于机械结构实测的随机载荷进行雨流计数，然后基于雨流计数结果来重建随机载荷，使得重建之后的随机载荷谱具有与实测随机载荷谱相同的雨流计数结果，但载荷序列完全不同。通过使用重建之后的多种随机载荷谱进行疲劳分析，可充分考虑随机载荷的不确定对机械结构疲劳可靠性的影响。

本文建立的载荷谱重建方法可视为雨流计数法的逆过程，即根据雨流计数获得的恒幅载荷，不断向部分重建的载荷谱中插入这些恒幅载荷，最终雨流计数获得的恒幅载荷全部插入到部分重建的随机载荷谱中，构成一段完整的随机载荷谱。为了保证重建之后的随机载荷谱也能获得与原载荷谱相同的雨流计数结果，在插入过程中必须满足特定的规则。

首先，我们需要根据雨流计数结果中的幅值和均值来构成一个二维雨流矩阵，这个二维雨流矩阵的行索引代表雨流计数获得的恒幅载荷的最大载荷，列索引代表雨流计数获得的恒幅载荷的最小载荷，矩阵元素的值代表雨流计数获得的该恒幅载荷的循环周次。这里，我们在重建过程中不考虑载荷方向的影响（即不区分从最小载荷→最大载荷还是最大载荷→最小载荷），此时雨流矩阵的非零元素将全部位于下三角部分。

在重建随机载荷谱的过程中，雨流矩阵的元素必须依次插入，并且只有当某一元素的全部循环均插入到载荷谱中之后，再开始插入下一个元素。插入的顺序为首先插入雨流矩阵中最左下角的元素，然后插入对应行的全部列的元素。然而，再插入上一行的全部元素，以此类推，直到所有元素均插入到重建的载荷谱中。在将一个载荷循环插入部分重建的载荷历程中时，只需让插入循环的峰值小于接收循环的峰值，插入循环的谷值大于接收循环的谷值即可。在插入过程中，插入循环可随机插入到任意满足上述条件的接收循环中。

下面，通过一个简单示例来解释上述重建过程。考虑如图1所示的一个4×4的二维雨流矩阵。

图1 二维雨流矩阵示例

基于上述原理，首先插入位于左下角的最大循环，即4-1-4，插入之后的结果如下图所示。

然后插入该行中其余列的循环，即插入4-2-4的2个循环，这2个循环既可以插入左侧，也可以插入右侧，这里选择将2个循环分别插入两侧，插入之后的结果如下图所示。

然后继续插入4-3-4的2个循环，不难看出，由于该循环的峰值均等于当前序列的峰值，谷值均大于当前序列的所有谷值，因此可以插入到任意位置。这里选择将2个循环同时插入到第2个循环的位置，插入之后的结果如下图所示。

在此之后，继续插入第3行的数据，首先为3-1-3（或者1-3-1）的2个循环，根据上面的规则，这2个循环只能插入到4-1-4循环中，这里选择将其插入到4-1-4循环的两侧，插入之后的结果如下图所示。

然后插入3-2-3（或者2-3-2）的1个循环，这里可以有多种插入方法，我们选择插入在第1个循环的位置，插入之后的结果如下图所示。

最后插入下一行的2-1-2，共有3个循环，这里选择将其同时插入到4-1-3的循环处，插入后的结果如下图所示。

将上述循环进行雨流计数，得到雨流计数的结果如下表所示。

可以发现，该计数结果与上述雨流矩阵中的结果完全一致。

基于上述原理，笔者编制了相应的载荷谱重建程序。该程序由主程序Rainflow_Reconstruction和函数Rainflow_Load_Insert构成。函数Rainflow_Load_Insert用于将某一段恒幅载荷按上述规则插入到部分重建的随机载荷谱中，并返回一个新的随机载荷谱。主程序用于读取雨流计数获得的幅值、均值及循环周次，并生成对应的二维雨流矩阵，然后主程序将遍历二维雨流矩阵中的元素，并调用函数Rainflow_Load_Insert，实现循环的插入过程，最终形成完整的重建载荷谱。完整的程序参见附录。

主程序在运行时只需要读取包含雨流计数结果的文本文件（Rainflow_Rst.txt），雨流计数结果的排列格式如下所示。

1   1   2.5
2   2   3
1   3   2.5

这里，第一列数据代表雨流计数获得的每一个恒幅载荷的循环周次，第二列代表对应的载荷幅值，第三列代表对应的载荷均值。

下面，本文通过一个简单的示例来验证程序的正确性。考虑如图2所示的一段载荷谱，下面需要对该段载荷谱进行雨流计数，并根据雨流计数结果重建载荷谱，要求重建之后的载荷谱必须与原始载荷谱具有相同的雨流计数结果。

图2 原始载荷谱数据

首先本文基于《基于简化雨流计数法处理随机载荷谱》中给出的程序对图2中的原始载荷谱数据进行处理，以剔除数据中的非峰谷值，并将载荷谱中的最大峰值移动到载荷谱的起始位置。处理完成之后的载荷谱如图3所示。

图3 剔除非峰谷值并重排序之后的原始载荷谱数据

对图3中的载荷谱数据进行雨流计数，计数结果如图4所示，由于雨流计数结果中包含了非常多的载荷循环，因此这里通过图形的形式给出每个循环的幅值和均值。

图4 雨流计数结果

基于图4中的雨流计数结果，采用本文建立的方法进行载荷重建，得到重建之后的载荷谱如图5所示。可以看到，由于每次运行程序时载荷序列是被随机插入到部分重建的载荷谱中的，因此可以获得不同的载荷序列。

图5 重建之后的载荷谱数据

如果分别对图5中的三组重建的载荷谱进行雨流计数，可以得到其雨流计数结果如图6所示。图中还给出了原始载荷谱的雨流计数结果作为对比。

图6 原始载荷谱和重建载荷谱的雨流计数结果对比

从图6中可以看到，原始载荷谱与重建载荷谱的雨流计数结果完全一致，但这些载荷谱具有完全不同的载荷序列。如果我们认为载荷谱的雨流计数结果从疲劳角度反映了载荷的变化规律，则基于本文重建的载荷谱，可以充分考虑载荷谱的随机性对疲劳寿命的影响。

%程序用于基于雨流计数结果重建载荷谱
%雨流计数结果以循环周次 幅值 均值的形式指定(不考虑循环方向)
clear,clc

%----------------------------读取雨流计数结果--------------------------------
%读取存放雨流计数结果的数据文件
Rainflow_Rst_Data=importdata("Rainflow_Rst.txt");
%检查数据文件有效性
%检查第一列的循环周次是否为正整数
if all(Rainflow_Rst_Data(:,1)>0)
else
    error("Error:The cycle count must be a positive integer.");
end
if all(Rainflow_Rst_Data(:,1)==floor(Rainflow_Rst_Data(:,1)))
else
    error("Error:The cycle count must be a positive integer.");
end

%根据雨流计数的幅值和均值计算最大值和最小值(分别存放至雨流计数结果的第4列和第5列)
%读取幅值和均值
Delta_S=Rainflow_Rst_Data(:,2);
Fm=Rainflow_Rst_Data(:,3);
%计算最大值和最小值
Fmax=(2*Fm+Delta_S)/2;
Fmin=(2*Fm-Delta_S)/2;
%存放最大值和最小值至雨流计数结果
Rainflow_Rst_Data(:,4)=Fmax;
Rainflow_Rst_Data(:,5)=Fmin;
%获取雨流计数结果的数据行数
Rainflow_Rst_Data_Row_Num=size(Rainflow_Rst_Data,1);
%--------------------------------------------------------------------------

%----------------------基于雨流计数结果构造二维雨流矩阵-----------------------
%统计雨流计数结果中出现的载荷值
Load_List=unique([Rainflow_Rst_Data(:,4);Rainflow_Rst_Data(:,5)]);
%计算雨流计数结果中出现的载荷数量
Load_Num=size(Load_List,1);
%初始化二维雨流计数矩阵
Rainflow_Matrix=zeros(Load_Num,Load_Num);
%逐循环写入二维雨流矩阵元素
for i=1:Rainflow_Rst_Data_Row_Num
    %读取当前行元素的最大值和最小值
    Fmax_Val=Rainflow_Rst_Data(i,4);
    Fmin_Val=Rainflow_Rst_Data(i,5);
    %根据最大值确定雨流矩阵的行索引
    Rainflow_Matrix_Row_index=find(Load_List==Fmax_Val);
    %根据最小值确定雨流矩阵的列索引
    Rainflow_Matrix_Column_index=find(Load_List==Fmin_Val);
    %根据行索引和列索引写入二维雨流矩阵元素(计算得到的循环周次)
    Rainflow_Matrix(Rainflow_Matrix_Row_index,Rainflow_Matrix_Column_index)=...
       Rainflow_Rst_Data(i,1); 
end

%以可视化的方式展示二维雨流矩阵
%创建表格的行标题和列标题
UITable_RowNames=num2str(Load_List);
UITable_ColumnNames=num2str(Load_List);
%创建表格用于显示二维雨流矩阵
%Table_Hanlde=figure;
uitable('Data',Rainflow_Matrix,...
    'RowName',UITable_RowNames,'ColumnName',UITable_ColumnNames);
%--------------------------------------------------------------------------

%-------------------------基于二维雨流矩阵重建载荷谱--------------------------
%初始化判断是否为首个元素的标志符
Is_First_Element=1;
%遍历二维雨流矩阵中的所有元素
for Row_Index=Load_Num:-1:1
    for Column_Index=1:Load_Num
        %判断矩阵中的当前元素是否存在循环周次
        if Rainflow_Matrix(Row_Index,Column_Index)==0
            %跳过当前循环 继续执行下一元素的判断
            continue;
        end

        %判断当前元素是否为首个被 插入载荷历程的元素
        if Is_First_Element==1
            %基于首个元素构建初始重建载荷谱
            %获取首个元素的最大载荷(行索引对应的载荷)
            Fmax_Val=Load_List(Row_Index);
            %获取首个元素的最小载荷(列索引对应的载荷)
            Fmin_Val=Load_List(Column_Index);

            %创建载荷历程(首个插入元素不考虑插入位置)
            %获取待插入的循环周次数量
            Inserted_Cycle_Num=Rainflow_Matrix(Row_Index,Column_Index);
            %初始化存放重建载荷历程的数组
            Reconstructed_Load_Hist=nan(2*Inserted_Cycle_Num+1,1);
            %插入元素至存放重建载荷历程的数组
            Reconstructed_Load_Hist(1:2*Inserted_Cycle_Num,1)=...
                repmat([Fmax_Val;Fmin_Val],Inserted_Cycle_Num,1);
            Reconstructed_Load_Hist(end,1)=Fmax_Val;
            %更新代表首个被 插入元素的标志符
            Is_First_Element=0;

            %跳过下面针对非首个插入元素执行的流程
            continue;
        end

        %针对非首个插入元素执行常规插入流程以构成重建载荷谱
        %获取待插入的循环周次数量
        Inserted_Cycle_Num=Rainflow_Matrix(Row_Index,Column_Index);
        %获取待插入循环的最大载荷(行索引对应的载荷)
        Fmax_Val=Load_List(Row_Index);
        %获取待插入循环的最小载荷(列索引对应的载荷)
        Fmin_Val=Load_List(Column_Index);
        %逐循环插入以构成重建载荷谱
        for i=1:Inserted_Cycle_Num
            Origin_Load_Hist=Reconstructed_Load_Hist;
            [Reconstructed_Load_Hist]=Rainflow_Load_Insert(Origin_Load_Hist,...
                Fmax_Val,Fmin_Val,1);
        end

    end
end
%--------------------------------------------------------------------------

%----------------------绘制基于二维雨流矩阵重建的载荷谱-----------------------
figure;plot(Reconstructed_Load_Hist,"-o");
title("Reconstructed load history");xlabel("Inversals");ylabel("Load");
%--------------------------------------------------------------------------