本文由 徐冬冬 撰寫，黃琳琳 排版。

在本文當中我們將給大家介紹如何去實現隨機疲勞分析當中的雨流計數法。

一：雨流計數法基礎

ASTM E1049-85介紹了疲勞分析當中用到的幾種計數方法（level crossing counting, peak counting, simple-range counting, rangepair counting, rain flow counting），并且論述了每種方法的具體實踐過程。因為雨流計數法得出的結果和實際加載歷程相同，結果也比較接近真實值，故而一般都用雨流計數法。ASTM E1049-85當中的rainflow counting 計出來的是整循環和半循環。這也是最初的雨流計數法。為了避免半循環，我們采用ASTM E1049-85當中的simplified rainflow counting，這種方法的特點是先把信號順序略作改變、使得最大值提前到第一個位置。然后開始循環計數。這樣計出來的循環沒有半循環，只有整循環。

 

信號的預處理包括2部分：

Part I：最大值提前

Part II：峰值谷值交替

黑圈中的點直接去除

下圖展示的是對信號完整的預處理過程：

 紅色信號是原始信號，灰色信號是最大值提前后的信號，藍色信號是提取的峰值谷值交替的信號。

下圖是處理數據的歷程圖

 經過了預處理的算法就可以開始正式的統計。

二、雨流計數法的程序及程序解釋

附件cpp文件當中包括了一個雨流計數法的函數，函數定義如下：

      voidRainFlowCounting(double*A,intL,int N1,int N2,double **R)

   A代表需要分析的應力信號段，L代表應力信號段的長度，N1是均值等分的份數，N2是幅值等分的份數，R用來存儲結果。R是個二維數組，R中存儲的數如下：

         R[0][0]存儲了統計過程中發現的全部載荷歷程的數目。

         R[1][1]存儲了統計歷程中所有載荷歷程的均值。

         R[0][1],R[1][0]全部賦值0，沒有意義。

         R[2][0],R[3][0],R[4][0]…….R[N1+1][0]分別代表每個均值組的下限，

         R[2][1],R[3][1],R[4][1]…….R[N1+1][1]分別代表每個均值組的上限。

        

         R[0][2],R[0][3],R[0][4]…….R[0][N2+1]分別代表每個幅值組的下限，

         R[1][2],R[1][3],R[1][4]…….R[1][N2+1]分別代表每個幅值組的上限。

 

         二維數組當中余下的值是代表當同時滿足

的應力循環的頻次。

         調用過程的具體過程參見cpp文件，以下做簡要說明。注意：頭文件必須包括<vector>,同時紅色標記的3行定義了二維數組。

        

         //從a.dat讀取應力數據，最后分析的結果輸出到b.dat

         ofstreammyfile1("b.dat");

         ifstreammyfile("a.dat");

        

         intL=500000;  //信號長度

         clock_tstart,finish;

         doubletotaltime;

         start=clock();

         double*A=new double[L];

         for(inti=0;i<L;i++)

                   myfile>>A[i];

 

         //First,Movethe signal

         intN1,N2;

         N1=16;N2=16;

         //N1代表均值分的個數

         //N2代表幅值分的個數

 

         double **R=newdouble* [N1+2];  //**R用來存儲結果，并且**R只能按照這一種格式來寫

         for(intk=0;k<N1+2;k++)

         R[k]=new double[N2+2];

 

         //調用雨流計數法的函數。

         RainFlowCounting(A,L,N1,N2,R);

 

         //以下是往外輸出結果的部分。

         cout<<"上兩行是幅值"<<endl;

         cout<<"左邊兩行是均值"<<endl;

         finish=clock();

         totaltime=double(finish-start)/CLOCKS_PER_SEC;

         myfile1<<"%%total time is :"<<totaltime<<endl;

         myfile1<<"%%"<<"第一行是幅值下限，第二行是幅值上限"<<endl;

         myfile1<<"%%"<<"第一列是均值下限，第二列是均值上限"<<endl;

         myfile1<<"%%"<<"第一行第一列的數字代表的是總載荷歷程數"<<endl;

for(int k1=0;k1<=N1+1;k1++)

       {

         for(intk2=0;k2<=N2+1;k2++)

         {

                   cout<<R[k1][k2]<<"   ";

                                     myfile1<<setw(12)<<R[k1][k2]<<"";

         }

                            myfile1<<endl;

         cout<<endl;

       }

 

         輸出結果如下：

三、正確性驗證

Matlab官網提供了RainFlowCounting的程序（評分4.5，評價數53）。利用這個程序的計算結果和C++的結果進行對比。備注: Matlab的結果第一行是幅值，第二行是均值，第三行是循環次數。最后2個0.5次應該合并成1次。C++結果第一列是幅值，第二列是均值，全是整循環，循環次數為1。

 

MATLAB code:

clc

clearall

length=20;

s=10*randn(length,1)+rand(length,1);

s(1)=max(s)+0.1;

s(length)=s(1);

savea.dats-ascii

clc

clearall

tic

loada.dat

tp=sig2ext(a);

rf=rainflow(tp)

toc

 

案例1：

20個點

結果完全一致。

案例2：

30個點

兩者結果完全一致。

來源：CAE仿真大視界公眾號