軌道座屈確率シミュレーションの概要[1,2]

はじめに

 軌道系には様々な不確実性が存在する.その中でも,軌道座屈に大きく影響を及ぼすものとして,①初期通り変位,②道床横抵抗力,③中立温度の3つが挙げられる[3].特に①と②は,軌道の座屈強度に直接関わるものである.軌道座屈強度が初期通り変位振幅に対して鋭敏性を有していることはⅠ.で既に述べた.なお,初期通り変位波形は従来の軌道座屈解析で用いられている様な所定の幾何形状を初めから有しているものでは無く,座屈後に起こる局所化現象の結果としてその様な波形が成長すると考えるのが合理的である.この様に考えると,座屈波形が成長する前の形状を初期通り変位として設定して,それを初期条件に座屈解析を行うのが適切であると言える.常時の軌道通り変位波形は,測定データを見る限り,概ねランダム性を有するものと考えられる.そこで本研究では,所定のパワースペクトル密度(PSD)に従うランダム波形として初期通り変位波形を生成し,これらに対して軌道座屈解析を行う.その下で,モンテカルロシミュレーション(MCS)を実施し,座屈強度の確率特性を評価する.以下では,その一連の過程について述べる.

軌道のモデル化

 左右レールとまくらぎから構成される軌きょうを図-1の様にモデル化する.レールは,横方向たわみと長手方向の伸縮を考慮したEulerばりにより表現する.まくらぎは等間隔Lで配置し,横方向と軌道長手方向変位を考慮した剛な棒として与える.レール締結部には,回転バネkRと横方向バネkTを設定する.また,まくらぎに作用する道床横抵抗力fTと縦抵抗力fLについては,次の非線形式で与える[4].

式(1) 図-1 軌道のモデル化
図-1 軌道のモデル化

ここで,uST,uSLはまくらぎの横変位と軌道長手方向変位,f0T,f0Lは最終道床横・縦抵抗力,aT,aLは作用力が最終道床抵抗力の1/2を与える時の変位である.ただし,道床横抵抗力はまくらぎ1本分に対応するのに対し,道床縦抵抗力は図-1に示す様にレール1本分(まくらぎ1/2本分)に対応している.
 なお,式(1)の道床横・縦抵抗力とまくらぎ変位との関係は,図-2 の様な骨格曲線で与えられ,最終道床横・縦抵抗力f0T,f0Lに漸近する.また,除荷時は図中に示した様に,作用力がゼロになるまで初期剛性(f0T/aT,f0L/aL)で戻り,その後骨格曲線に沿って逆方向に抵抗力が作用するように設定した.

図-2 道床横・縦抵抗力の載・除荷曲線
図-2 道床横・縦抵抗力の載・除荷曲線

初期通り変位原波形の生成

 軌道検測装置により過去に取得された通り変位の10m弦正矢データ[5]から距離相関を直接評価すると,原波形の距離に関する自己相関関数(距離相関)R(x)は概ね次式で近似することができる.

式(2)

ここで,σdは通り変位波形の標準偏差と相関長(相関性が存在する距離の代表値)である.なお,ある区間において得られた既存の検測データに関しては,標準偏差と相関長がそれぞれ1cm以下,および1.7m程度と推定された.
 また,同一データでも,10m弦正矢測定データのPSDから原波形のPSDを求め,その逆Fourier変換より概ね2(1/m)以下の波数域特性を反映した初期通り変位原波形の距離相関を求めると,それは次式で近似できる.

式(3)

なお,式(3)による場合の相関長はd=2.84mと推定された.
 長さlの軌道区間をM個のはり要素によりx軸方向に等分割する場合を考える.その際に,i番節点のx座標xiを次式で与えるものとする.

式(4)

レール初期通り変位波形のxiにおける値をw0iとし,その離散データを成分に持つベクトルを{W0}とおく.すると,当該ベクトルに関する分散・共分散行列[C]は次式で与えられる.

式(5)

ここで,E( )は期待値を,また( )Tは転置を意味する.式(2)より,行列[C]の成分cijは次式で与えられる.

式(6)

[C]に関して次の固有値問題を設定する.

式(7)

 初期通り変位波形ベクトル{W0}を,期待値がゼロであり,且つ式(5)の分散・共分散行列で与えられる正規確率過程に従うものとすると,{W0}は次式により生成することができる.

式(8)

ここで,[Λ1/2]は固有値の平方根√(λ_i )を対角項に持つ対角行列,[Φ]はそれに対応する固有ベクトル{ϕi}を縦ベクトル成分に持つ行列である.また{ξ}は,期待値ゼロ,分散1の標準正規乱数を成分に持つベクトルである.
 以上の手順により,左右レールそれぞれに対してランダムな通り変位波形を設定する.

モンテカルロ法における軌道座屈解析

 本研究では,式(2), (3)に従うランダムな通り変位波形を式(8)に基づき多数生成し,各々の波形の下で軌道座屈解析を行い飛び移り座屈温度を求める,モンテカルロシミュレーションを実施する.その際の軌道座屈解析では,有限変位理論に基づきレールをはり要素(有限要素)で離散化する.なお,通り変位をランダムな波形で設定する場合,座屈発生箇所が事前に特定できないため,所定のはり要素節点のたわみを制御変数とした変位増分解析は適用できない.そこで,節点変位ベクトルとレール温度とで与えられる一般座標空間において弧長増分法[6]を適用し,つり合い経路上の飛び移り座屈点をレール温度の極大点として探索する手法を構築した.ただし,1ステップ目のみ温度増分を規定し,それ以降におけるつり合い経路の探索方向を設定した.

確率密度関数の作成

 得られた座屈温度から,次のカーネル密度推定[7]により確率密度関数p(t)を求める.

式(9)

ここで,tはレール温度,σtは座屈温度データtiから求めた標準偏差,NはMCSのサンプル数である.
 式(9)に対応する確率分布関数P(t)は次式により評価する.

式(10)

ここで,Erf(t)は誤差関数であり,本研究では次式により定義する.

式(11)

異なる軌道長における座屈確率の評価

 本研究では一定の軌道長の下で座屈温度の確率分布を求める.なお,軌道長が異なれば,当然全区間にわたる座屈発生確率は異なるものとなる.ただし,ある軌道長の下で座屈確率が一旦得られたならば,異なる軌道長におけるそれは理論的に導出することができる.以下では,レール軸力が概ね一定値をとる不動区間において軌道座屈が発生するものと仮定して議論を進める.
 長さlRの不動区間において,温度tまでに飛び移り座屈を生ずる確率PlR(t)は次式で与えられる.

式(12)

ここで,PlRは不動区間長がlRの時に得られた座屈温度の確率密度関数である.
 すると,当該軌道が温度tまで座屈しない確率P ̅_(l_R ) (t)は次式で与えられる.

式(13)

 同様の力学状態にある長さn×lRの軌道が温度tまで座屈しない確率は,n個ある長さlRの区間全てが座屈しない確率であるので,それは次式で与えられる.

式(14)

ただし,各区間における座屈発生は互いに独立であるものと仮定する.
 したがって,不動区間長n×lRの軌道が温度t以下で座屈する確率は次式で与えられる.

式(15)

また確率密度関数は,式(15)を微分することで次式によって与えることができる.

式(16)

参考文献

  1. 阿部和久,水野雄太,紅露一寛 : 通り変位波形におけるバラツキが軌道座屈強度の確率特性に及ぼす影響,鉄道工学シンポジウム論文集,No.24, 167-174, 2020.
  2. 岩井 翔,阿部和久,紅露一寛 : 通り変位と道床横抵抗力のバラツキを考慮した軌道座屈余裕度の確率的評価,鉄道工学シンポジウム論文集,第25号,69-76, 2021.
  3. U.S. Department of Transportation : Track buckling prevention: theory, safety, concepts, and applications, Chap.4, National Technical Information Service, 2013.
  4. 宮井 徹 : エネルギー法による軌道座屈の数値解析,鉄道技術研究報告, No.1271, 1984.
  5. 千葉颯兵,阿部和久,小松佳弘,紅露一寛 : 通り変位測定データに基づくレール軸力推定法に関する理論的検討,J-RAIL2017, CD-ROM, S2-14-4, 2017.
  6. 岩崎英治,松野純一,長井正嗣 : 弧長法のための一反復解法と弧長自動設定法,応用力学論文集,Vol.5, pp.207-216, 2002.
  7. Bowman, A.W. and Azzalini, A. : Applied smoothing techniques for data analysis, The kernel approach with S-plus illustrations, p.31, Oxford Science Publications, 1997.