JP6672122B2 - 自励振動評価方法 - Google Patents

Description

本開示は、構造物の振動解析の分野に関し、より具体的には、流体中に配置された管群に生じる自励振動を評価するための自励振動評価方法に関する。

流体中に配置される管群で生じる振動現象について、コンピュータのような電子演算装置を利用して解析を行う、いわゆる流動解析技術が知られている。このような流動解析技術は、例えば、ボイラや原子力機器などの蒸気発生器を構成する伝熱管群のように、熱交換媒体である流体の励振力を受けて振動する管群の挙動を解析するために応用されている。

例えば加圧水型原子炉（ＰＷＲ）に使用される蒸気発生器では、原子炉から供給される一次冷却水を流す伝熱管を並列配置して管群とし、管群の伝熱面の外表面に二次冷却水を流して熱交換している。このような蒸気発生器では、二次冷却水の流れの流速を高めることにより熱交換効率を促進することが試みられるが、流速がある限界流速を超えてしまうと、自励振動（流力弾性振動）が生じることがある。自励振動は、管群の動きと流体流れが相互に影響しあう不安定振動であり、管群の損傷を招くおそれがある。

このような自励振動の発生を評価するための技術として、例えば特許文献１には、多数のパラメータを用いながらも時間とコストを削減可能なＣＦＤ（数値流体力学）による数値シミュレーションを用いて限界流速を予測することが開示されている。

特開２０１５−０２６２５９号公報

近年、Ｕ字管のようなＵベンド部を有する管群において、流体の流れ方向に沿って、流力弾性振動のような自励振動現象が発生する可能性が指摘されている。Ｕ字管は、各管の隙間に配置された振れ止め部材（振動抑制部材）で支持されており、このような流れ方向に沿った振動現象は、管と振れ止め部材との間の摩擦力によって抑制されている。特許文献１のような従来手法では、線形の減衰（構造減衰・二相減衰など）については考慮されているものの、各Ｕ字管と振止め部材とが接触する全支持点で押付力がゼロであるとの仮定を前提としており、流れ方向に沿った振動現象を抑制するための摩擦力については考慮がなされていない。そのため、流れ方向に沿った自励振動を適切に評価することができていない。

本発明の少なくとも１実施形態は上記事情に鑑みなされたものであり、管群と支持部材との間の摩擦減衰を考慮することにより、自励振動を適切に評価可能な自励振動評価方法を提供することを目的とする。

（１）本発明の少なくとも１実施形態に係る自励振動評価方法は上記課題を解決するために、流体中に支持部材に支持されて配置される管群の自励振動を評価するための自励振動評価方法であって、前記管群の一つ以上の固有モードの各々について、
前記流体の励振力に対応する負減衰比を変化させながら、前記管群の振動振幅の変化を模擬する時刻歴応答解析を実行する時刻歴応答解析ステップと、
前記時刻歴応答解析において前記管群の振動振幅の変化が発散する最小の負減衰比に基づいて、前記流体の限界流速を算出する限界流速算出ステップと、
前記流体の想定流速を入力する入力ステップと、
前記流体の想定流速を前記限界流速と比較することにより、前記固有モードごとに前記管群の自励振動を評価する評価ステップと、を備える。

上記（１）の方法では、流体の励振力に対応する負減衰比を変化させながら、管群の振動振幅の変化を模擬する時刻歴応答解析を実行し、振動振幅の変化が発散する最小の負減衰比に基づいて、流体の限界流速を算出している。ここで、管群の振動振幅の変化が発散する最小の負減衰比は、管群を表現する振動系が自励振動を起こさずに耐え得る最大の負減衰比であり、自励振動を抑制するために最大限作用させることが可能な摩擦減衰比に相当する。その結果、上記（１）の方法によれば、流体中に配列された複数の管から成る管群が流体の励振力に抗して支持部材との間の摩擦力によって支持されている場合に、管群に作用する摩擦減衰効果を考慮して自励振動を適切に評価することができる。

（２）幾つかの実施形態では、上記（１）の方法において、前記時刻歴応答解析は、前記管群の振動解析モデルに対して前記負減衰比に対応する励振力を外力項として作用させた際に生じる振動振幅を時系列的に模擬する演算を含み、
前記振動解析モデルは、前記管群と前記支持部材との間に作用する接触荷重の分布を仮定することによって、前記管群と前記支持部材との間の摩擦力の大きさを規定する。

上記（２）の方法では、管群と支持部材との間の摩擦力の大きさを規定する振動解析モデルを構築した上で、当該振動解析モデルに対して上記負減衰比に対応する励振力を作用させた際に生じる振動振幅を時系列的に模擬している。従って、上記（２）の方法によれば、上記負減衰比に対応する励振力が管群と支持部材との間の摩擦力によって減殺される効果を考慮して、管群の振動振幅の変化が発散する最小の負減衰比を求めることが可能となる。

（３）幾つかの実施形態では、上記（１）または（２）の方法において、前記時刻歴応答解析は、
前記管群と前記支持部材との間の摩擦力に応じて消散される前記自励振動のエネルギー消散量に対応する第１減衰比と前記負減衰比との間の相殺関係に基づいて前記管群の実効減衰比を算出し、
前記算出された実効減衰比に基づいて、前記管群の振動振幅を時系列的に推定する演算を含む。

上記（３）の方法では、管群と支持部材との間の摩擦力に応じて消散される自励振動のエネルギー消散量に対応する第１減衰比が上記負減衰比との間で相殺関係にある点に着目して管群全体の実効減衰比を算出している。そして、上記（３）の方法では、上記実行減衰比に基づいて管群の振動振幅を時系列的に推定している。従って、上記（３）の方法によれば、自励振動のエネルギーが管群と支持部材との間の摩擦力に応じて消散される効果を、上記エネルギーの消散量に対応する第１減衰比と上記負減衰比との間の相殺効果として評価することができる。その上で、上記（３）の方法によれば、このような相殺効果を考慮して、管群の振動振幅の変化が発散する最小の負減衰比を求めることが可能となる。

（４）幾つかの実施形態では、上記（３）の方法において、前記時刻歴応答解析は、
前記管群の振動振幅の変化に伴って、前記負減衰比が前記第１減衰比と一致した時点において前記振動振幅が発散すると判定する。

上記（３）の方法では、自励振動のエネルギー消散量に対応する第１減衰比と流体の励振力に対応する負減衰比との間の相殺効果として管群の振動特性を評価し、当該相殺効果を考慮して、管群の振動振幅の変化が発散する最小の負減衰比を求めている。その上で、上記（４）の方法では、管群の振動振幅の変化に伴って、上記負減衰比が上記第１減衰比と一致した時点において管群の振動振幅が発散すると判定する。その結果、上記（４）の方法によれば、限界流速時における流体の励振力に対応する負減衰比を、自励振動のエネルギー消散量に対応する第１減衰比と均衡する負減衰比として推定することが可能となる。

（５）本発明の少なくとも１実施形態に係る自励振動評価方法は上記課題を解決するために、流体中に支持部材に支持されて配置される管群の自励振動を評価するための自励振動評価方法であって、
前記流体の想定流速を取得する想定流速取得ステップと、
前記想定流速を限界流速と仮定し、前記限界流速と前記管群全体の負減衰比との間の相関に基づいて、前記想定流速に対応する負減衰比を算出する負減衰比算出ステップと、
前記管群の自励振動を模擬する演算に前記負減衰比を入力して実行した際に、前記管群の振動振幅が発散するか否かに基づいて前記管群の自励振動を評価する評価ステップと、
を備える。

上記（５）の方法では、取得された想定流速を仮の限界流速と仮定し、仮定された当該仮の限界流速に対応する負減衰比を入力して管群の自励振動を模擬する演算を実行した際に、管群の振動振幅が発散するか否かに基づいて自励振動が生じるか否かを判定する。言い換えると、上記（５）の方法では、仮の限界流速に基づいて管群の自励振動を模擬する演算を実行した際に、管群の振動振幅が発散するか否かに応じて、当該仮の限界流速が真の限界流速を上回るか否かを検査している。従って、上記（５）の方法によれば、仮の限界流速として仮定された流速で流体を流した場合に、管群の自励振動が実際に起きるか否かを、管群の自励振動を模擬する模擬演算によって高精度に予測することができる。

（６）幾つかの実施形態では、上記（５）の方法において、前記想定流速取得ステップは、前記管群に含まれる各管に作用する前記流体の動圧、前記各管の密度または前記各管の振幅の少なくとも１つについて前記各管の長さ方向に沿った分布に基づいて前記流体の有効流速を算出する有効流速算出ステップを含み、
前記負減衰比算出ステップでは、前記有効流速を前記想定流速として前記負減衰比を算出する、ことを特徴とする。

上記（６）の方法では、管群に含まれる各管に作用する流体の動圧、各管の密度または各管の振動振幅が各管の長さ方向に沿って異なる場合、上述した動圧、密度または振動振幅の長さ方向に沿った分布に基づいて流体の有効流速を算出する。その上で、上記（６）の方法では、当該有効流速を仮の限界流速と仮定して負減衰比を算出する。従って、上記（６）の方法によれば、管群に含まれる各管に作用する流体の動圧、各管の密度または各管の振動振幅が各管の長さ方向に沿って異なる場合でも、管内の場所による流速の違いを考慮して負減衰比を算出するための単一の流速値を求めることができる。

（７）幾つかの実施形態では、上記（１）〜（６）の方法において、前記管群は、同一の平面内において延在するとともに、互いに曲率中心を共有し、且つ、互いに曲率半径が異なる曲り部を有する複数のＵ字管によって形成される管列を少なくとも一つ含み、
前記支持部材は、前記管列を挟むように前記平面に沿って延在する形で前記管列の両側に配置された
少なくとも一対の振止め部材を含み、
前記管群を通る流体の励振力に抗して、振止め部材との摩擦力により支持される前記管群の前記平面に沿った方向における流力弾性振動の安定判別を行う。

一般的な熱交換器において、管群は、Ｕ字形状の曲り部をそれぞれ有する複数のＵ字管により構成され、振止め部材が、当該曲り部を含む平面に直交する面外方向において、隣接するＵ字管の曲り部の間に挿入される構造を有していても良い。その場合、面外方向において、隣接する管列の間に挿入される振止め部材が、面外方向への各Ｕ字管（曲り部）の動きを規制しているので、面外方向に作用する加振力に対しては、管群全体が一体となって振動する。しかしながら、当該曲り部を含む平面に沿った面内方向に配列された一連のＵ字管は、両側の振止め部材との間の摩擦力によってのみ規制されている。従って、上記（１）〜（６）の方法において、各管が振動する方向は面内方向とほぼ一致し、各管が隣接する振止め部材と衝突することにより受ける接触荷重も主として面内方向の摩擦力である。

そこで、上記（７）の方法では、曲り部を有するＵ字管から成る管列が延在する平面内（面内方向）において、管列に隣接する振止め部材から受ける摩擦力が各Ｕ字管に作用する励振力に抗して作用するものとして、管群の面内方向における流力弾性振動の安定判別を適切に行うことができる。

（８）幾つかの実施形態では、上記（１）〜（７）の方法において、前記管群は、ＰＷＲ原子力発電設備の蒸気発生器の伝熱管群である。

これにより、上記（８）の方法では、蒸気発生器などの熱交換器が加圧水型原子炉を含む原子力関連施設に設けられるものである場合、熱交換用の流体中に配置された当該管群がどの程度の流速まで自励振動を起こさずに耐えられるかを事前に評価することができる。その結果、当該伝熱管群の構造を、防振性を考慮して設計することが可能となる。

本発明の少なくとも１実施形態によれば、管群と支持部材との間の摩擦減衰を考慮することにより、自励振動を適切に評価可能な自励振動評価方法を提供できる。

一実施形態に係る伝熱管群のＵベンド部の斜視図である。 振止め部材による支持構造の一例を面内方向から見た図である。 振止め部材による支持構造の一例を面外方向から見た図である。 一実施形態に係る自励振動評価方法を実行するためのコンピュータ装置を示す図である。 図４Ａに示すコンピュータ装置の演算部の内部構成を示す図である。 時刻歴応答解析によって求めた負減衰比と振動振幅の関係を示す曲線グラフである。 一実施形態に係る自励振動評価方法の実行手順を示すフローチャートである。 さらに別の一実施形態に係る自励振動評価方法を実行するためのコンピュータ装置を示す図である。 さらに別の一実施形態に係る自励振動評価方法の実行手順を示すフローチャートである。 自励振動の発生を判定するための安定判別マップを示す図である。 安定限界における換算流速と減衰比との間の相関関係を示す図である。

以下、添付図面を参照して本発明の幾つかの実施形態に係る自励振動評価方法について説明する。ただし、実施形態として記載されている又は図面に示されている構成部品の寸法、材質、形状、その相対的配置等は、本発明の範囲をこれに限定する趣旨ではなく、単なる説明例にすぎない。本発明の幾つかの実施形態に係る自励振動評価方法は、流体中に配置された複数の管が流体力に抗して支持部材との摩擦力によって支持されて成る管群構造であれば、どのような管群構造であっても適用対象とすることが可能である。以下、本発明の幾つかの実施形態に係る自励振動評価方法の適用対象となり得る管群構造の一例として、図１〜図３に示す伝熱管群の構造について説明する。続いて、当該自励振動評価方法の処理内容について図４乃至図８を参照して説明する。

図１は、一実施形態に係る伝熱管群１０のＵベンド部１０ａの斜視図を示す。図２は、伝熱管群１０を、図１の面内方向Ｄ２（図１の列方向ｄ２）から見た側面図であり、図３は、伝熱管群１０を、図１の面外方向Ｄ１（図１の列方向ｄ１）から見た側面図である。なお、図１においては、図面の見やすさを考慮して一部の構成部材を省略して記載している。図１において当該省略された構成部材の一部は、図１における伝熱管群１０を側面から見た図２および図３において示されている。

幾つかの実施形態において、伝熱管群１０は、複数の伝熱管３と、複数の伝熱管３が挿通される管支持板７と、を備えており、複数の伝熱管３内を流れる流体との熱交換によって蒸気を生成するように構成される。複数の伝熱管３は、それぞれ、流体の入口側に位置する第１直管部４と、流体の出口側に位置する第２直管部５と、第１直管部４と第２直管部５との間に位置する曲り部６と、を有している。管支持板７には、第１直管部４及び第２直管部５が挿通される複数の貫通孔が形成されている。

伝熱管群１０は、Ｕ字形状の曲り部６をそれぞれ有する複数の伝熱管３により構成される。複数の伝熱管３の曲り部６によってＵベンド部１０ａが形成される。図１に示す構造においては、曲り部６の半径方向における外側（図１における上側）に向かうに従って曲り部６の曲率半径が大きな伝熱管３が互いに曲率中心を共有する形で同一平面に沿って（面内方向Ｄ２に沿って）配列される（図１の管列８）。図３は、このように伝熱管３が面内方向Ｄ２に沿って配列された管列８が複数列存在し、これら複数の管列８が、曲り部６を含む平面に直交する方向（図１の面外方向Ｄ１）において並んでいる様子を示す。

図１および図３に示すように、複数の管列８のそれぞれにおいて最外周側に位置する伝熱管３の曲り部６の曲率半径は、各管列８の面外方向Ｄ１における位置に応じて異なる。このように複数の管列８を面外方向Ｄ１に重ねながら曲り部６の曲率半径を変化させることで、伝熱管群１０の上端部において半球形状のＵベンド部１０ａが形成される。その結果、図１に示すように、曲率半径の異なる複数の曲り部６ａ_１，６ａ_２，６ａ_３，…が面内方向Ｄ２に沿って並び、曲率半径の同一な複数の曲り部６ａ_１，６ｂ_１，６ｃ_１，…が面外方向に沿って並ぶように配列される。

伝熱管群１０においては、振止め部材１２が、曲り部６を含む平面に直交する面外方向Ｄ１において、隣接する伝熱管３の曲り部６の間に挿入され、面外方向Ｄ１への複数の伝熱管３（曲り部６）の動きを規制している。例えば、図１においては、面外方向Ｄ１において並ぶ各々の管列８の両側に、面内方向Ｄ２に沿って複数の振止め部材１２が挿入され、各管列８に属する複数の伝熱管３の曲り部６の面外方向Ｄ１への動きを規制している。

図１に示すように、第１保持部材１１は、Ｕベンド部１０ａの外周、すなわち、Ｕベンド部１０ａの半球状の外周に沿って取り付けられた円弧状の棒状部材である。上述した振止め部材１２は、第１保持部材１１からＵベンド部１０ａの半球形状の径方向における内側に向かって延在している。振止め部材１２の端部１２ａには、図１に示すように第１保持部材１１が溶接されて、複数の振止め部材１２の端部１２ａを連結している。第１保持部材１１は、複数の伝熱管３が面内方向Ｄ２に沿って重ねられる管列８と直交しながらＵベンド部１０ａの半球面に沿って延在している。

図２および図３に示すように、複数の第１保持部材１１は、第２保持部材（ブリッジ）１４によって連結されていてもよい。第２保持部材１４は、Ｕベンド部１０ａの外周、すなわち、Ｕベンド部１０ａの半球状の外周に沿って配置された、円弧形状かつ板状の部材である。第２保持部材１４は、Ｕベンド部１０ａにおいて伝熱管３の曲り部６が延在する方向に沿って延在している。第２保持部材１４は、面外方向Ｄ１において並ぶように複数配置されていてもよい。

伝熱管群１０では、面外方向において、隣接する伝熱管３の曲り部６の間に挿入される振止め部材１２が、面外方向Ｄ１への複数の伝熱管３（曲り部６）の動きを規制しているので、面外方向Ｄ１に作用する加振力に対しては、伝熱管群１０全体が一体となって振動する。しかしながら、曲り部６を含む平面に沿った面内方向Ｄ２に配列された一連の伝熱管３（図１の管列８）は、両側の振止め部材１２とは接続されておらず、両側の振止め部材１２との間の摩擦力によってのみ規制されている。その結果、各伝熱管３が振動する方向は面内方向Ｄ２とほぼ一致し、各伝熱管３が隣接する振止め部材１２と衝突することにより受ける接触荷重も主として面内方向Ｄ２の摩擦力である。

例示的な一実施形態では、図１〜図３を用いて上述した伝熱管群１０は、加圧水型（ＰＷＲ）原子力発電設備において、１次冷却水と２次冷却水の間で熱交換を行うための蒸気発生器の伝熱管群として構成されてもよい。その場合、当該２次冷却水は、図１に示す面外方向Ｄ１と面内方向Ｄ２に直交する方向Ｇに沿ってＵベンド部１０ａの真上から真下に流れることにより、伝熱管３内を流れる１次冷却水との間で熱交換を行う。そのため、当該２次冷却水の流れは、Ｕベンド部１０ａの最上部において伝熱管３の曲り部６に直交する直交流となる。以上より、本発明の幾つかの実施形態に係る自励振動評価方法は、伝熱管群１０に自励振動を生じさせる限界流速を、上述した蒸気発生器において、Ｕベンド部１０ａに直交して流れる熱交換用の２次冷却水の流れの限界流速として事前に評価するために実施されてもよい。

以上のように、伝熱管群１０が加圧水型原子炉に備えられる蒸気発生器を構成する場合を例として説明するならば、原子炉から供給される一次冷却水を流す伝熱管３を並列配置して伝熱管群１０とし、この伝熱管群１０の伝熱面の外表面に二次冷却水を流して熱交換している。この蒸気発生器においては、二次冷却水の流れの流速を高めて熱交換を効率化する必要があるが、この流速がある限界流速を超えると、伝熱管群１０に自励振動が生じる場合がある。当該自励振動は、伝熱管群１０の動きと流体流れが相互に影響しあう構造的な不安定挙動であり、時間の経過に伴って振動振幅が増大し、伝熱管群１０の損傷に至る大きな問題である。

そこで、上述した蒸気発生器における伝熱管群が自励振動を起こすのを防止するため、下端部が管支持板７により支持された複数の伝熱管３は、上部のＵベンド部１０ａにて、複数の振止め部材１２が挿入されて支持される構造となっている。つまり、蒸気発生器のＵベンド部１０ａにて、同一平面に沿って配置された複数の伝熱管３から成る管列８は、その間に振止め部材１２が挿入されて支持されている。この場合、振止め部材１２と伝熱管３の曲り部６との間に作用する接触荷重は、二次冷却水の流体力により生じる自励振動のエネルギーを減衰させる制振力として働くが、与えられた伝熱管群１０の構造から当該制振力の大きさに応じて限界流速を事前に評価することは有益である。

なお、以下に述べる幾つかの実施形態では、自励振動の評価は、専ら伝熱管群１０のＵベンド部１０ａを構成する各伝熱管３の曲り部６を対象として行われる。従って、以下の幾つかの実施形態においては、伝熱管群１０のＵベンド部１０ａを指して単に伝熱管群１０と呼び、各伝熱管３の曲り部６を指して単に伝熱管６または管６と呼ぶ。

次に、本発明の幾つかの実施形態に係る自励振動評価方法と当該自励振動評価方法を実施するためのコンピュータ装置について図４〜図６を参照しながら説明する。図４Ａは、幾つかの実施形態に係る自励振動評価方法を実施するためのコンピュータ装置２０の全体構成を示す図である。コンピュータ装置２０は、演算部２１、記憶部２２、出力部２３および入力部２４を含んで構成される。例示的な一実施形態では、演算部２１は、記憶部２２に記憶されたプログラム２２ａを読み込んで実行することにより、流体ｆｌ中に振止め部材１２に支持されて配置される伝熱管群１０の自励振動を評価するための自励振動評価方法を実行する演算回路として構成されてもよい。また、この実施形態では、当該自励振動評価方法の実行に際して演算部２１が読み書きする必要のあるデータは記憶部２２内においてデータ２２ｂとして記憶されてもよい。

また、出力部２３は、演算部２１による演算結果や記憶部２２に記憶されているデータ２２ｂの一部を利用者に提示するための出力装置である。例示的な一実施形態では、出力部２３は、ディスプレイ装置などの画面表示手段を出力手段として備えるものであってもよい。また、入力部２４は、利用者が操作することによって演算部２１に各種情報やパラメータ類を表す外部データを入力するための入力装置である。例示的な一実施形態では、入力部２４は、キーボードやマウスなどを入力手段として備えるものであってもよい。

図４Ｂは、コンピュータ装置２０が備える演算部２１の内部構成を示す図である。図４Ｂを参照すると、演算部２１は、後述する限界流速Ｕｃｒを算出するための限界流速算出部２１１と、限界流速算出部２１１から限界流速Ｕｃｒの算出結果を受け取り、伝熱管群１０の自励振動の評価を行う自励振動評価部２１３とを含んで構成される。また、演算部２１は、時刻歴応答解析を行うために、限界流速算出部２１１から反復的に呼び出される時刻歴応答解析部２１２をさらに含んで構成される。一例においては、演算部２１は、汎用プロセッサにより実現されてもよい。その場合、限界流速算出部２１１、時刻歴応答解析部２１２および自励振動評価部２１３は、演算部２１が記憶部２２からプログラム２２ａを読み込むことにより演算部２１内に生成されるプログラム・モジュールとして実現されてもよい。

なお、一般的に、伝熱管群１０は、多自由度振動系としてモデル化されるので、伝熱管１０は複数の固有振動数ｆ（ｉ）（１≦ｉ≦Ｉ）を有し、それに対応して伝熱管群１０の振動は、複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の合成として表される。従って、限界流速算出部２１１による限界流速Ｕｃｒの算出は、固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について個別に実行され、固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について限界流速Ｕｃｒ（ｉ）（１≦ｉ≦Ｉ）が算出される。言い換えると、複数の固有振動数ｆ（ｉ）（１≦ｉ≦Ｉ）を有する伝熱管群１０の自励振動の評価は、複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について個別に行われる。

図４Ｂに示す限界流速算出部２１１は、以下のようにして、一つの固有モードφ（ｉ’）に対応する限界流速Ｕｃｒ（ｉ’）を算出する。まず、伝熱管群１０を通過する流体ｆｌの想定流速の値に加え、限界流速Ｕｃｒ（ｉ’）を算出するために必要な入力パラメータを入力部２４から受け取る。例えば、限界流速算出部２１１は、想定流速の値の他に、記憶部２２内のデータ２２ｂのうち、伝熱管群１０の振動解析モデルを定義するデータを識別する情報を入力部２４から受け取る。次に、振止め部材１２に支持されて配置される伝熱管群１０に対して励振力Ｆ_ｅｘを加える流体ｆｌの流速を増加させていった場合に、伝熱管群１０が固有振動数ｆ（ｉ’）に対応する自励振動を起こすようになる最小の流速を限界流速Ｕｃｒ（ｉ’）として算出する。

その際、上記限界流速Ｕｃｒ（ｉ’）の算出過程において、限界流速算出部２１１は、負減衰比ζｎ（ｉ’）の値を変化させながら、負減衰比ζｎ（ｉ’）の各値について時刻歴応答解析部２１２を反復的に呼び出す。負減衰比ζｎ（ｉ’）の各値を入力として渡されて限界流速算出部２１１から呼び出された時刻歴応答解析部２１２は、伝熱管群１０の振動振幅の変化を模擬する時刻歴応答解析を実行する。つまり、上記時刻歴応答解析は、負減衰比ζｎ（ｉ’）の値を入力として、負減衰比ζｎ（ｉ’）に対応する励振力Ｆ_ｅｘが伝熱管群１０に作用した場合の伝熱管群１０の振動振幅を算出するパラメトリック・スタディ演算である。

このようにして、限界流速算出部２１１は、負減衰比ζｎ（ｉ’）の値を変化させながら、負減衰比ζｎ（ｉ’）の各値について時刻歴応答解析の結果を時刻歴応答解析部２１２から受け取り、当該時刻歴応答解析において伝熱管群１０の振動振幅の変化が発散する最小の負減衰比である臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）を求める。最後に、限界流速算出部２１１は、上記のようにして求めた臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）に基づいて限界流速Ｕｃｒ（ｉ’）を算出し、自励振動評価部２１３に出力する。限界流速算出部２１１から固有モードφ（ｉ’）に対応する限界流速Ｕｃｒ（ｉ’）を受け取った自励振動評価部２１３は、入力部２４から入力された流体ｆｌの想定流速を限界流速Ｕｃｒ（ｉ’）と比較することにより、固有モードごとに伝熱管群１０の自励振動を評価する。つまり、この実施形態では、負減衰比ζｎを徐々に大きくしながら、負減衰比ζｎを入力とする時刻歴応答解析を反復実行することにより、疑似的に振動振幅が大きくなる過程を模擬することとなる。

なお、図４Ｂに示す限界流速算出部２１１が時刻歴応答解析部２１２を反復的に呼び出すことにより、伝熱管群１０の振動振幅の変化が発散する最小の負減衰比である臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）を求める処理を、図５を参照しながら詳しく説明すると以下のとおりである。図５の曲線グラフは、固有モードφ（ｉ’）について、負減衰比ζｎ（ｉ’）の値の変化に応じた伝熱管群１０の振動振幅の変化を表す。つまり、図５の曲線グラフは、限界流速算出部２１１が、負減衰比ζｎ（ｉ’）の値を徐々に大きくしながら、負減衰比ζｎ（ｉ’）の各値を入力として伝熱管群１０の振動振幅を求める時刻歴応答解析を反復実行し、負減衰比ζｎの各値について算出された振動振幅をプロットしたものである。

図５の曲線グラフを参照すると、負減衰比ζｎ（ｉ’）の値が１０以下であるときは、伝熱管群１０の振動振幅は負減衰比ζｎ（ｉ’）の増加に従って微増しながらもほぼ横ばいであるのに対して、負減衰比ζｎ（ｉ’）の値が１１になると、伝熱管群１０の振動振幅は急激に上昇する。つまり、負減衰比ζｎ（ｉ’）の値を徐々に大きくしながら、負減衰比ζｎ（ｉ’）の各値を入力として伝熱管群１０の振動振幅を求める時刻歴応答解析を反復実行する過程で、負減衰比ζｎ（ｉ’）の値が１１になると伝熱管群１０の振動振幅が発散するのである。このように、複数の固有モードφ（ｉ）の各々について、負減衰比ζｎ（ｉ）の値を徐々に大きくしていったときに、伝熱管群１０の振動振幅が発散する臨界点に対応する臨界負減衰比ζ_ｎ ^ｃｒ（ｉ）を求めることができる。そして、複数の固有モードφ（ｉ）の各々について、臨界負減衰比ζ_ｎ ^ｃｒ（ｉ）から限界流速Ｕｃｒ（ｉ）を算出することが可能となる。

次に、本発明の幾つかの実施形態に従い、図４Ａおよび図４Ｂに示すコンピュータ装置２０によって実行される自励振動評価方法の実行手順を図６のフローチャートに沿って説明する。図６に示すフローチャートは、伝熱管群１０が有する複数の固有振動数ｆ（ｉ）（１≦ｉ≦Ｉ）の各々に対応して、複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について伝熱管群１０の自励振動を評価する手順を示す。図６のフローチャートに示す自励振動評価方法は、伝熱管群１０の自励振動に関する以下のような基本的性質に着目して実現されている。すなわち、伝熱管群１０の振動振幅の変化が発散する最小の負減衰比ζ_ｎ ^ｃｒは、伝熱管群１０を表現する振動系が自励振動を起こさずに耐え得る最大の負減衰比であり、自励振動を抑制するために最大限作用させることが可能な摩擦減衰比ζ_ｐ ^ｍａｘに相当するという点に着目している。一の固有モードφ（ｉ）について、図６に示すフローチャートの実行が開始すると、ステップＳ５１において、限界流速算出部２１１は、負減衰比ζｎ（ｉ’）の初期値を設定し、当該負減衰比の初期値を時刻歴応答解析部２１２に渡す。

続いて、処理はステップＳ５２に進み、負減衰比ζｎ（ｉ’）の初期値を限界流速算出部２１１から受け取った時刻歴応答解析部２１２は、時刻歴応答解析部２１２が負減衰比ζｎ（ｉ’）を入力パラメータとして、伝熱管群１０の振動振幅の変化を模擬する時刻歴応答解析を実行する。例示的な一実施形態では、時刻歴応答解析部２１２が実行する上記時刻歴応答解析は、伝熱管群１０の振動解析モデルＨ（φ，ｔ）に対して負減衰比ζｎ（ｉ’）に対応する励振力Ｆ_ｅｘを外力項として作用させた際に生じる振動振幅を時系列的に模擬する演算を含んでいてもよい。また、振動解析モデルＨ（φ，ｘ）は、伝熱管群１０と振止め部材１２との間に作用する接触荷重の分布を仮定することによって、伝熱管群１０と振止め部材１２との間の摩擦力の大きさを規定するモデルであってもよい。例えば、振動解析モデルＨ（φ，ｘ）は、伝熱管群１０と振止め部材１２との間の摩擦力の大きさを、伝熱管群１０の振動を記述する運動方程式に反映させるべき付加剛性や付加減衰の成分としてモデル化するものであってもよい。

続いて、処理はステップＳ５３に進み、限界流速算出部２１１は、上述した時刻歴応答解析に負減衰比ζｎ（ｉ’）を入力した結果得られる伝熱管群１０の振動振幅の大きさを時刻歴応答解析部２１２から受け取る。続いて、限界流速算出部２１１は、伝熱管群１０の振動振幅の大きさが発散するか否かを判定する。ステップＳ５３において、伝熱管群１０の振動振幅の大きさが発散すると判定された場合には、限界流速算出部２１１は、負減衰比ζｎ（ｉ’）の現在の値を臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）の値に設定する。続いて、限界流速算出部２１１は、臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）の値に基づいて限界流速Ｕｃｒ（ｉ’）を算出して、限界流速Ｕｃｒ（ｉ’）の値を自励振動評価部２１３に渡す。続いて、自励振動評価部２１３は、限界流速算出部２１１から限界流速Ｕｃｒ（ｉ’）の値を受け取り、入力部２４から流体ｆｌの想定流速の値を受け取ると、限界流速Ｕｃｒ（ｉ’）の値を想定流速の値と比較する。最後に、自励振動評価部２１３は、当該比較の結果から、固有モードφ（ｉ’）について伝熱管群１０の自励振動の発生の有無を判定し、図６に示すフローチャートの実行を終える。

ステップＳ５３において、伝熱管群１０の振動振幅の大きさが発散しないと判定された場合には、図６のフローチャートの実行はステップＳ５１に戻り、ステップＳ５１では、負減衰比ζｎ（ｉ’）の値を微小変位量だけ増加させる。続いて、微小変位量だけ増加された負減衰比ζｎ（ｉ’）の値を入力として、ステップＳ５２およびステップＳ５３が再度実行される。つまり、微小変位量だけ増加された負減衰比ζｎ（ｉ’）の値を入力として時刻歴応答解析を再度実行し、その結果得られた伝熱管群１０の振動振幅が発散しないと判定されれば、ステップＳ５１に戻り、負減衰比ζｎ（ｉ’）の値がさらに増加させられ、以下、同様の処理が繰り返し実行される。逆に、微小変位量だけ増加された負減衰比ζｎ（ｉ’）の値について、伝熱管群１０の振動振幅が発散すると判定されれば、負減衰比ζｎ（ｉ’）の現在の値を臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）の値に設定し、臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）の値に基づいて限界流速Ｕｃｒを算出する。

ここで、固有モードφ（ｉ’）について、流体の流速が限界流速Ｕｃｒ（ｉ’）に等しい時に作用する励振力Ｆ_ｅｘに対応する臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）は、伝熱管群１０内における伝熱管６と振止め部材１２との間の摩擦減衰に対応する正の減衰比ζｐ（ｉ’）と均衡していると見なすことが可能である。そこで、例示的な一実施形態では、ステップＳ５３において、限界流速算出部２１１は、以下のようにして、臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）の値から限界流速Ｕｃｒ（ｉ’）を算出するようにしてもよい。まず、限界流速算出部２１１は、臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）が伝熱管群１０の構造によって定まる減衰比ζｐ（ｉ’）に等しいと仮定する。続いて、減衰比ζｐ（ｉ’）に対応する対数減衰率δを算出し、対数減衰率δを以下の式に代入することにより、限界流速Ｕｃｒ（ｉ’）を算出する。

ここで、上記式（１）は、流体ｆｌ中の伝熱管群１０が流体力によって自励振動（流力弾性振動）を生じる最小の流速である限界流速Ｕｃｒと伝熱管群１０の構造によって定まる対数減衰率δとの間の関係を表している。上記式（１）において、ｆは、伝熱管群１０の固有モードに対応する固有振動数であり、Ｄは、伝熱管６の直径であり、ｍは、伝熱管６の単位長さあたりの質量であり、ρは、流体の質量密度であり、Ｋは、限界係数である。要するに、この実施形態では、限界流速算出部２１１が臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）の値から限界流速Ｕｃｒ（ｉ’）を算出するための計算式を、管群の流力弾性振動を対象とするＣｏｎｎｏｒｓの安定判別式に基づいて定義している。

以上より、図４〜図６を参照しながら上述した自励振動評価方法では、固有モードφ（ｉ’）について、以下のようにして流体ｆｌの限界流速Ｕｃｒ（ｉ’）を算出している。すなわち、流体ｆｌの励振力Ｆ_ｅｘに対応する負減衰比ζｎ（ｉ’）を変化させながら、伝熱管群１０の振動振幅の変化を模擬する時刻歴応答解析を実行し、振動振幅の変化が発散する最小の負減衰比ζ_ｎ ^ｃｒ（ｉ’）に基づいて、流体の限界流速Ｕｃｒ（ｉ’）を算出している。ここで、伝熱管群１０の振動振幅の変化が発散する最小の負減衰比ζ_ｎ ^ｃｒ（ｉ’）は、伝熱管群１０を表現する振動系が自励振動を起こさずに耐え得る最大の負減衰比であり、自励振動を抑制するために最大限作用させることが可能な摩擦減衰比ζ_ｐ ^ｍａｘに相当する。その結果、この自励振動評価方法によれば、流体ｆｌ中に配列された複数の伝熱管６から成る伝熱管群１０が流体ｆｌの励振力Ｆ_ｅｘに抗して振止め部材１２との間の摩擦力によって支持されている場合に、伝熱管群１０に作用する摩擦減衰効果を考慮して自励振動を適切に評価することができる。

また、この実施形態では、流体の流速に応じて変化させられるパラメータは、流体ｆｌの励振力Ｆ_ｅｘに対応する負減衰比だけである。ここで、伝熱管群１０を多自由度振動系としてモデル化するならば、流体ｆｌの流速に応じて変化させるべき負減衰比の個数は、当該多自由度振動系の一つ以上の固有振動数ｆ（ｉ）に対応する一つ以上の固有モードφ（ｉ）の個数に等しい。しかも、振動解析においては、所定次数を上回る高次の固有モードに対応する負減衰比の自励振動への寄与は無視できるのが一般的である。従って、この実施形態によれば、伝熱管群１０を構成する伝熱管６の本数や伝熱管群１０の構造的複雑さと無関係に、伝熱管群１０の振動解析モデルのパラメータのうち、非常に僅かな一定個数のパラメータのみを流速に応じて変化させることで限界流速を求めることができる。

例示的な一実施形態では、時刻歴応答解析部２１２は、負減衰比ζｎ（ｉ’）を入力パラメータとして以下のような時刻歴応答解析を実行することにより、伝熱管群１０の振動振幅の変化を模擬してもよい。例えば、負減衰比ζｎ（ｉ’）に対応する励振力Ｆ_ｅｘは、伝熱管群１０を取り巻く流体が伝熱管６の各々に対して作用する流体力であり、以下のように算出されてもよい。すなわち、当該流体力は、伝熱管群１０を取り巻く流体について、圧力に関するポアソン方程式を解いて伝熱管群１０内の圧力場を求め、続いて、ナビエ・ストークスの方程式（Ｎ−Ｓ式）を解いて伝熱管群１０内の流速場を求めた結果を用いて算出されてもよい。

また、この実施形態では、負減衰比ζｎ（ｉ’）に対応する励振力Ｆ_ｅｘを外力項として作用させる伝熱管群１０の振動解析モデルＨ（φ，ｘ）は、伝熱管群１０の振動を模擬する多自由度振動系の運動方程式に対して、以下の等価特性を付加したモデルとして記述されてもよい。まず、伝熱管群１０の振動を記述する運動方程式に対して付加してもよい第１の等価特性は、伝熱管群１０を取り巻く流体によって伝熱管群１０の振動特性に付加される流体付加質量、流体付加剛性および流体付加減衰である。また、伝熱管群１０の振動を記述する運動方程式に対して付加してもよい第２の特性は、伝熱管６が振止め部材１２から摩擦力を受けることによって伝熱管６に加わる励振力Ｆ_ｅｘを減殺する摩擦減衰効果に対応する付加剛性および付加減衰である。

例えば、伝熱管群１０の気中振動を模擬する多自由度振動系の運動方程式に対して上述した付加質量、付加剛性および付加減衰を反映させた運動方程式は、以下のように定義することが可能である。

ここで、ベクトルｘは、伝熱管群１０の振動による変位を表す変位ベクトルであり、伝熱管群１０の自由度に対応した次元を有するベクトルである。また、Ｍ_０、Ｃ_０およびＫ_０は、伝熱管群１０に含まれる複数の伝熱管６についての気中単位長さ質量、気中単位長さ構造減衰および気中単位長さ安定剛性をそれぞれ表す質量マトリクス、減衰マトリクスおよび剛性マトリクスである。また、Ｍ、ＣおよびＫは、伝熱管群１０を取り巻く流体および伝熱管６が振止め部材１２から摩擦力を受けることによって生じる摩擦減衰効果に対応して伝熱管群１０の振動特性に付加される付加質量、付加減衰および付加剛性をそれぞれ表すマトリクスである。

以上のような振動解析モデルＨ（φ，ｘ）を使用して、固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について入力として与えられた負減衰比ζｎ（ｉ）から伝熱管群１０の振動振幅を算出するためには、例えば、以下のようにすればよい。まず、複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）に対応する複数の負減衰比ζｎ（ｉ）（１≦ｉ≦Ｉ）を要素として含む負減衰ベクトルΖｎを以下のように定義する。

続いて、振動解析モデルＨ（φ，ｘ）に対して負減衰比ベクトルΖｎに対応する励振力Ｆ_ｅｘを外力項として作用させた場合に生じると予測される強制振動を模擬実行し、伝熱管群１０の振動変位を表す変位ベクトルｘのノルムの大きさを求めればよい。

例えば、時刻歴応答解析部２１２は、上述した時刻歴応答解析を以下のように実装してもよい。まず、負減衰比ベクトルΖｎに対応する励振力Ｆ_ｅｘを、負減衰比ベクトルΖｎをパラメータとする関数Ｆ_ｅｘ（Ζｎ）で表す。なお、一例においては、関数Ｆ_ｅｘ（Ζｎ）は、固有モード座標系を変位ベクトルｘの座標系に変換する基底変換を含んでいてもよい。続いて、負減衰比ベクトルΖｎに対応する励振力Ｆ_ｅｘ（Ζｎ）を外力項とし、上記式（２）で表される伝熱管群１０の振動特性に当該外力項を作用させて生じる強制振動を想定する。これにより、負減衰比ベクトルΖｎと変位ベクトルｘとの間の相関関係を以下の式（４）のようにモデル化する。

その上で、時刻歴応答解析部２１２は、上記式（４）を複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）にわたってモード展開し、モード展開した式をモデル式として、固有モード毎に入力された負減衰比ζｎ（ｉ）から変位ベクトルｘのノルムを逆算するパラメトリック・スタディ演算処理を実行する。これにより、時刻歴応答解析部２１２は、固有モード毎に入力として与えられた負減衰比ζｎ（ｉ’）から伝熱管群１０の振動振幅を算出する時刻歴応答解析を実現することができる。

以上より、この実施形態では、伝熱管群１０内の各伝熱管６と振止め部材１２との間の摩擦力の大きさを規定する振動解析モデルＨ（φ，ｘ）を構築した上で、振動解析モデルＨ（φ，ｘ）に対して負減衰比ベクトルΖｎに対応する励振力Ｆ_ｅｘを作用させた際に生じる振動振幅を時系列的に模擬している。従って、この実施形態によれば、負減衰比ベクトルΖｎに対応する励振力Ｆ_ｅｘが伝熱管群１０内の各伝熱管６と振止め部材１２との間の摩擦力によって減殺される効果を考慮して、各固有モードφ（ｉ’）について、伝熱管群１０の振動振幅の変化が発散する最小の負減衰比ζ_ｎ ^ｃｒ（ｉ’）を求めることが可能となる。

また、別の例示的な一実施形態では、時刻歴応答解析部２１２により実行される上述した時刻歴応答解析は、以下の演算を含むようにしてもよい。まず、複数の伝熱管６と振止め部材１２との間の摩擦力に応じて消散される自励振動のエネルギー消散量Ｅ_ｒｄに対応する第１減衰比（正の減衰比）ζｐと負減衰比ζｎとの間の相殺関係に基づいて伝熱管群１０の実効減衰比ζ_ｅｆｆを算出する。続いて、実効減衰比ζ_ｅｆｆに基づいて、伝熱管群１０の振動振幅を推定する。

この実施形態では、管群と支持部材との間の摩擦力に応じて消散される自励振動のエネルギー消散量に対応する第１減衰比ζｐが負減衰比ζｎとの間で相殺関係にある点に着目して伝熱管群１０全体の実効減衰比ζ_ｅｆｆを算出している。そして、この実施形態では、実行減衰比ζ_ｅｆｆに基づいて伝熱管群１０の振動振幅を時系列的に推定している。従って、この実施形態によれば、自励振動のエネルギーが伝熱管群１０内の複数の伝熱管６と振止め部材１２との間の摩擦力に応じて消散される効果を、上記エネルギーの消散量に対応する第１減衰比ζｐと負減衰比ζｎとの間の相殺効果として評価することができる。その上で、この実施形態によれば、このような相殺効果を考慮して、伝熱管群１０の振動振幅の変化が発散する最小の負減衰比である臨界負減衰比ζ_ｎ ^ｃｒを求めることが可能となる。

また、別の例示的な一実施形態では、時刻歴応答解析部２１２により実行される上述した時刻歴応答解析は、以下のように実現されてもよい。すなわち、上述した時刻歴応答解析は、伝熱管群１０の振動振幅の変化に伴って、負減衰比ζｎが第１減衰比ζｐと一致した時点において振動振幅が発散すると判定するようにしてもよい。

なお、この実施形態では、負減衰比ζｎは、伝熱管群１０の振動振幅に応じて非線形に増加するが、第１減衰比ζｐは、振動振幅に応じて非線形に低下する性質を有するものとしてモデル化されている。つまり、振動系における減衰比ζとは、励振エネルギーの消散量を振動の振幅に応じたエネルギーで除算して得られる比であるから、伝熱管６の振動振幅が大きくなるほど振止め部材１２との摩擦力に対応する第１減衰比ζｐは小さくなる。逆に、振動振幅の増大とともに互いに隣接する伝熱管６同士の間で相手の振動が相互に励振力Ｆ_ｅｘとして作用するので、負減衰比ζｎは振動振幅の増大とともに非線形に大きくなる。また、限界流速算出部２１１が実行する反復処理は、負減衰比ζｎを徐々に大きくしながら、負減衰比ζｎを入力とする時刻歴応答解析を反復実行することにより、疑似的に振動振幅が大きくなる過程を模擬するものである。そこで、この実施形態では、伝熱管群１０の振動振幅の増加に伴って非線形に増加する負減衰比ζｎが振動振幅の増加に伴って非線形に低下する第１減衰比ζｐと一致した時点において振動振幅が発散すると判定するようにしている。

以上より、この実施形態においても、自励振動のエネルギー消散量に対応する第１減衰比ζｐと流体ｆｌの励振力Ｆ_ｅｘに対応する負減衰比ζｎとの間の相殺効果に基づいて振動特性を評価し、当該相殺効果を考慮して、伝熱管群１０の振動振幅の変化が発散する最小の負減衰比ζ_ｎ ^ｃｒを求めている。その上で、この実施形態では、伝熱管群１０の振動振幅の変化に伴って、負減衰比ζｎが第１減衰比ζｐと一致した時点において伝熱管群１０の振動振幅が発散すると判定する。その結果、この実施形態によれば、限界流速時における流体の励振力Ｆ_ｅｘに対応する負減衰比ζ_ｎ ^ｃｒを、限界流速時における自励振動のエネルギー消散量に対応する第１減衰比ζ_ｐ ^ｍａｘと均衡する負減衰比として推定することが可能となる。

図４〜図６を用いて上述した実施形態に係る自励振動評価方法は、流体ｆｌ中に配置された伝熱管群１０が自励振動を起こす限界流速を予測可能であった。他方、本発明のさらに別の実施形態では、伝熱管群１０の自励振動評価方法は、流体ｆｌの流速を想定流速として入力した場合に、伝熱管群１０が自励振動を起こすか否かを判定するようにしてもよい。図７は、この実施形態に係る自励振動評価方法を実行するコンピュータ装置２０’が備える演算部２１’の内部構成を示す図である。なお、コンピュータ装置２０’の構成は、演算部２１が演算部２１’で置き換えられている点を除いて図４Ａに示すコンピュータ装置２０と同一である。

図７を参照すると、演算部２１’は、有効流速算出部２１５、負減衰比算出部２１６、時刻歴応答解析部２１７および自励振動評価部２１８を含んで構成される。また、演算部２１’が備える有効流速算出部２１５、負減衰比算出部２１６、時刻歴応答解析部２１７および自励振動評価部２１８は、図８に示すフローチャートに従って、この実施形態に係る自励振動評価方法を実行する。一例においては、演算部２１’は、汎用プロセッサにより実現されてもよい。その場合、有効流速算出部２１５、負減衰比算出部２１６、時刻歴応答解析部２１７および自励振動評価部２１８は、演算部２１’が記憶部２２からプログラム２２ａを読み込むことにより演算部２１内に生成されるプログラム・モジュールとして実現されてもよい。以下、固有モードφ（ｉ’）に対応する自励振動の評価を行う場合を仮定し、図８に示すフローチャートに沿って、有効流速算出部２１５、負減衰比算出部２１６、時刻歴応答解析部２１７および自励振動評価部２１８の動作を説明する。

まず、図８のステップＳ８１において、有効流速算出部２１５は、入力部２４から後述する流速算出用パラメータを受け取り、想定流速を取得する処理を開始する。続いて、ステップＳ８２において、有効流速算出部２１５は、伝熱管６の内外を流れる流体ｆｌの有効流速Ｕｅを想定流速として算出する。例示的な一実施形態では、伝熱管群１０に含まれる各伝熱管６に作用する流体ｆｌの動圧、各伝熱管６の密度または各伝熱管６の振動振幅の少なくとも１つについて各伝熱管６の長さ方向ｙに沿った分布に基づいて流体ｆｌの有効流速Ｕｅを算出するようしてもよい。

この実施形態では、上記の流速算出用パラメータは、伝熱管６の長さ方向ｙに沿って分布する流体密度ρ（ｙ）、流体の流速Ｕ（ｙ）、伝熱管６の単位長さあたりの質量密度ｍ（ｙ）、モード形状ψ（ｙ）の値の分布データを含んでいてもよい。ここで、流体密度ρ（ｙ）は、伝熱管６内部を流れる流体の長さ方向ｙに沿った密度分布および伝熱管６の外側の流体が伝熱管６によって押し除けられることによる押し除け質量の密度分布の両者を加味して推定される。また、モード形状ψ（ｙ）とは、伝熱管６の長さ方向位置ｙにおいて、伝熱管６が振動によって基準形状から変位した場合の変位量を、複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）間での相対振幅比率を考慮して定量化した値である。つまり、伝熱管６の振動振幅の評価に影響を与える流速Ｕ、流体密度ρおよび伝熱管密度ｍなどは、伝熱管６の長さ方向位置ｙによって異なる。そこで、伝熱管６の長さ方向ｙに沿ったこれらの値の分布を考慮して、伝熱管６に作用する励振力Ｆ_ｅｘに対応する実効的な流速を有効流速Ｕｅとして求めている。

例えば、例示的な一実施形態では、有効流速算出部２１５は、上述した流速算出用パラメータを使用して、以下の式に基づいて有効流速Ｕｅを算出するようにしてもよい。

上記式（５）において、“ρ_０”および“ｍ_０”は所定の定数であり、“ｌ”は、伝熱管６の長さである。また、上記式（５）の分子は、以下のように算出される。まず、伝熱管６の長さ方向位置ｙにおける流体の流速Ｕ（ｙ）の２乗値に流体密度ρ（ｙ）を乗算して、長さ方向位置ｙにおいて伝熱管６に作用する流体ｆｌの動圧を求める。続いて、長さ方向位置ｙにおける動圧にモード形状の２乗値を乗算し、当該乗算の結果を長さ方向ｙに沿って伝熱管６の長さｌにわたって線積分している。つまり、上記式（５）の分子は、伝熱管６に作用する流体ｆｌの動圧を長さ方向位置ｙにおけるモード形状ψ（ｙ）の２乗値で重みづけ平均した値に等しい。また、上記式（５）の分母は、伝熱管６の長さ方向位置ｙにおける質量密度をモード形状の２乗値で重み付けして長さ方向ｙに沿って線積分することにより平均化したものである。

続いて、図８のフローチャートの実行はステップＳ８３に進み、負減衰比算出部２１６は、有効流速算出部２１５から有効流速Ｕｅを受け取り、有効流速Ｕｅを想定流速として負減衰比ζｎ（ｉ’）を算出する処理を実行する。具体的には、負減衰比算出部２１６は、上述した想定流速を仮の限界流速Ｕｃｒ（ｉ’）と仮定し、限界流速Ｕｃｒ（ｉ’）と伝熱管群１０全体の負減衰比ζｎ（ｉ’）との間の相関に基づいて、想定流速に対応する負減衰比ζｎ（ｉ’）の具体的な値を算出する。

例示的な一実施形態では、負減衰比算出部２１６は、以下の式に基づいて、仮の限界流速として仮定された有効流速Ｕｅから負減衰比ζｎ（ｉ’）を算出するようにしてもよい。

上記式（６）は、上記式（１）に基づいて導出された関係式である。具体的には、上記式（１）を、対数減衰率δを左辺とする式に書き換え、対数減衰率δに２πで除算した値が減衰比ζｐ（ｉ’）に等しい点に着目し、対数減衰率δを２π×ζｐ（ｉ’）で置き換える。その結果、流体ｆｌ中の伝熱管群１０が流体力によって自励振動（流力弾性振動）を生じる最小の流速である限界流速Ｕｃｒ（ｉ’）と伝熱管群１０の摩擦減衰構造によって定まる摩擦減衰比ζｐとの間の関係を規定する関係式が得られる。その上で、限界流速Ｕｃｒ（ｉ’）と摩擦減衰比ζｐ（ｉ’）との間の関係式に無次元化された流速（Ｕｅ／Ｕｃｒ（ｉ’））^２を乗じることで上記式（６）が導かれる。

つまり、流体ｆｌの流速を増加させていった場合に、流体ｆｌ中の伝熱管群１０が自励振動を開始する瞬間の状態を安定限界と呼ぶならば、上記式（１）は、安定限界時における流速と、伝熱管群１０の摩擦減衰構造によって定まる摩擦減衰比ζｐ（ｉ’）との間の関係を規定している。言い換えると、上記式（１）は、安定限界時における流速を、安定限界時における伝熱管群１０の摩擦減衰比ζ_ｐ ^ｍａｘ（ｉ’）に換算する換算式でもある。さらに、上述したように、流体ｆｌの流速が限界流速Ｕｃｒ（ｉ’）に等しい時に作用する励振力Ｆ_ｅｘに対応する臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）は、伝熱管群１０内における伝熱管６と振止め部材１２との間の摩擦減衰に対応する正の減衰比ζｐ（ｉ’）と均衡している。従って、安定限界時における流速に対する有効流速Ｕｅの比の２乗を上記関係式にさらに乗じることで、有効流速Ｕｅが安定限界時における流速に近づくにつれて増大する励振力Ｆ_ｅｘに対応して、負減衰比ζｎ（ｉ’）を求めるための式を得ることができる。

続いて、図８のフローチャートの実行はステップＳ８４に進み、時刻歴応答解析部２１７は、伝熱管群１０の自励振動を模擬する演算に負減衰比ζｎ（ｉ’）の具体的な値を入力して実行する。具体的には、ステップＳ８３において、時刻歴応答解析部２１７は、入力として受け取った負減衰比ζｎ（ｉ’）の具体的な値を用いて時刻歴応答解析を実行し、伝熱管群１０の振動振幅の大きさを算出する。時刻歴応答解析部２１７によって実行される上記時刻歴応答解析は、図４〜図６を用いて上述した時刻歴応答解析と同様の演算である。つまり、上記時刻歴応答解析は、負減衰比ζｎ（ｉ’）の値を入力として、負減衰比ζｎ（ｉ’）に対応する励振力Ｆ_ｅｘが伝熱管群１０に作用した場合の伝熱管群１０の振動振幅を算出するパラメトリック・スタディ演算である。

続いて、図８のフローチャートの実行はステップＳ８５に進む。ステップＳ８４では、自励振動評価部２１８は、負減衰比ζｎ（ｉ’）の具体的な値を用いて時刻歴応答解析部２１７が時刻歴応答解析を実行した結果に基づいて、伝熱管群１０の振動振幅が発散するか否かを判定する。そして、自励振動評価部２１８は、当該判定の結果に基づいて伝熱管群１０の自励振動を評価する。具体的には、まず、自励振動評価部２１８は、負減衰比ζｎ（ｉ’）の具体的な値を用いて時刻歴応答解析を実行した時刻歴応答解析部２１７から伝熱管群１０の振動振幅を受け取る。続いて、自励振動評価部２１８は、当該振動振幅が発散しているか否かの判定を行う。例えば、自励振動評価部２１８は、当該振動振幅が所定の閾値を上回るか否かを判定し、閾値を上回る場合には、当該振動振幅が発散していると判定するようにしてもよい。最後に、自励振動評価部２１８は、当該振動振幅が発散しているならば、与えられた有効流速Ｕｅに対応する励振力Ｆ_ｅｘが伝熱管群１０に作用した場合に伝熱管群１０が自励振動を起こすと予測し、図フローチャートの実行を終了する。

このように、図７および図８を用いて上述した実施形態では、図４〜図６を用いて上述した実施形態のように計算コストが高いパラメトリック・スタディ演算を多数回反復実行するのではなく、伝熱管群１０がある流速（有効流速Ｕｅ）で自励振動（不安定振動）を起こすかどうかを確認するために実施することができる。具体的には、この実施形態では、ある流速Ｕｅが限界流速だと仮定すると、そのときの減衰比ζｎ（ｉ’）を上記式（６）から逆算し、その減衰比ζｎ（ｉ’）を時刻歴応答解析（パラメトリック・スタディ演算）に入力として与え、伝熱管群１０の振動振幅が発散しなければ流速Ｕｅで安定であることを確認することができる。これを別の観点から説明すると、この実施形態では、実際の運転条件に見合った有効流速Ｕｅを仮の限界流速であると仮定した上で、その運転条件（流速Ｕｅ）での負減衰比ζｎ（ｉ’）を求めている。そして、その時の運転条件（流速Ｕｅ）でその振動系が実際には安定であれば、その時の負減衰比ζｎ（ｉ’）は、本来よりも過大に見積もった負減衰比である。従って、その時の負減衰比ζｎ（ｉ’）をパラメトリック・スタディ演算に入力として与えても振動振幅が発散しないので、その時の運転条件（流速Ｕｅ）では自励振動を起こさないことの確認ができる。

以上より、図７および図８を用いて上述した実施形態では、取得された想定流速を仮の限界流速と仮定し、仮定された当該仮の限界流速に対応する負減衰比ζｎ（ｉ’）を入力して伝熱管群１０の自励振動を模擬する演算を実行した際に、伝熱管群１０の振動振幅が発散するか否かに基づいて自励振動が生じるか否かを判定する。言い換えると、この実施形態では、仮の限界流速に基づいて伝熱管群１０の自励振動を模擬する演算を実行した際に、伝熱管群１０の振動振幅が発散するか否かに応じて、当該仮の限界流速が真の限界流速Ｕｃｒ（ｉ’）を上回るか否かを検査している。従って、この実施形態によれば、仮の限界流速として仮定された流速で流体を流した場合に、伝熱管群１０の自励振動が実際に起きるか否かを、伝熱管群１０の自励振動を模擬する模擬演算によって高精度に予測することができる。

この実施形態では、伝熱管群１０に含まれる各伝熱管６に作用する流体ｆｌの動圧、各伝熱管６の質量密度または各伝熱管６の振動振幅が各伝熱管６の長さ方向ｙに沿って異なる場合、上述した動圧、密度または振動振幅の長さ方向ｙに沿った分布に基づいて流体ｆｌの有効流速Ｕｅを算出する。その上で、この実施形態では、有効流速Ｕｅを仮の限界流速と仮定して負減衰比ζｎ（ｉ’）を算出する。従って、この実施形態によれば、伝熱管群１０に含まれる各伝熱管６に作用する流体ｆｌの動圧、各伝熱管６の質量密度または各伝熱管６の振動振幅が各伝熱管６の長さ方向ｙに沿って異なる場合でも、伝熱管６内の場所による流速の違いを考慮して負減衰比ζｎ（ｉ’）を算出するための単一の流速値を求めることができる。

また、さらに別の実施形態では、時刻歴応答解析（パラメトリック・スタディ演算）に入力として与える負減衰比ζｎ（ｉ’）を有効流速Ｕｅから算出するのではなく、以下のように算出してもよい。すなわち、この実施形態では、伝熱管群１０を自由振動系で表したモデルの複素固有値分解によって、伝熱管群１０内において伝熱管６と振止め部材１２との間の摩擦力によって生じる摩擦減衰に対応する正の減衰比ζｐ（ｉ）を複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について算出する。そして、伝熱管群１０が構造上持っている摩擦減衰効果に対応する正の減衰比ζｐ（ｉ）の絶対値から負減衰比ζｎ（ｉ）を複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について求めている。何故なら、流体ｆｌの流速が限界流速Ｕｃｒ（ｉ’）に等しい時に作用する励振力Ｆ_ｅｘに対応する臨界負減衰比ζ_ｎ ^ｃｒ（ｉ’）は、伝熱管群１０内における伝熱管６と振止め部材１２との間の摩擦減衰に対応する正の減衰比ζｐ（ｉ’）と均衡しているからである。

具体的には、この実施形態では、伝熱管群１０に作用する励振力Ｆ_ｅｘに対応する外力項を０とした自由振動系として振動解析モデルを構築している。同時に、この実施形態では、伝熱管群１０を取り巻く流体および伝熱管６が振止め部材１２から摩擦力を受けることによって生じる摩擦減衰効果に対応して伝熱管群１０の振動特性に付加される付加質量、付加減衰および付加剛性を上記振動解析モデルに含めている。つまり、この実施形態では、伝熱管群１０が、ポアソン方程式で規定される圧力場とナビエ・ストークスの方程式（Ｎ−Ｓ式）で規定される流速場から成る流力場内に配置された場合に、伝熱管群１０が流体ｆｌから受ける励振力Ｆ_ｅｘを直接考慮せずに負減衰比ζｎ（ｉ’）を求める。

続いて、この実施形態では、複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について上記のように求めた負減衰比ζｎ（ｉ）の具体的な値を用いて時刻歴応答解析を実行し、伝熱管群１０の振動振幅の大きさを算出する。上記時刻歴応答解析は、図４〜図６を用いて上述した時刻歴応答解析と同様の演算である。つまり、上記時刻歴応答解析は、負減衰比ζｎ（ｉ’）の値を入力として、負減衰比ζｎ（ｉ’）に対応する励振力Ｆ_ｅｘが伝熱管群１０に作用した場合の伝熱管群１０の振動振幅を算出するパラメトリック・スタディ演算である。

一例において、この実施形態では、例えば、以下のようにして、伝熱管群１０を自由振動系で表したモデルの複素固有値分解によって、負減衰比ζｎ（ｉ）を複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々について算出してもよい。まず、上記式（２）の運動方程式を状態空間表現に書き換えて、以下の式に変形する。

そして、上記式（７）によって定義される一般固有値問題を解くと、上記式（２）の運動方程式で表される伝熱管群１０の振動特性について複数の固有値λ（ｉ）（１≦ｉ≦Ｉ）が複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々に対応して得られる。続いて、以下の式により、複数の固有値λ（ｉ）（１≦ｉ≦Ｉ）から複数の固有モードφ（ｉ）（１≦ｉ≦Ｉ）の各々に対応する負減衰比ζｎ（ｉ）と固有周波数ω（ｉ）が得られる。

以上より、この実施形態では、伝熱管群１０に作用する励振力Ｆ_ｅｘに対応する外力項を０とした自由振動系として振動解析モデルを構築することにより、伝熱管群１０に対して外力として作用する流力場の詳細な解析を行わなくても、当該自由振動系の運動方程式を複素固有値分解するだけで、時刻歴応答解析に入力すべき負減衰比ζｎ（ｉ）の適切な値を算出することができる。

３ 伝熱管
４ 第１直管部
５ 第２直管部
６（６ａ１，６ａ２，６ａ３，６ｂ１，６ｃ１） 曲り部
７ 管支持板
８ 管列
１０ 伝熱管群
１０ａ ベンド部
１１ 第１保持部材
１２ 止め部材
１２ａ 端部
１４ 第２保持部材
２０ コンピュータ装置
２１ 演算部
２２ 記憶部
２２ａ プログラム
２２ｂ データ
２３ 出力部
２４ 入力部
２１１ 限界流速算出部
２１２，２１７ 時刻歴応答解析部
２１３，２１８ 自励振動評価部
２１５ 有効流速算出部
２１６ 負減衰比算出部
Ｄ１ 面内方向
Ｄ２ 面外方向
Ｅｒｄ エネルギー消散量
Ｆｅｘ 関数
Ｆｅｘ 励振力
Ｈ 振動解析モデル
Ｕ，Ｕｅ 流速
Ｕｃｒ 限界流速
Ｕｅ 有効流速
ｄ１，ｄ２ 列方向
ｆ 固有振動数
ｆｌ 流体
ｍ 伝熱管の質量密度
ｘ 変位ベクトル
ｙ 長さ方向位置

Claims (9)

1. 流体中に支持部材に支持されて配置される管群の自励振動をコンピュータ装置を用いて評価するための自励振動評価方法であって、前記管群の一つ以上の固有モードの各々について、
   前記コンピュータ装置が備える時刻歴応答解析部によって前記流体の励振力に対応する負減衰比を変化させながら、前記管群の振動振幅の変化を模擬する時刻歴応答解析を実行する時刻歴応答解析ステップと、
   前記コンピュータ装置が備える限界流速算出部によって、前記時刻歴応答解析部で実行された前記時刻歴応答解析において前記管群の振動振幅の変化が発散する最小の負減衰比に基づいて、前記流体の限界流速を算出する限界流速算出ステップと、
   前記コンピュータ装置が備える入力部によって、前記流体の想定流速を入力する入力ステップと、
   前記コンピュータ装置が備える自励振動評価部によって、前記流体の想定流速を前記限界流速算出部で算出された前記限界流速と比較することにより、前記固有モードごとに前記管群の自励振動を評価する評価ステップと、
   を備える、自励振動評価方法。
2. 前記時刻歴応答解析部は、前記管群の振動解析モデルに対して前記負減衰比に対応する励振力を外力項として作用させた際に生じる振動振幅を時系列的に模擬し、
   前記振動解析モデルは、前記管群と前記支持部材との間に作用する接触荷重の分布を仮定することによって、前記管群と前記支持部材との間の摩擦力の大きさを規定する、
   ことを特徴とする、請求項１記載の自励振動評価方法。
3. 前記時刻歴応答解析部は、
   前記管群と前記支持部材との間の摩擦力に応じて消散される前記自励振動のエネルギー消散量に対応する第１減衰比と前記負減衰比とに基づいて前記管群の実効減衰比を算出し、
   前記算出された実効減衰比に基づいて、前記管群の振動振幅を時系列的に推定する
   ことを特徴とする、請求項１または請求項２に記載の自励振動評価方法。
4. 前記時刻歴応答解析部は、
   前記管群の振動振幅の変化に伴って、前記負減衰比が前記第１減衰比と一致した時点において前記振動振幅が発散すると判定する、
   ことを特徴とする、請求項３記載の自励振動評価方法。
5. 流体中に支持部材に支持されて配置される管群の自励振動をコンピュータ装置を用いて評価するための自励振動評価方法であって、
   前記コンピュータ装置が備える有効流速算出部によって前記流体の想定流速を取得する想定流速取得ステップと、
   前記コンピュータ装置が備える負減衰比算出部によって、前記有効流速算出部で取得された前記想定流速を限界流速と仮定し、前記限界流速と前記管群全体の負減衰比との間の相関に基づいて、前記想定流速に対応する負減衰比を算出する負減衰比算出ステップと、
   前記コンピュータ装置が備える自励振動評価部によって、前記管群の自励振動を模擬する演算に前記負減衰比算出部で算出された前記負減衰比を入力して実行した際に、前記管群の振動振幅が発散するか否かに基づいて前記管群の自励振動を評価する評価ステップと、
   を備える、自励振動評価方法。
6. 前記想定流速取得ステップは、
   前記管群に含まれる各管に作用する前記流体の動圧、前記各管の密度または前記各管の振幅の少なくとも１つについて前記各管の長さ方向に沿った分布に基づいて前記流体の有効流速を算出する有効流速算出ステップを含み、
   前記負減衰比算出ステップでは、前記有効流速を前記想定流速として前記負減衰比を算出する、
   ことを特徴とする、請求項５記載の自励振動評価方法。
7. 前記管群は、同一の平面内において延在するとともに、互いに曲率中心を共有し、且つ、互いに曲率半径が異なる曲り部を有する複数のＵ字管によって形成される管列を少なくとも一つ含み、
   前記支持部材は、前記管列を挟むように前記平面に沿って延在する形で前記管列の両側に配置された
   少なくとも一対の振止め部材を含み、
   前記管群を通る流体の励振力に抗して、振止め部材との摩擦力により支持される前記管群の前記平面に沿った方向における流力弾性振動の安定判別を行う、
   ことを特徴とする、請求項１乃至請求項６の何れか一項に記載の自励振動評価方法。
8. 前記管群は、ＰＷＲ原子力発電設備の蒸気発生器の伝熱管群であることを特徴とする請求項１乃至７の何れか一項に記載の自励振動評価方法。
9. 流体中に支持部材に支持されて配置される管群の自励振動をコンピュータ装置を用いて評価するための自励振動評価プログラムであって、前記コンピュータ装置に、前記管群の一つ以上の固有モードの各々について、
   前記コンピュータ装置が備える時刻歴応答解析部によって前記流体の励振力に対応する負減衰比を変化させながら、前記管群の振動振幅の変化を模擬する時刻歴応答解析を実行する時刻歴応答解析ステップと、
   前記コンピュータ装置が備える限界流速算出部によって、前記時刻歴応答解析部で実行された前記時刻歴応答解析において前記管群の振動振幅の変化が発散する最小の負減衰比に基づいて、前記流体の限界流速を算出する限界流速算出ステップと、
   前記コンピュータ装置が備える入力部によって、前記流体の想定流速を入力する入力ステップと、
   前記コンピュータ装置が備える自励振動評価部によって、前記流体の想定流速を前記限界流速算出部で算出された前記限界流速と比較することにより、前記固有モードごとに前記管群の自励振動を評価する評価ステップと、
   を実行させるための自励振動評価プログラム。
   

Images