做疲劳计算的同行应该都有体会拿到一段实测载荷谱后最头疼的不是算应力大小而是怎么把一条杂乱无章的时间-载荷曲线整理成能直接喂给Miner疲劳损伤公式的“循环列表”。雨流计数法就是这个环节的行业标准而它的MATLAB实现里围绕三点法和四点法的选择一直有不少讨论。这里我把两种思路的完整代码、推导逻辑、验证过程和踩坑点一起整理出来希望对正在做结构疲劳分析、有限元后处理或者寿命预测的朋友有些帮助。先说结论三点法和四点法本质上是同一套雨流思想的两种表达方式。三点法更贴近“雨流”的物理图像适合理解原理四点法更接近工程标准ASTM E1049的推荐流程适合直接输出全循环。下文我会用同一个算例把两种方法跑通并对比两者的输出差异最后给出一个可以直接抄作业的MATLAB函数封装。1. 雨流计数法到底在数什么迟滞回线与循环统计1.1 为什么不能直接“数峰谷”很多人第一次接触疲劳分析时第一反应是把载荷时间历程里的波峰波谷数一遍认为一个峰到相邻谷就算一个循环。这个直觉在等幅正弦载荷下没错但用在随机载荷谱上就会出大问题。看一个简单的例子载荷序列从0拉升到10然后回落到-5再升到10最后回到0。如果只数相邻峰谷会得到三段小幅变程0→10、10→-5、-5→10然后剩下一个10→0的半段。这会造成两个后果一是把本来连续的应力路径割裂了二是忽略了一个跨越整段的大循环——材料实际上经历了一次从-5到10的大幅往返这个大幅度循环对疲劳损伤的贡献远大于其中任何一个小循环。雨流计数法之所以叫“雨流”是因为它把载荷历史想象成一系列屋顶雨滴从每个峰顶或谷底开始向下流只有形成封闭回线的那部分才被计为一个循环。物理上对应的是材料的迟滞回线只有闭合的应力-应变回线才产生疲劳损伤。所以雨流法数出来的循环不是简单的时间先后顺序而是材料应力-应变响应的闭合顺序。1.2 半循环与全循环三点法和四点法的本质在雨流计数体系里一个完整的循环由两个半循环组成一个上升半循环和一个下降半循环。当两个半循环的变程range和均值mean相同时它们在材料层面就会闭合出一个完整的迟滞回线可以合并为一个全循环。三点法的思路是每次取连续的三个峰谷点A、B、C判断A→B这段变程能否被提取为一个半循环。如果A→B的跨度小于等于B→C的跨度说明A→B这段回线已经闭合就记录一个半循环然后删掉A点回溯一格继续判断。三点法的名字就来自“每次看三个点”这个操作。四点法来自ASTM E1049标准的实现思路每次取连续的四个峰谷点A、B、C、D比较中间那段变程X|B−C|与两侧变程Y|A−B|、Z|C−D|的大小。只有当X同时小于等于Y和Z时B→C才是一个独立闭合的循环直接记录为一个全循环并删除B、C两个点。四点法直接输出全循环不需要事后做半循环配对。两者的计数结果在理想数据上是一致的区别在实现粒度和工程便利性。2. 三点法用三个点做半循环计数的逻辑2.1 三点判据的来历三点法的核心判据就一句话当峰谷序列中连续三点A、B、C满足|A−B| ≤ |B−C|时A→B这段可以单独提取为半循环。这个判据的物理含义值得琢磨。峰谷序列是交替变化的A到B是一个完整的上升或下降沿B到C是紧接着的反向沿。如果A→B的跨度比B→C小那么A→B这段折返已经被后续的B→C“覆盖”住了它在材料内部会先于外层大回线闭合。雨流法有一个基本原则先闭合的小回线先计数大回线跨过它们继续延伸。用数字说话。假设A0、B5、C-3那么A→B变程为5B→C变程为8。5≤8成立说明0→5这段上升沿是一个可以闭合的半循环。但如果A0、B5、C2虽然这不是标准的峰谷交替但用于理解判据则5≤3不成立说明0→5这段上升沿还没有闭合需要继续往后找更大的反向回线。Downing和Socie在1982年提出的简化算法就是基于这个判据它把原本需要反复扫描的雨流过程变成了一个线性扫描过程满足条件就提取并删除不满足就滑向下一个窗口。这也是目前多数教学版MATLAB代码的底子。2.2 MATLAB核心代码与逐行说明下面这段代码是我在实际项目里用的三点法版本为了可读性我保留了数组删除操作function cycles rainflow_3point(ext) % 输入: ext 峰谷值序列, 列向量 % 输出: cycles [range, mean, count] % count 为 0.5 表示半循环, 1.0 表示配对后的全循环 ext ext(:); n length(ext); cycles zeros(n, 3); nc 0; i 1; while i length(ext) - 2 A ext(i); B ext(i 1); C ext(i 2); rangeAB abs(B - A); rangeBC abs(C - B); if rangeAB rangeBC % 提取 A-B 半循环, 方向由 B-A 的符号决定 nc nc 1; if B A direction 1; % 上升半循环 else direction -1; % 下降半循环 end cycles(nc, :) [rangeAB, (A B) / 2, 0.5]; % 删除A点, 指针回退一格 ext(i) []; if i 1 i i - 1; end else % 不满足条件, 窗口右移 i i 1; end end % 剩余不足三点的序列, 逐段按半循环记录 for k 1:length(ext) - 1 nc nc 1; cycles(nc, :) [abs(ext(k 1) - ext(k)), ... (ext(k) ext(k 1)) / 2, 0.5]; end cycles cycles(1:nc, :); end有一行很多人容易写错删除A点后指针要回退一格。因为删除后原来的B点移到了A的位置而B与它前面的点之间可能又满足提取条件了必须重新检查一次。如果不回退就会漏计循环我最早写的时候就是漏了这个回溯导致结果比标准实现少了整整一个半循环。2.3 半循环配对策略三点法输出的是一堆半循环它们需要按“方向相反、变程和均值相等”的原则配成对。配对之后才是一般疲劳计算里用的全循环。function fullCycles pairHalfCycles(halfCycles) % halfCycles: [range, mean, count, direction] count size(halfCycles, 1); used false(count, 1); fullCycles []; for i 1:count if used(i) continue; end for j i 1:count if used(j) continue; end if halfCycles(i, 4) * halfCycles(j, 4) 0 ... % 方向相反 abs(halfCycles(i, 1) - halfCycles(j, 1)) 1e-8 ... % 变程相等 abs(halfCycles(i, 2) - halfCycles(j, 2)) 1e-8 % 均值相等 fullCycles(end 1, :) [halfCycles(i, 1), halfCycles(i, 2), 1.0]; used(i) true; used(j) true; break; end end end % 剩余未配对的半循环保留为0.5个循环 for i 1:count if ~used(i) fullCycles(end 1, :) halfCycles(i, 1:3); end end end关于浮点容差建议用变程量级的1e-8而不是固定值。如果载荷谱数值很大比如上千兆帕固定容差1e-8会无法配对如果数值很小固定容差又可能误配。用max(半循环变程) * 1e-8做容差更稳。3. 四点法直接输出全循环的工程实现3.1 四点法的判定含义四点法每次看四个点A、B、C、D设X |B − C|中间变程Y |A − B|左侧变程Z |C − D|右侧变程判据是当X ≤ Y 且 X ≤ Z时B→C构成一个完整循环。这个判据的含义用一句话概括中间的折返是最小的它被两侧的更大变程包裹住了。只有这样的循环才是“已经闭合”的迟滞回线可以直接从载荷路径里摘除。相比之下三点法每次提取的是一个半循环一个半循环到全循环还要经过配对四点法则通过同时看左右两侧的变程确认了中间这段回线已经闭合所以一步到位输出全循环。这也是为什么ASTM E1049的推荐实现里四点是真正的主流。工程上拿到雨流计数结果就是为了直接算损伤谁都不想在循环配对这件事上多花精力。3.2 MATLAB核心代码删除与重组function cycles rainflow_4point(ext) % 输入: ext 峰谷值序列, 列向量 % 输出: cycles [range, mean, count], count 始终为 1.0 或 0.5 ext ext(:); n length(ext); cycles zeros(n, 3); nc 0; while length(ext) 4 A ext(1); B ext(2); C ext(3); D ext(4); X abs(B - C); Y abs(A - B); Z abs(C - D); if X Y X Z % 提取 B-C 全循环 nc nc 1; cycles(nc, :) [X, (B C) / 2, 1.0]; % 删除 B 和 C ext [ext(1); ext(4:end)]; else % 删除 A, 窗口整体右移 ext ext(2:end); end end % 剩余不足四点的序列, 按半循环逐段记录 for k 1:length(ext) - 1 nc nc 1; cycles(nc, :) [abs(ext(k 1) - ext(k)), ... (ext(k) ext(k 1)) / 2, 0.5]; end cycles cycles(1:nc, :); end注意一个细节提取B→C循环后删除的是B和C两个点原来的A和D留了下来A与D变成相邻点。下一次循环会重新组合A、D、后续点判断它们之间是否还能形成循环。这一点与三点法不同——三点法删除A后需要回退指针四点法删除B、C后本来就应该重新组合相邻关系所以不需要额外回溯逻辑这也是四点法代码更简洁的原因。另一个细节是为什么删除B、C后要保留A而不是也删掉A和D之间可能存在尚未闭合的更大回线的一部分A作为这根大回线的端点必须继续参与后续判断。删多了会把大循环的信息丢掉。3.3 三点法配完对和四点法结果真的等价吗我在同一组数据上做过对比结论是在理想情况下等价但工程实现里有差异。三点法加配对等于先把所有可闭合的小回线都以半循环形式抽出来再寻找配对。问题在于配对阶段如果存在两个以上变程、均值相同的半循环配对顺序会影响最终结果是0.5还是1.0。而四点法通过“摘除闭合循环”的方式直接给出了全循环避免了半循环配对时的归属歧义。实际载荷谱很少出现完全相同的半循环所以差异通常不影响最终损伤值。但如果你的载荷谱包含大量对称等幅段建议直接用四点法结果更稳定。这也是我在经历过一次配对歧义问题后把主计算函数换成了四点法的原因。4. 实测验证一个手算例子全流程复现4.1 数据准备与峰谷提取先准备一个代表性测试序列手工也能算的那种s [0, 5, -3, 4, -2, 6, -4, 2, -5, 3, 0];峰谷提取函数如下function ext extractPeaks(series) % 去除相邻重复值并提取峰谷序列 series series(:); d diff(series); keep find(d ~ 0); if isempty(keep) ext series(1); return; end series series([1; keep 1]); dd diff(series); signs sign(dd); turns find(diff(signs) ~ 0) 1; ext series([1; turns; end]); end这组数据本身已经是峰谷交替所以extractPeaks不会改变它。这里要提醒一句如果原始数据是等间隔采样的连续信号必须先用这个函数做数据压缩否则雨流算法会在单调变化的中间点上做无意义判断既慢又会产出错误循环。4.2 三点法与四点法的手算对照我对这组数据做了手算推演。用四点法完整过程如下表步骤当前前四点判断动作10, 5, -3, 4X8, Y5, Z7XY丢弃025, -3, 4, -2X7, Y8, Z6XZ丢弃53-3, 4, -2, 6X6, Y7, Z8X≤Y且X≤Z记录循环(4,-2)range6mean1删除4和-24-3, 6, -4, 2X10, Y9, Z6XY丢弃-356, -4, 2, -5X6, Y10, Z7X≤Y且X≤Z记录循环(-4,2)range6mean-1删除-4和266, -5, 3, 0X8, Y11, Z3XZ丢弃67剩余-5, 3, 0无法形成四点记录(-5,3)半循环range8 mean-1记录(3,0)半循环range3 mean1.5最终输出两个全循环和两个半循环。用三点法加配对跑同样数据最后配对出来的也是两个range6的全循环配对不上的是range8和range3两个半循环。两者完全一致这个例子可以作为验证代码的基准测试用例。4.3 交叉验证与成熟参考实现的对比代码写完后一定要做交叉验证。我在File Exchange上找到了Adam Nieslony公开的Rainflow Counting Algorithm把他的输出和自己的四点法结果做了统计对比。方法是把结果矩阵按range和mean聚类再对比循环计数两者完全吻合时才算通过验证。这里分享一个验证技巧不要只对比最终损伤值损伤值只是把循环代入S-N曲线后的加权求和两个错误的实现也可能算出接近的损伤。要对比循环计数矩阵也就是range、mean的分布才能确认雨流计数本身是对的。5. 实战中的那些坑边界条件、起点处理和性能5.1 首尾不闭合问题与起点选择雨流法假设载荷谱是闭合回线但实际测量的时间序列是断开的。如果直接拿原始序列跑算法第一个上升沿和最后一个下降沿会因为“缺少另一侧的回线”而无法配对可能漏掉一个大幅循环。标准做法是先把峰谷序列做环形重排让绝对值最大的极值点最大峰或最小谷变成序列起点。这样最外层的大回线以极值点为中心展开首尾连接时自然闭合。重排操作在MATLAB里很简单[~, maxIdx] max(abs(ext)); ext circshift(ext, -(maxIdx - 1));但注意重排后仍然需要把残余的尾段按半循环处理因为有些回线本身就是不闭合的。工程上对半循环的处理有两种习惯一种是保留0.5次循环计入损伤另一种是把它翻折补全成一个完整循环。保守起见我按0.5循环计算因为半循环代表材料经历了一次完整的单程应变确实产生了接近半个循环的损伤。5.2 浮点误差与重复值处理实测数据的浮点噪声会导致三种典型问题等值平台、接近等值的抖动、以及变程比较时差一点点导致的误判。等值平台的处理已经在extractPeaks里做了用diff ~ 0把相邻重复压缩成一个点。这里有个细节平台值的保留位置应该取平台最后一个值因为平台终点才是极值点位置。用series series([1; keep 1])得到的恰好是每个变化段的终点符合要求。接近等值的抖动如果不去除雨流算法会数出一堆幅值极小的无效循环。我的做法是在提取峰谷前先做一个幅度滤波把与前一极值差小于tol的点丢弃。这个tol通常取载荷满量程的1%即可或者根据传感器精度来定。但要注意这个滤波必须在提取峰谷之后再做否则会破坏峰谷交替的序列结构。变程比较的浮点容差也很关键。三点法判据rangeAB rangeBC在两端几乎相等时会受浮点误差影响导致一次误判。建议改成tol max(abs(ext)) * 1e-12; if rangeAB rangeBC tol5.3 几万点长序列的性能优化前面给出的代码为了可读性用了动态数组删除这个写法在MATLAB里性能很差。一万点的峰谷序列还好几十万点就会慢到难以接受。原因是ext(i) []每次删除都会触发数组整体搬移复杂度是O(n²)。优化方案有两种。第一种是用逻辑索引分批删除但雨流算法内部删除和判断交错进行批量删除并不好写。第二种是用一个“栈”结构只在数组末尾做删除和追加操作。四点法天然适合栈实现每次从顶部取四个点满足条件就弹出B、C不满足就弹出A。MATLAB里可以用ext(1) []; ext(2) [];这样的前删操作性能同样一般我用的是索引指针加预分配数组的方式核心思想是先记录循环和存活点索引最后统一重建数组。内存不是主要问题循环次数才是。十万点随机载荷谱四点法的判断次数大约在二十万次以内MATLAB本身跑得动。如果要用于实时处理或者嵌入式验证建议把核心循环写成MEX函数或C代码。我在项目中实测纯MATLAB四点法处理十万点数据大约零点几秒MEX版本可以到几十毫秒看你项目需求决定。6. 从循环统计到寿命估算雨流结果的使用姿势6.1 用Miner线性损伤准则计算损伤雨流计数的最终目的通常是算疲劳损伤。Miner法则说损伤量是每个循环损伤的线性叠加$$D \sum_{i1}^{N} \frac{n_i}{N_i}$$其中$n_i$是雨流计数得到的循环次数$N_i$是第$i$个循环对应应力幅值下的疲劳寿命通常由S-N曲线给出。工程近似常用$$N_i C \cdot S_i^{-m}$$其中$S_i$是应力幅值C和m是材料参数。MATLAB实现function D minerDamage(cycles, C, m) % cycles: [range, mean, count], 这里用 range 代替应力幅 S % 注意: 如果 S-N 曲线用幅值定义, 应传入 range/2 S cycles(:, 1) / 2; N C .* S .^ (-m); D sum(cycles(:, 3) ./ N); end这里有一个容易踩的坑S-N曲线里的应力幅值到底是range还是range/2不同材料手册习惯不同。ASTM文档里通常用应力范围ΔS有些书用应力幅值Sa。用错一个倍率寿命预测可能差出十来倍。接入项目前先把手册里的公式定义搞清楚。6.2 生成雨流矩阵用于后处理除了损伤值工程上还常把循环结果整理成range-mean矩阵也称为雨流矩阵。它直观展示了载荷谱在不同幅值和均值区间的分布是后续频域疲劳分析的基础。% cyclesAll: 所有循环(含半循环补全后)的 [range, mean, count] rangeBins 0:5:100; meanBins -50:5:50; flowMatrix histcounts2(cyclesAll(:, 2), cyclesAll(:, 1), ... meanBins, rangeBins);注意histcounts2第一个输入对应行方向第二个输入对应列方向参数顺序反了矩阵会转置。生成的矩阵可以用imagesc可视化横轴均值、纵轴变程颜色深浅代表循环次数。6.3 封装成完整函数的建议我在项目里的最终封装是这样的接口function result fatigueCountFromTimeSeries(loadHistory, params) % loadHistory: 原始载荷时间历程 % params: 包含材料常数、滤波容差、计数方式等 ext extractPeaks(loadHistory); ext amplitudeFilter(ext, params.tol); ext reorderToMaxExtreme(ext); if strcmpi(params.method, three) half rainflow_3point(ext); cycles pairHalfCycles(half); else cycles rainflow_4point(ext); end result.cycles cycles; result.damage minerDamage(cycles, params.C, params.m); result.matrix buildRainflowMatrix(cycles, params.bins); end这样的接口好处是换计数方法、换材料参数都只改一处不会牵动整个疲劳分析流程。我现在的疲劳分析工具链基本就是这个结构前处理提取峰谷、雨流计数、损伤计算三层分离方便出问题和同事核对时快速定位。最后分享一个小技巧写完这个模块后我给自己的代码配了一组带标准解的回归测试集包括正弦载荷、方波载荷、随机载荷以及上面那个手算的七点序列。每次改完代码先跑回归确保计数输出不回归。疲劳计算不像数值算法有直接的收敛性判断一个隐蔽的循环计数错误会直接变成错误的寿命预测有个回归测试兜底心里踏实很多。