レインフロー法 #1（準備編）

実装するための準備（勉強）です。

変動荷重の応力波形計数法

変動荷重が加わった場合、負荷の頻度処理*1をすることで損傷率計算が可能になる。
なお今回時間軸は無視しているため、変動荷重つまり応力波形の周期は損傷率に寄与しないものと仮定している。

計数法の中でも応力 - 歪線図との対応からレインフロー法 (Rainflow counting method) がよく用いられてる。*2

レインフロー法

レインフロー法とは

レインフロー法は次の条件を満たす。
＜前提＞

応力（またはひずみ、以下では応力と言う）の変化を横軸に、時間経過を下向きに取り、応力波形のピーク値を直線で結んだグラフを考える。
このグラフを多重になった屋根と見なし、屋根上の流れを停止または軒先からの落下を停止した雨滴の出発点～停止点までの応力変化量（必要に応じて平均応力や応力振幅）を1/2カウントとして計上する。

＜条件＞

雨滴は屋根の付け根から番号順に流れ始め、停止条件が満たされるまで下の屋根に流れ落ち続ける。
軒先から落下中の雨滴は次の2つの条件の一報を満たした際に落下を停止する。
- 右向きの流れの場合：右向きに流れる雨滴の出発点より左側に他の屋根の軒先が現れたとき。
- 左向きの流れの場合：左向きに流れる雨滴の出発点より右側に他の屋根の軒先が現れたとき。
屋根の一部をすでに雨滴が流れていたら、その流れは停止する。

レインフロー法の図は例えばレインフロー法 - Wikipedia 参照。

P/V 差法

とある時刻 $X_i$ での応力値 $T_i$ のデータ列 $\{(X_i, T_i)\}_{i=0}^D$ *3に対して、 $\Delta T_i \equiv T_i-T_{i-1}$ と隣接する応力変化量を定義する。ここで、以下でもデータ列となる応力値は極値であるとする。
とある点 j で小ループ除去可能とは、 $|\Delta T_{j-2}|>|\Delta T_{j-1}|\leq |\Delta T_j|$ であることを言い、データ列 $\{(X_i, T_i)\}_{i=0}^D$ から順番を変更せずに点 j-2、j-1 だけを取り除いたデータ列 $\{(X'_i, T'_i)\}_{i=0}^{D-2}$ を作ることを点 j に対し小ループ除去すると言う。なお、除去された小ループは応力波形として $T_{j-2}\rightarrow T_{j-1}$ を受けた（1サイクル）として計上する。

さて、小ループ除去しきったデータ列 $\{(x_i, t_i)\}_{i=0}^{m+k+1}$ *4 は次の性質を満たす。（証明は参考文献参照）
$0\equiv |\Delta t_0|\leq |\Delta t_1|(=t_1) \leq |\Delta t_2|\leq \cdots \leq |\Delta t_m| \leq |\Delta t_{m+1}|> |\Delta t_{m+2}|> \cdots >|\Delta t_{m+k}|$
以下ではこれを $b_0(\equiv 0)\leq b_1\leq b_2\leq \cdots \leq b_m \leq a_0 > a_1 > a_2 >\cdots >a_{k-1}$ と書く。

つまり、小ループ除去しきったデータ列は隣接する応力差（の絶対値）が徐々に（単調に）大きくなり、あるところで応力差（の絶対値）が最大をとった後に徐々に（単調に）小さくなる特徴を持った、応力値に対するジグザグ波形となる。

小ループ除去しきったデータ列に新しいデータ $(\alpha, \beta)$ を後尾に追加*5したあとで、その新たなデータ列を再び小ループ除去しきったデータ列にする処理を実施する。
このとき、データ点の極値性*6から次の2パターンが考えられる。
なお、 $a_k\equiv |\beta -t_{m+k+1}|$ と定義する。

単に追加されるのみ： $a_{k-1}>a_k$ が成り立っている場合に該当。新たなデータ列は既に小ループ除去しきったデータ列となっているので処理は不要。
小ループ除去を繰り返す： $a_{k-1}>a_k$ が成り立たない場合に該当。 $a_{k-2}>a_{k-1}\leq a_k$ であるので点 $(\alpha, \beta)$ にて小ループ除去可能。以下、「単に追加されるのみ」となるまでこれを繰り返す。

ここで注意がある。小ループ除去を繰り返した結果、データ列が $\{(x_0,t_0),\cdots, (x_m,t_m), (x_{m+1},t_{m+1}),(\alpha, \beta)\}$ まで達した場合*7である。
このとき、既に小ループ除去しきったデータ列ではあるが、 $b_0(\equiv 0)\leq b_1\leq b_2\leq \cdots \leq b_m \leq a_0(=|t_{m+1}-t_m|) \leq |\beta -t_{m+1}|$ が成り立つので新たに、 $b_{m+1}:=a_0, a_0':= |\beta -t_{m+1}|$ としておける。これは応力差（の絶対値）が最大となる箇所がずれたことを示している。

つまり、小ループ除去工程において $a_0\sim a_{k-1}$ の部分をどれだけ小ループ除去するかだけを考えれば十分である。この部分を記憶系列と呼ぶ。

上記の点をまとめると、

実波形をピークカウント法にて極値のみのデータ列に変換する*8。ここで必要に応じて間引きを行う。
データ0点目として原点を定義し、 $b_0(=0)$ を定める。
データ1点目を読み込み、 $a^{(0)}_0:=t_1-b_0$ を定める。
データ2点目を読み込み、 $a^{(0)}_1:=t_2-a^{(0)}_0$ を定め、記憶系列を小ループ除去し、必要に応じ $b_1$ を定義する。小ループ除去した応力波形を記録する。
データ3点目を読み込み、 $a^{(1)}_k:=t_3-a^{(1)}_{k-1}$ (k : 先の処理による）を定め（以下略）

として逐次的に処理することで

小ループ除去され記録された1サイクル分の応力波形（複数）
小ループ除去しきったジグザグな応力波形（1つ、各”ジグ”・各”ザグ”が1/2サイクル）

が得られる。これを P/V 差法（P/V Difference Method）と呼ぶ。

レインフロー法と P/V 差法

P/V 差法を実施した後に小ループ除去しきった応力波形を各ピーク間ごとに読み取り、小ループ除去され記録された応力波形の集団に加える。
これによりレインフロー法が実行される。

数式で表すとちょっとわかりにくいが P/V 差法で行っている小ループ除去は条件2、3で途中で遮られる（小さな）流れそのものである。
また、小ループ除去しきった応力波形は左から右、もしくは右から左に大きく流れる雨滴の流れの総計である。

ゆえに P/V 差法によって

逐次的な手法でレインフロー法を実施できる
全データに対して更新することなく、記憶系列のみ処理で良い
新規にデータ入力があった際も必要な範囲のみの処理にて算出可能

が可能となった。

参考文献

簡明にされたレインフローアルゴリズム「Ｐ／Ｖ差法」について | CiNii Research：レインフロー法の実装アルゴリズムの1つである P/V 差法について書かれた論文。本記事はこれを簡潔にまとめた。
「Rain Flow Method」の提案とその応用 | CiNii Research：レインフロー法を提唱した論文。他の計数法では何が不足でレインフロー法だとどういうふうに取り入れられるのかが解説されている。
疲労設計便覧：「これは持ってて損はない」と言われて速攻買った。辞書的にも勿論使えるが各話題に対しても下手な本より詳しい。正直諭吉1枚でこのレベルが手に入るなら安い。
製品開発のための疲労破壊事故の解析と強度対策：このシリーズは結構お気に入り。実際の設計でどうか、過去どうだったのかといったところに力点が置かれている。
疲労解析のためのレインフローカウント - MATLAB rainflow - MathWorks 日本：MathWorksのヘルプ。割りとよくまとまってたので備忘録として。中身は見てない（笑）