(57)【要約】【課題】回転三角配列に配置された管群について、多くの時間とコストを要することなく、多数のパラメータについて臨界流速を予測する。【解決手段】管群１０の軸方向に直交する面内において、各管１１１〜１１１１をそれぞれ含む複数のセル１２１〜１２１１に分割するステップと、各セルの外周において流体の流速及び圧力をモード展開した外周のモード座標を用いて外周の少なくとも一部の流体の流速及び圧力を表すステップと、管の中心を始点としてセルに含まれる管との境界をなす内周と外周を結ぶ動径について、外周のモード座標と、流体の流速及び圧力を周方向にモード展開した周方向のモード座標とを用いて、各セル内の流体の流速及び圧力を表すステップと、を含み、各セルは、隣接する管の中心を結ぶ直線の少なくとも１つを用いて分割したものである。【選択図】図２
2015026259.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いて、回転三角配列に配置された所定方向に延びる複数の所定径の管を含む管群の軸方向に直交して流れる流体の流速及び圧力を予測する方法であって、
前記管群の軸方向に直交する面内において、各管をそれぞれ含む複数の領域に分割するステップと、
各領域の外周において流体の流速及び圧力をモード展開した外周のモード座標を用いて外周の少なくとも一部の流体の流速及び圧力を表すステップと、
管の中心を始点として前記領域に含まれる管との境界をなす内周と外周を結ぶ動径について、前記外周のモード座標と、流体の流速及び圧力を周方向にモード展開した周方向のモード座標とを用いて、各領域内の流体の流速及び圧力を表すステップと、を含み、
各領域は、隣接する管の中心を結ぶ直線の少なくとも１つを用いて分割したものであることを特徴とする管群振動予測方法。
【請求項２】
前記周方向のモード座標は、前記面内において管の中心を軸とした極座標
【数１】
2015026259.tif 000102
について、管の半径をａとして、径方向の
【数２】
2015026259.tif 000103
及び周方向の
【数３】
2015026259.tif 000104
において、時刻ｔにおける流速の径方向及び周方向の成分、圧力
【数４】
2015026259.tif 000105
を各領域の外周における流速の径方向及び周方向の成分、圧力
【数５】
2015026259.tif 000106
が寄与する第１項
【数６】
2015026259.tif 000107
と、整数ｅを領域の番号、整数ｍ及びｋを周方向及び径方向のモード次数、
【数７】
2015026259.tif 000108
として、各領域の内部における流速の径方向及び周方向の成分、圧力を直交基底で展開した第２項
【数８】
2015026259.tif 000109
との和として表したときに、第２項の展開係数
【数９】
2015026259.tif 000110
として与えられることを特徴とする請求項１記載の管群振動予測方法。
【請求項３】
各領域は、前記管群を含む領域について各管を含むように分割した正六角形の少なくとも一部を含み、前記外周のモード座標は、正六角形の外周において流体の流速及び圧力をモード展開したものであることを特徴とする請求項１又は２に記載の管群振動予測方法。
【請求項４】
前記外周のモード座標は、前記面内における直交座標
【数１０】
2015026259.tif 000111
について、整数ｉ，ｊを正六角形の辺の頂点又は中点の番号、
【数１１】
2015026259.tif 000112
をこれら頂点又は中点における座標として、時刻ｔにおける辺上での流速のｘ，ｙ成分、圧力
【数１２】
2015026259.tif 000113
を頂点又は中点における流速のｘ，ｙ成分、圧力
【数１３】
2015026259.tif 000114
が辺上の位置
【数１４】
2015026259.tif 000115
に寄与する第１項及び第２項
【数１５】
2015026259.tif 000116
と、整数ｃを辺の番号、整数ｎをモード次数として、辺上における流速のｘ，ｙ成分及び圧力を辺上の直交基底で展開した第３項
【数１６】
2015026259.tif 000117
との和として表したときに、第３項の展開係数
【数１７】
2015026259.tif 000118
として与えられることを特徴とする請求項２に記載の管群振動予測方法。
【請求項５】
時刻ｔにおける流速の径方向及び周方向の成分、圧力
【数１８】
2015026259.tif 000119
の外周における値
【数１９】
2015026259.tif 000120
時刻ｔにおける外周での流速のｘ，ｙ成分、圧力
【数２０】
2015026259.tif 000121
は、方向余弦を用いて
【数２１】
2015026259.tif 000122
なる関係を有することを特徴とする請求項４に記載の管群振動予測方法。
【請求項６】
ナビエ・ストークスの式及び連続条件に基づいて、前記モード座標を一般化座標として、前記流速及び圧力の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記流体の定常流による流れ場を決定することを特徴とする請求項１〜５のいずれか一項に記載の管群振動予測方法。
【請求項７】
流体の流速及び圧力にそれぞれ振動外乱を重畳し、これら流速、圧力及び振動外乱についてナビエ・ストークスの式、連続条件及び前記定常流による流れ場に基づいて、前記振動外乱の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記管群の振動を予測することを特徴とする請求項６に記載の管群振動予測方法。
【請求項８】
前記振動外乱の一般化座標の時間に関する常微分方程式に基づいて、前記管群の不安定振動が発生する前記流体の臨界流速を決定することを特徴とする請求項７に記載の管群振動予測方法。
【請求項９】
前記常微分方程式は、ナビエ・ストークスの式及び連続条件に基づいて、ガレルキン法を用いて求めることを特徴とする請求項６〜８のいずれか一項に記載の管群振動予測方法。
【請求項１０】
前記管群は加圧水型原子炉の蒸気発生器における伝熱管群を構成し、前記流体は前記蒸気発生器における二次冷却水であることを特徴とする請求項１〜９のいずれか一項に記載の管群振動予測方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、回転三角配列に配置された管群に直交して流れる流体により発生する管群の振動を予測する管群振動予測方法に関する。
【背景技術】
【０００２】
従来、加圧水型原子炉（ＰＷＲ）に備えられる蒸気発生器（ＳＧ）においては、原子炉から供給される一次冷却水を流す伝熱管を並列配置して管群とし、この管群を伝熱面としてこの管群の軸方向と直交する向きに二次冷却水を流して熱交換している。この伝熱管の管群の配置には、回転三角配列、ノーマル三角配列、正方配列等が使用されている。この内、回転三角配列は、高い熱交換の効率が期待できるので広範に使用されている。
【０００３】
蒸気発生器においては、熱交換の効率化を図るため、管群に直交した流れの流速を高めることが求められている。しかしながら、この直交流の流速がある値に達すると、不安定振動が発生することが知られている（下記非特許文献１参照）。原子炉の安全性を確保するためにはこのような不安定振動を制御することが必要であり、不安定振動を引き起こす臨界流速を予測する研究がなされてきた。
【０００４】
この臨界流速の予測には、ＣＦＤ（数値流体力学）による数値シミュレーション（下記非特許文献２参照）が提供されている。また、１本の管の振動によって隣接する複数の管に作用する力を実験により求め、この実験データに基づいて不安定振動解析モデルを作成する方法が提供されている（下記非特許文献３参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】日本機械学会編、事例に学ぶ流体関連振動、ｐ．６６
【非特許文献２】市岡丈彦、川田裕、中村友道、泉元、小林敏雄、高松洋、「流動解析を用いた管群の流力弾性振動に関する研究（第１報、並列２円柱および１列管群の解析）」、日本機械学会論文集（Ｂ編）、１９９５年、６１巻、５８２号、ｐ．１４５−１５１
【非特許文献３】田中博喜、「円柱群の流力弾性振動に関する研究（流水中の円柱群）」、昭和５６年、日本機械学会論文集（Ｂ編）、４７巻、４１９号、ｐ．１２２４−１２３３
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、ＣＦＤによるシミュレーションは、モデルの入力、計算処理等に多くの時間とコストを要し、多数のパラメータについて臨界流速を予測することは困難であった。また、実験データに基づいて不安定振動解析モデルを作成する方法についても、多数のパラメータについて臨界流速を予測することは困難であった。
【０００７】
この出願に係る発明は、上述の実情に鑑みて提案されるものであって、三角配列で配置された管群に対する直交流について、多くの時間とコストを要することなく、多数のパラメータについて臨界流速を予測することができるような管群振動予測方法を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上述の課題を解決するために、本発明に係る管群振動予測方法の発明は、コンピュータを用いて、回転三角配列に配置された所定方向に延びる複数の所定径の管を含む管群の軸方向に直交して流れる流体の流速及び圧力を予測する方法であって、前記管群の軸方向に直交する面内において、各管をそれぞれ含む複数の領域に分割するステップと、各領域の外周において流体の流速及び圧力をモード展開した外周のモード座標を用いて外周の少なくとも一部の流体の流速及び圧力を表すステップと、管の中心を始点として前記領域に含まれる管との境界をなす内周と外周を結ぶ動径について、前記外周のモード座標と、流体の流速及び圧力を周方向にモード展開した周方向のモード座標とを用いて、各領域内の流体の流速及び圧力を表すステップと、を含み、各領域は、隣接する管の中心を結ぶ直線の少なくとも１つを用いて分割したものである。
【０００９】
前記周方向のモード座標は、前記面内において管の中心を軸とした極座標
【数１】
2015026259.tif 000003
【００１０】
について、管の半径をａとして、径方向の
【数２】
2015026259.tif 000004
【００１１】
及び周方向の
【数３】
2015026259.tif 000005
【００１２】
において、時刻ｔにおける流速の径方向及び周方向の成分、圧力
【数４】
2015026259.tif 000006
【００１３】
を各領域の外周における流速の径方向及び周方向の成分、圧力
【数５】
2015026259.tif 000007
【００１４】
が寄与する第１項
【数６】
2015026259.tif 000008
【００１５】
と、整数ｅを領域の番号、整数ｍ及びｋを周方向及び径方向のモード次数、
【数７】
2015026259.tif 000009
【００１６】
として、各領域の内部における流速の径方向及び周方向の成分、圧力を直交基底で展開した第２項
【数８】
2015026259.tif 000010
【００１７】
との和として表したときに、第２項の展開係数
【数９】
2015026259.tif 000011
【００１８】
として与えられることが好ましい。
【００１９】
各領域は、前記管群を含む領域について各管を含むように分割した正六角形の少なくとも一部を含み、前記外周のモード座標は、正六角形の外周において流体の流速及び圧力をモード展開したものであることが好ましい。
【００２０】
前記外周のモード座標は、前記面内における直交座標
【数１０】
2015026259.tif 000012
【００２１】
について、整数ｉ，ｊを正六角形の辺の頂点又は中点の番号、
【数１１】
2015026259.tif 000013
【００２２】
をこれら頂点又は中点における座標として、時刻ｔにおける辺上での流速のｘ，ｙ成分、圧力
【数１２】
2015026259.tif 000014
【００２３】
を頂点又は中点における流速のｘ，ｙ成分、圧力
【数１３】
2015026259.tif 000015
【００２４】
が辺上の位置
【数１４】
2015026259.tif 000016
【００２５】
に寄与する第１項及び第２項
【数１５】
2015026259.tif 000017
【００２６】
と、整数ｃを辺の番号、整数ｎをモード次数として、辺上における流速のｘ，ｙ成分及び圧力を辺上の直交基底で展開した第３項
【数１６】
2015026259.tif 000018
【００２７】
との和として表したときに、第３項の展開係数
【数１７】
2015026259.tif 000019
【００２８】
として与えられることが好ましい。
【００２９】
時刻ｔにおける流速の径方向及び周方向の成分、圧力
【数１８】
2015026259.tif 000020
【００３０】
の外周における値
【数１９】
2015026259.tif 000021
【００３１】
時刻ｔにおける外周での流速のｘ，ｙ成分、圧力
【数２０】
2015026259.tif 000022
【００３２】
は、方向余弦を用いて
【数２１】
2015026259.tif 000023
【００３３】
なる関係を有することが好ましい。
【００３４】
ナビエ・ストークスの式及び連続条件に基づいて、前記モード座標を一般化座標として、前記流速及び圧力の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記流体の定常流による流れ場を決定することが好ましい。
【００３５】
流体の流速及び圧力にそれぞれ振動外乱を重畳し、これら流速、圧力及び振動外乱についてナビエ・ストークスの式、連続条件及び前記定常流による流れ場に基づいて、前記振動外乱の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記管群の振動を予測することが好ましい。
【００３６】
前記振動外乱の一般化座標の時間に関する常微分方程式に基づいて、前記管群の不安定振動が発生する前記流体の臨界流速を決定することが好ましい。
【００３７】
前記常微分方程式は、ナビエ・ストークスの式及び連続条件に基づいて、ガレルキン法を用いて求めることが好ましい。
【００３８】
前記管群は加圧水型原子炉の蒸気発生器における伝熱管群を構成し、前記流体は前記蒸気発生器における二次冷却水であることが好ましい。
【発明の効果】
【００３９】
この発明は、回転三角配列により配置された管群の正六角形のセルについて、モード展開により解析的に得られた常微分方程式を使用しているので、多くの時間とコストを要することなく、多数のパラメータについて臨界流速を予測することができる。
【図面の簡単な説明】
【００４０】
【図１】管群振動予測方法の一連の流れを示す図である。
【図２】回転三角配列の管群を示す断面図である。
【図３】１セルにおける座標設定を示す図である。
【図４】三角配列の数値例題を示す図である。
【図５】換算流速の臨界値の質量減衰パラメータ依存性を示す図である。
【図６】正方配列の管群を示す断面図である。
【図７】管群振動を予測する計算装置の構成を示す図である。
【図８】計算装置における一連の動作の流れを示す図である。
【図９】粘性力の導出を示す図である。
【発明を実施するための形態】
【００４１】
以下、本発明に係る管群振動予測方法について、図面を参照して詳細に説明する。
【００４２】
〔１．本実施の形態の概要〕
図１は、本実施の形態の管群振動予測方法の一連の流れを示す図である。最初のステップＳ１１では、管群の軸の方向に直交する平面を管の属する正六角形領域に分割する。この領域を以下、セルと称する。ステップＳ１２では、セル外周（各辺）での未知量（流速成分、圧力）を、セルの各頂点での値の１次補間と補間誤差を補うモード展開式の和で表す。ステップＳ１３では、セルの外周と内周（管周囲）を結ぶ任意周方向座標の動径上で、任意の動径座標での未知量を、ステップＳ１２で設定した外周での解に適合する１次式とモード展開式の和で表す。ステップＳ１４では、ステップＳ１２〜Ｓ１３によって得た解の表示式を流体の支配方程式系に代入して、変分直接法（ガレルキン法）で時間の常微分方程式系に離散化し、管群の運動方程式と連成させて解く。
【００４３】
ステップＳ１４に続いて、次の双方のステップ（Ａ）、（Ｂ）をさらに実行することもできる。ここで、ステップ（Ａ）はステップ（Ｂ）の前段階として必要である。
【００４４】
（Ａ）各管静止条件下での定常（平均）流れ場の決定（流れ解析）
（Ｂ）振動外乱に関する安定性解析
本実施の形態では、セルを正六角形の形状にすることにより、回転三角配列に配置した管群に適用するものである。ここで、回転三角配列とは、三角配列の内、管群の軸方向に直交する平面内において隣接する管を結ぶ直線の一つの方向に流体が流れるものをいう。
【００４５】
〔２．解析方法〕
図２のように流れを受ける管群１０について、不安定振動が発生する臨界流速を求める。この管群１０は、それぞれ半径ａの第１の管１１１から第１１の管１１１１までの１１本の管１１１〜１１１１が、軸間の間隔（ピッチ）をＴとして、各管１１１〜１１１１の軸が延びるｚ方向について回転三角配列に配置され、Ｘ方向に沿って３列に並んでいる。以下では簡単のため、各管１１１〜１１１１を特定する添字１〜１１を省略する。
【００４６】
管１１の軸方向に流速、圧力、管の変位の変化しない２次元問題を考え、この軸に直交する２次元平面を想定する。この面内において、最初に各管１１の静止条件下での定常（平均）流れ場を決定し（２．１節）、次に振動外乱を重畳させて安定性解析を行う（２．２節）。
【００４７】
〔２．１  流れ場の決定〕
〔２．１．１  解の表示〕
図３において、各管１１を含む正六角形領域である各セル１２に関して座標系
【数２２】
2015026259.tif 000024
【００４８】
を定義する。
【００４９】
また、セル１２の頂点又は中点の番号ｉを定義し、頂点又は中点ｉでのｘ，ｙ方向の流速、圧力をそれぞれ
【数２３】
2015026259.tif 000025
【００５０】
とする。さらに、各外周辺に沿って、座標ｓをとり、各外周辺上でのｘ，ｙ方向の流速、圧力をそれぞれ
【数２４】
2015026259.tif 000026
【００５１】
とする。
【００５２】
これら各外周辺上での流速、圧力は、辺の頂点又は中点ｉ，ｊでの値の直線補間項に精度向上のための外周辺上の直交基底についてのモード展開項を付加して、次のように表すことができる。ここで、ｃは辺の番号を表している。
【数２５】
2015026259.tif 000027
【００５３】
ここで
【数２６】
2015026259.tif 000028
【００５４】
である。
【００５５】
また、
【数２７】
2015026259.tif 000029
【００５６】
である。
【００５７】
流体運動解析は、各セル１２に関する極座標を用いて行う。このため、式（１）の流速成分を、次のように極座標に変換する。
【数２８】
2015026259.tif 000030
【００５８】
ここで、方向余弦は
【数２９】
2015026259.tif 000031
【００５９】
である。
【００６０】
周方向座標
【数３０】
2015026259.tif 000032
【００６１】
での内周と外周を結ぶ動径上での流速、圧力を、次のように表す。
【数３１】
2015026259.tif 000033
【００６２】
ここで
【数３２】
2015026259.tif 000034
【００６３】
である。
【００６４】
式（５）右辺第１項は、内周ｒ＝ａ、外周ｒ＝ｒｍａｘでそれぞれ０，１となり、外周で解（３）に等しくなるための項である。式（５）右辺第２項は、精度向上のための直交関数展開項で、
【数３３】
2015026259.tif 000035
【００６５】
である。このような一般化座標は、モード次数で展開したものであり、モード座標と称することもある。内周での流速は、境界層が非常に薄いため、法線方向成分のみが０となる境界条件を満たすように設定している。
【００６６】
本実施の形態の解析法は、
（１）式（１）のｉ，ｊに相当する点をセルの頂点だけでなく辺の中点にも設定でき、周方向角
【数３４】
2015026259.tif 000036
【００６７】
の区間を各管ごとに変えられる、
（２）周方向モード関数を流速成分、圧力
【数３５】
2015026259.tif 000037
【００６８】
について個別に設定できる、という２つの条件を満たすものである。
【００６９】
このような条件（１），（２）を満たすことにより、
【数３６】
2015026259.tif 000038
【００７０】
について、例えば
【数３７】
2015026259.tif 000039
【００７１】
のように、回転三角配列に適用することを可能にしている。
【００７２】
式（３）に式（１）を代入した結果を、式（５）右辺第１項に代入する。これによって、各セル１２内の任意位置での流速と圧力が、頂点又は中点でのｘ，ｙ方向の流速と圧力、式（１）右辺第３項のモード展開係数（モード座標）、式（５）右辺第２項のモード展開係数（モード座標）とで表される。すなわち、流体の無限自由度が、これらの時間変数に縮小する。そこで、このようにして得た式を自由度縮小式と称する。
【００７３】
〔２．１．２  ガレルキン法〕
自由度縮小式を、下記の連続条件とナビエ・ストークス方程式の変分重み付き表示式に代入して、未知の時間変数に関する常微分方程式を導く。
【数３８】
2015026259.tif 000040
【数３９】
2015026259.tif 000041
【数４０】
2015026259.tif 000042
【００７４】
ここで
【数４１】
2015026259.tif 000043
【００７５】
である。
【００７６】
常微分方程式は、次のマトリックス方程式の形に導く。
【数４２】
2015026259.tif 000044
【００７７】
ここで
【数４３】
2015026259.tif 000045
【００７８】
である。
【００７９】
式（１０）を時間積分し、時間変動の小さい定常近似解に達した解によって流れ場を決定する。このとき数値計算の安定化を行った。多用されている人工粘性は流れ場を変えてしまうので使わず、運動方程式の時間微分を陽に含む非定常慣性項について係数ρｆを
【数４４】
2015026259.tif 000046
【００８０】
とする人工慣性を導入した。粘性項と異なり非定常慣性項は定常状態で０に近づくため、求めるべき定常近似解への影響は小さいことに留意した工夫である。
【００８１】
〔２．２  振動解析〕
〔２．２．１  解の表示〕
解を次のように設定する。
【数４５】
2015026259.tif 000047
【００８２】
これは、前節で求めた定常成分
【数４６】
2015026259.tif 000048
【００８３】
に振動外乱
【数４７】
2015026259.tif 000049
【００８４】
が重畳した形である。
【００８５】
２．１．１節の自由度低減手法を、振動成分に関する方程式系の導出にも用いる。まず、セル１２の外周辺上でのｘ，ｙ方向の流速、圧力に関する振動成分を、頂点又は中点ｉ，ｊでの値
【数４８】
2015026259.tif 000050
【００８６】
の直線補間と、モード展開項の和で表し、式（１）に対応して次式を得る。
【数４９】
2015026259.tif 000051
【００８７】
流体運動解析は、図２のような、各セル１２に関する極座標を用いて行うため、流速成分を極座標に変換し、周方向座標
【数５０】
2015026259.tif 000052
【００８８】
での内周と外周を結ぶ動径上での流速、圧力（振動成分）を、式（６）に対応して次のように表す。
【数５１】
2015026259.tif 000053
【００８９】
ここで
【数５２】
2015026259.tif 000054
【００９０】
である。
【００９１】
このようにして、各セル１２内の任意位置での流速と圧力（振動成分）を、頂点又は中点でのｘ，ｙ方向の流速と圧力、式（１２）右辺第３項のモード展開係数、式（１４）右辺第２項のモード展開係数で表す。すなわち、流体の無限自由度を、これらの時間変数に縮小する。
【００９２】
〔２．２．２  ガレルキン法〕
ナビエ・ストークス方程式は、変分重み付きの形で次のように表される。
【数５３】
2015026259.tif 000055
【数５４】
2015026259.tif 000056
【００９３】
運動学的条件は、管１１表面のｘ，ｙ方向の速度（ｄ／ｄｔ）ｑｘ，（ｄ／ｄｔ）ｑｙの法線方向成分を管１１表面での湧き出しと考え、湧き出し項を流体連続条件式に導入することによって、次のようになる。
【数５５】
2015026259.tif 000057
【００９４】
各管の運動方程式は、次式によって与えられる。
【数５６】
2015026259.tif 000058
【００９５】
ここで
【数５７】
2015026259.tif 000059
【００９６】
である。
【００９７】
式（１４）を式（１５），（１６），（１７），（１８）に代入することにより、未知の時間変数に関する常微分方程式系を下記のマトリックス方程式の形に導く。
【数５８】
2015026259.tif 000060
【００９８】
ここで
【数５９】
2015026259.tif 000061
【００９９】
は一般化座標によって構成される列ベクトルであり、マトリックス
【数６０】
2015026259.tif 000062
【０１００】
の固有値を求めることによって安定性を判定する。
【０１０１】
〔３．計算結果〕
〔３．１  計算条件〕
次の緒元について計算を行った。
【数６１】
2015026259.tif 000063
【０１０２】
流れ場の決定では、粘性係数の代わりに乱流粘性係数をμの６０倍の値に設定した（ＣＦＤにおいては２０〜１００倍の値を用いることに基づく）。
【０１０３】
管群１０の問題では、管１１が多数あるため、図４の太い実線の直線で仕切られ、境界Ａ，Ｂにそれぞれ正、負の、大きさ
【数６２】
2015026259.tif 000064
【０１０４】
の湧き出しを与えた流体領域を縮小モデルとした。正および負の湧き出しの影響の少ない管Ｄの
【数６３】
2015026259.tif 000065
【０１０５】
での流れ場を図４の１１本の管の
【数６４】
2015026259.tif 000066
【０１０６】
の領域に設定し、これらの各管の
【数６５】
2015026259.tif 000067
【０１０７】
での流れ場は、
【数６６】
2015026259.tif 000068
【０１０８】
を結ぶ直線に関する対称性より決定した。
【０１０９】
〔３．２  臨界流速の従来データとの比較〕
質量減衰パラメータと呼ばれる次の量
【数６７】
2015026259.tif 000069
【０１１０】
を、減衰比ζbを変えることによって変え、換算流速（reduced velocity）と称されている下記の量の臨界値を求めた。
【数６８】
2015026259.tif 000070
【０１１１】
ただし、
【数６９】
2015026259.tif 000071
【０１１２】
である。
【０１１３】
管１１の間の隙間の流速Ｖｃは、流れ場を決めるときの式（９）中の湧き出し、吸込みソースの強さＶｓｒｃのＴ／（Ｔ−２ａ）倍である。臨界値を従来データ（非特許文献１）と比較した結果を図５に示す。図中の“△”は従来データであり、図中の曲線Ａは本解析により得られた値である。本解析により従来データに近い値が得られていることが確かめられる。
【０１１４】
上記は、流体密度の大きい液体の場合を対象とした、質量減衰パラメータの小さい範囲の結果である。しかし工学上は、蒸気発生器（ＳＧ）のように、ボイド率９０％、すなわち流体密度が１／１０に低下し、式（２０）の質量減衰パラメータがより大きい場合が重要になる。この場合の結果を図中の曲線Ｂで同様に示した。この場合、
【数７０】
2015026259.tif 000072
【０１１５】
である。
【０１１６】
〔３．３自由度低減効果〕
本解析では、自由度数がＣＦＤに比べると少なくて済むことを確かめる。解の表示式（１２），（１３）より、流体解析に関する自由度数Ｎは次式によって求められる。
【数７１】
2015026259.tif 000073
【０１１７】
ここで、Ｎ１は頂点の個数、Ｎ２は辺の個数、Ｎ１３は管の個数（セルの個数）であり、係数３は、ｘ，ｙ方向の流速成分と圧力とで３個の変数があることによる。図４で計算に用いた配列の場合について自由度数を求めると、Ｎ１＝３６，Ｎ２＝３６，ｎｍａｘ＝２，Ｎ３＝１１であり、ｍｍａｘ＝１，ｋｍａｘ＝１で上のように従来データを説明し得る解析結果が得られる。
【０１１８】
この場合の自由度数はＮ＝３９０である。ＣＦＤの数百万オーダに比べるとかなり少なくて済む。安定判別の固有値を求めるために要する計算時間は１ケースについて２分程度である（ただし臨界流速を求めるためには、固有値計算を１０回程度繰り返す必要がある）。このような低次元の固有値問題に帰着できることは、ＣＦＤ使用の際に流体と管１１の周波数差により流体解析に合わせて時間きざみを非常に小さくして時間積分する必要が生じる問題の回避に有効である。
【０１１９】
〔４．比較例〕
本実施の形態は回転三角配列に配置された管群に適用されるが、例えば正方配列により配置された管群に適用される場合には、対称性を考慮してモード座標の設定や流れの計算を簡略化することができる。ここでは、回転三角配列の場合との比較のため、正方配列により配置された管群への適用について簡単に説明する。
【０１２０】
図６は、正方配列により配置された管群１０を示す図である。この図においては、各管１１が正方配列により配置されている他は回転三角配列により配置された図２と同様の構成であるので、対応する符号を付して説明を省略する。
【０１２１】
このような正方配列の場合、内周と外周を結ぶ動径上での流速、圧力の式（５）は、次のように簡略化することができる。
【数７２】
2015026259.tif 000074
【０１２２】
ここで
【数７３】
2015026259.tif 000075
【０１２３】
である。
【０１２４】
そして、流れ場の計算は、図６において、３個の領域
【数７４】
2015026259.tif 000076
【０１２５】
での流れ場は同じである。そこで、１個の領域のみについて、流れ場を決定する計算を行うことができる。次に、Ｘ方向についても管１１が多数であるため周期対称な問題である。そこで、出入り口での湧き出し、吸込みの影響の少ない領域
【数７５】
2015026259.tif 000077
【０１２６】
で求まった流れ場を９本の管１１の周りのセル１２に設定して安定性解析を行うことができる。
【０１２７】
本実施の形態の回転三角配列の場合は、このような対称性を利用することができないので、式（５）のように周方向の角度に対する依存性を含む流速、圧力を導出する必要がある。
【０１２８】
〔５．計算装置〕
前述したような臨界流速の予測は、図７に示すような計算装置２０によって実現することができる。この計算装置２０は、ＣＰＵ、ＤＳＰの如き演算部２１、ＲＡＭ、ＲＯＭ、ハードディスクの如き記憶部２２、ＬＣＤ、プリンタの如き出力部２３、キーボード、マウスの如き入力部２４を含み、例えばパーソナルコンピュータを利用することができる。
【０１２９】
図８に示す計算装置２０の一連の動作は、記憶部２２に格納された臨界流速算定プログラム２２ａを演算部２１が読み出して実行することにより実現される。この臨界流速算定プログラム２２ａは、前述のような手順によって得られた流速及び圧力の一般化座標の時間に関する常微分方程式と振動外乱の一般化座標の時間に関する常微分方程式の表現を含んでいる。
【０１３０】
最初のステップＳ２１においては、モデルを設定する。ここでは、前述したような管群１０のモデルを設定するものとする。入力部２４は、このモデルについて、管１１の径ａ、管１１を配置した間隔Ｔ、管１１の単位長さ当たりの質量、管１１の固有振動数、流体の密度、流体の粘性係数の少なくとも１つの値を入力値として受け取る。演算部２１は、入力部２４が受け取った入力値を記憶部２２に格納する。
【０１３１】
ステップＳ２２においては、演算部２１は、記憶部２２に格納された入力値を読み出し、その数値計算部２１ａにおいて、この入力値に基づいて前記常微分方程式の表現を用いてこの管群における臨界流速について数値計算する。演算部２１は、得られた臨界流速の値を記憶部２２に格納する。
【０１３２】
ステップＳ２３においては、演算部２１は、記憶部２２に格納された臨界流速の値と、同じく記憶部２２に格納された所定の閾値２２ｂとを読み出す。演算部２１は、その判定部２１ｂにおいて、臨界流速値が閾値２２ｂを超えた場合にはＯＫと判定して一連のステップを終了する。一方、臨界流速が閾値を越えない場合にはＮＧとして判定して前のステップＳ２１のモデル設定に手順を戻す。なお、閾値２２ｂは、入力部２４を介して設定することができる。
【０１３３】
このような一連の工程において、臨界流速の値が所定の閾値を超えるまでモデル設定、数値計算、判定のループを繰り返すことにより、臨界流速が閾値を超えるようなモデル設定を可能としている。また、前述のような常微分方程式の表現を利用することにより、精度の高いモデル設定を可能としている。なお、これに限られず、臨界流速が閾値内に収まるようにモデル設定することもできる。
【０１３４】
なお、このような臨界流速の算定は、記憶部２２に格納した臨界流速算定プログラム２２ａのような、前述の常微分方程式の表現を含み、モデル設定、数値計算、判定のステップを有するプログラムによっても提供することができる。
【０１３５】
〔６．付録：式（１８）における粘性力の導出〕
図９に示すように液体とｚ＝０の平面３１で接する物体３２のｘ方向の速度が次式
【数７６】
2015026259.tif 000078
【０１３６】
のようにある周波数ωで振動する場合の境界層流れを考える。ｘ方向の流速をｖｘ、圧力をｐとし、これらの境界層のすぐ外側の主流での値を、上添字（ｍ）を付けて表わす。
【０１３７】
主流に関する運動方程式は
【数７７】
2015026259.tif 000079
【０１３８】
となる。
【０１３９】
この運動方程式の解を
【数７８】
2015026259.tif 000080
【０１４０】
とおき、式（Ａ２）に代入すると、次式が得られる。
【数７９】
2015026259.tif 000081
【０１４１】
境界層内流れを解析するため、ｚ方向の流速は非常に小さいと仮定すると、ｚ方向のナビエ・ストークス方程式は 
【数８０】
2015026259.tif 000082
【０１４２】
に帰着する。
【０１４３】
従って、条件
【数８１】
2015026259.tif 000083
【０１４４】
によって次式が要求される。
【数８２】
2015026259.tif 000084
【０１４５】
境界層は非常に薄いので、近似
【数８３】
2015026259.tif 000085
【０１４６】
をｘ方向のナビエ・ストークス方程式に適用できる。このようにして次式が得られる。
【数８４】
2015026259.tif 000086
【０１４７】
式（Ａ７）を下記の境界条件の下で解く。
【数８５】
2015026259.tif 000087
【数８６】
2015026259.tif 000088
【０１４８】
解を
【数８７】
2015026259.tif 000089
【０１４９】
と表わし、式（Ａ１０）と式（Ａ３）の第２式を式（Ａ７）〜（Ａ９）に代入すると、次のようになる。
【数８８】
2015026259.tif 000090
【数８９】
2015026259.tif 000091
【数９０】
2015026259.tif 000092
【０１５０】
式（Ａ４）を式（Ａ１１）に代入して次式を得る。
【数９１】
2015026259.tif 000093
【０１５１】
式（Ａ１４）に関する同次方程式の特性根αは、
【数９２】
2015026259.tif 000094
【０１５２】
によって定まり、解が有限である条件より、 
【数９３】
2015026259.tif 000095
【０１５３】
となる。
【０１５４】
従って、式（Ａ１４）の一般解は次のようになる。
【数９４】
2015026259.tif 000096
【０１５５】
Ｃｘ（ｘ）は任意関数で、式（Ａ１３）の条件を課すことにより、次のように物体速度と主流速度で表わせる。
【数９５】
2015026259.tif 000097
【０１５６】
式（Ａ１６）を式（Ａ１５）に代入すると、
【数９６】
2015026259.tif 000098
【０１５７】
となる。粘性によってｚ＝０で接線方向に作用する力は、
【数９７】
2015026259.tif 000099
【０１５８】
である。
【０１５９】
本文の円筒管の問題では、ｘを周方向とみて、周方向に下記の力が作用する。
【数９８】
2015026259.tif 000100
【０１６０】
虚数単位を陽に含むので
【数９９】
2015026259.tif 000101
【０１６１】
と変形し、第２項目ではｉωを時間微分作用素とみなす。また、ωには近似的に流体付加質量を考慮しない管の固有振動数を与える。
【０１６２】
なお、上述の実施の形態は、本発明の一具体例を示すものであり、本発明を限定するものではない。例えば、この管群振動の予測が適用されるのは加圧水型原子炉の蒸気発生器に限らず、例えば一般の熱交換器やボイラに適用することもできる。また、管は中空なものに限らず、一般の構造物に適用することができ、例えば建物を風が通過する際の振動の予測にも適用することができる。
【符号の説明】
【０１６３】
１０  管群
１１  管
１２  セル
２０  計算装置
【図１】
2015026259.tif 000123
【図２】
2015026259.tif 000124
【図３】
2015026259.tif 000125
【図４】
2015026259.tif 000126
【図５】
2015026259.tif 000127
【図６】
2015026259.tif 000128
【図７】
2015026259.tif 000129
【図８】
2015026259.tif 000130
【図９】
2015026259.tif 000131