【有限元分析】风力发电机塔筒筒体校核——matlab雨流计数法

全部程序下载链接
https://download.csdn.net/download/yjw0911/20865581
1、定义sn曲线

%定义SN曲线并求出对应应力幅值的N，包含厚度修正系数及安全系数
          Smax = 500;
%          if bb<=3%DC80  
%            S1 = 29.4723*2.*thicknessSF ;SE = 0.6445*2.*thicknessSF ;
%          else %DC90
           S1 = 33.1563*2.*thicknessSF ;SE = 0.7250*2.*thicknessSF ;
%          end
          N1 = 5E6; NE = 1E15;
%调用雨流统计程序        
          [LoadSpectrum,Marcov,Range,Level,EquivalentLoad]=RainFlow(X,m,bin,T);
          RangeRot = Range';
                 for Tn = 1:bin
                 S = RangeRot(Tn,1);
%调用雨流统计后的应力谱，并求取相应损伤                    
                 LoadSpectrum(Tn,3) = RangeRot(Tn,1);
                       if  S > S1 
                       LoadSpectrum(Tn,2) = 10.^(-3.*log10(S/S1)+log10(N1));
                       else 
                       LoadSpectrum(Tn,2) = 10.^(-5.*log10(S/SE)+log10(NE));
                       end
                   LoadSpectrum(Tn,4) = LoadSpectrum(Tn,1)./LoadSpectrum(Tn,2);
                end
%计算单个工况损伤值
      Data(ii,2) = sum(LoadSpectrum(:,4));
%提取循环次数并计算单个工况损伤值，，对应ii=1:fre
      Data(ii,3) = Data(ii,1).*Data(ii,2);
end
      damage=sum(Data(:,3),'double');

2、Matlab雨流统计方法

function [LoadSpectrum,Marcov,Range,Level,EquivalentLoad]=RainFlow(X,m,bin,T)
%本程序用于进行雨流计数
%作者：文茂诗 日期：2013.11.25
%理论：基于阎楚良,王公权，《雨流计数法及其统计处理程序研究》,农业机械学报，1982年第4期
%处理结果已通过GH Bladed4.2后处理对比验证 详见"Matlab雨流统计处理结果对比2013.11.25.xls"
%[LoadSpectrum,Marcov,Range,Level,EquivalentLoad]=RainFlow(X,m,bin,T)
%输入：X,载荷序列      m,SN曲线的反斜率系数        bin,载荷区间        T，载荷循环次数
%输出：LoadSpectrum,载荷谱        Marcov,马科夫矩阵        Range,载荷幅值
%      Level,载荷均值             EquivalentLoad，当量载荷（1Hz）
[SA,SM]=RainFlowCount(X);
[LoadSpectrum,Marcov,Range,Level,EquivalentLoad]=RainFlowCountStatistics(SA,SM,m,bin,T);

function [SA,SM,X]=RainFlowCount(X)
X=MaxMinX(X);
SA=[];
SM=[];
i=1;
timeuse=300;
ht=tic;
gettime=toc(ht);
    while i<=length(X)-3 && gettime<timeuse
        [SA,SM,X,i,STS]=DoCount(X,SA,SM,i);
        if STS==0
            i=i+1;
        end
        gettime=toc(ht);
    end
index=find(X==max(X));
X=X([index:end 1:index]);
dx=[0;diff(X)];
dy=[dx(2:end);0];
dd=dx.*dy;
X=X(dd<=0);
X(diff(X)==0)=[];
由于篇幅一下省略，最后有下载链接

下载链接
https://download.csdn.net/download/yjw0911/20865581
CAE有限元仿真服务