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