6252011 管群振動予測方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6252011

(24)【登録日】2017年12月8日

(45)【発行日】2017年12月27日

(54)【発明の名称】管群振動予測方法

(51)【国際特許分類】

   G06F 17/50 20060101AFI20171218BHJP

   G06F 19/00 20110101ALI20171218BHJP

   F22B 1/16 20060101ALI20171218BHJP

【ＦＩ】

   G06F17/50 612J

   G06F17/50 650Z

   G06F19/00 110

   F22B1/16 A

【請求項の数】10

【全頁数】28

(21)【出願番号】特願2013-155501(P2013-155501)

(22)【出願日】2013年7月26日

(65)【公開番号】特開2015-26259(P2015-26259A)

(43)【公開日】2015年2月5日

【審査請求日】2016年5月26日

(73)【特許権者】

【識別番号】000000099

【氏名又は名称】株式会社ＩＨＩ

(74)【代理人】

【識別番号】100083806

【弁理士】

【氏名又は名称】三好 秀和

(74)【代理人】

【識別番号】100100712

【弁理士】

【氏名又は名称】岩▲崎▼ 幸邦

(74)【代理人】

【識別番号】100101247

【弁理士】

【氏名又は名称】高橋 俊一

(74)【代理人】

【識別番号】100095500

【弁理士】

【氏名又は名称】伊藤 正和

(74)【代理人】

【識別番号】100098327

【弁理士】

【氏名又は名称】高松 俊雄

(72)【発明者】

【氏名】内海 雅彦

【審査官】 合田 幸裕

(56)【参考文献】

【文献】 特開２０１２−１３７１１６（ＪＰ，Ａ）

【文献】 特開平０３−０６３８７９（ＪＰ，Ａ）

【文献】 木村 憲明，直交流による管群の流力弾性振動（ポテンシャル流理論による不安定振動の解析），日本機械学会論文集Ｃ編，日本，社団法人日本機械学会，１９８２年，第48巻／第429号，第662-671頁

【文献】 原 文雄，流体関連ダイナミクス（その４） 管群の流力弾性振動，ターボ機械，日本，一般社団法人ターボ機械協会，１９９３年，第21巻／第4号，第254-258頁

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １７／５０

Ｆ２２Ｂ １／１６

Ｇ０６Ｆ １９／００

ＩＥＥＥ Ｘｐｌｏｒｅ

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

コンピュータが実行する、回転三角配列に配置された所定方向に延びる複数の所定径の管を含む管群の軸方向に直交して流れる流体の流速及び圧力を予測する方法であって、
前記管群の軸方向に直交する面内において、各管をそれぞれ含む複数の領域に分割するステップと、
各領域の外周において流体の流速及び圧力をモード展開した外周のモード座標を用いて外周の少なくとも一部の流体の流速及び圧力を表すステップと、
管の中心を始点として前記領域に含まれる管との境界をなす内周と外周を結ぶ動径について、前記外周のモード座標と、流体の流速及び圧力を周方向にモード展開した周方向のモード座標とを用いて、各領域内の流体の流速及び圧力を表すステップと、を含み、
各領域は、隣接する管の中心を結ぶ直線の少なくとも１つを用いて分割したものであることを特徴とする管群振動予測方法。

【請求項2】

前記周方向のモード座標は、前記面内において管の中心を軸とした極座標

【数1】

について、管の半径をａとして、径方向の

【数2】

及び周方向の

【数3】

において、時刻ｔにおける流速の径方向及び周方向の成分、圧力

【数4】

を各領域の外周における流速の径方向及び周方向の成分、圧力

【数5】

が寄与する第１項

【数6】

と、整数ｅを領域の番号、整数ｍ及びｋを周方向及び径方向のモード次数、

【数7】

として、各領域の内部における流速の径方向及び周方向の成分、圧力を直交基底で展開した第２項

【数8】

との和として表したときに、第２項の展開係数

【数9】

として与えられることを特徴とする請求項１記載の管群振動予測方法。

【請求項3】

各領域は、前記管群を含む領域について各管を含むように分割した正六角形の少なくとも一部を含み、前記外周のモード座標は、正六角形の外周において流体の流速及び圧力をモード展開したものであることを特徴とする請求項１又は２に記載の管群振動予測方法。

【請求項4】

前記外周のモード座標は、前記面内における直交座標

【数10】

について、整数ｉ，ｊを正六角形の辺の頂点又は中点の番号、

【数11】

をこれら頂点又は中点における座標として、時刻ｔにおける辺上での流速のｘ，ｙ成分、圧力

【数12】

を頂点又は中点における流速のｘ，ｙ成分、圧力

【数13】

が辺上の位置

【数14】

に寄与する第１項及び第２項

【数15】

と、整数ｃを辺の番号、整数ｎをモード次数として、辺上における流速のｘ，ｙ成分及び圧力を辺上の直交基底で展開した第３項

【数16】

との和として表したときに、第３項の展開係数

【数17】

として与えられることを特徴とする請求項２に記載の管群振動予測方法。

【請求項5】

時刻ｔにおける流速の径方向及び周方向の成分、圧力

【数18】

の外周における値

【数19】

時刻ｔにおける外周での流速のｘ，ｙ成分、圧力

【数20】

は、方向余弦を用いて

【数21】

なる関係を有することを特徴とする請求項４に記載の管群振動予測方法。

【請求項6】

ナビエ・ストークスの式及び連続条件に基づいて、前記モード座標を一般化座標として、前記流速及び圧力の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記流体の定常流による流れ場を決定することを特徴とする請求項１〜５のいずれか一項に記載の管群振動予測方法。

【請求項7】

流体の流速及び圧力にそれぞれ振動外乱を重畳し、これら流速、圧力及び振動外乱についてナビエ・ストークスの式、連続条件及び前記定常流による流れ場に基づいて、前記振動外乱の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記管群の振動を予測することを特徴とする請求項６に記載の管群振動予測方法。

【請求項8】

前記振動外乱の一般化座標の時間に関する常微分方程式に基づいて、前記管群の不安定振動が発生する前記流体の臨界流速を決定することを特徴とする請求項７に記載の管群振動予測方法。

【請求項9】

前記常微分方程式は、ナビエ・ストークスの式及び連続条件に基づいて、ガレルキン法を用いて求めることを特徴とする請求項６〜８のいずれか一項に記載の管群振動予測方法。

【請求項10】

前記管群は加圧水型原子炉の蒸気発生器における伝熱管群を構成し、前記流体は前記蒸気発生器における二次冷却水であることを特徴とする請求項１〜９のいずれか一項に記載の管群振動予測方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、回転三角配列に配置された管群に直交して流れる流体により発生する管群の振動を予測する管群振動予測方法に関する。

【背景技術】

【0002】

従来、加圧水型原子炉（ＰＷＲ）に備えられる蒸気発生器（ＳＧ）においては、原子炉から供給される一次冷却水を流す伝熱管を並列配置して管群とし、この管群を伝熱面としてこの管群の軸方向と直交する向きに二次冷却水を流して熱交換している。この伝熱管の管群の配置には、回転三角配列、ノーマル三角配列、正方配列等が使用されている。この内、回転三角配列は、高い熱交換の効率が期待できるので広範に使用されている。

【0003】

蒸気発生器においては、熱交換の効率化を図るため、管群に直交した流れの流速を高めることが求められている。しかしながら、この直交流の流速がある値に達すると、不安定振動が発生することが知られている（下記非特許文献１参照）。原子炉の安全性を確保するためにはこのような不安定振動を制御することが必要であり、不安定振動を引き起こす臨界流速を予測する研究がなされてきた。

【0004】

この臨界流速の予測には、ＣＦＤ（数値流体力学）による数値シミュレーション（下記非特許文献２参照）が提供されている。また、１本の管の振動によって隣接する複数の管に作用する力を実験により求め、この実験データに基づいて不安定振動解析モデルを作成する方法が提供されている（下記非特許文献３参照）。

【先行技術文献】

【非特許文献】

【0005】

【非特許文献1】日本機械学会編、事例に学ぶ流体関連振動、ｐ．６６

【非特許文献2】市岡丈彦、川田裕、中村友道、泉元、小林敏雄、高松洋、「流動解析を用いた管群の流力弾性振動に関する研究（第１報、並列２円柱および１列管群の解析）」、日本機械学会論文集（Ｂ編）、１９９５年、６１巻、５８２号、ｐ．１４５−１５１

【非特許文献3】田中博喜、「円柱群の流力弾性振動に関する研究（流水中の円柱群）」、昭和５６年、日本機械学会論文集（Ｂ編）、４７巻、４１９号、ｐ．１２２４−１２３３

【発明の概要】

【発明が解決しようとする課題】

【0006】

しかしながら、ＣＦＤによるシミュレーションは、モデルの入力、計算処理等に多くの時間とコストを要し、多数のパラメータについて臨界流速を予測することは困難であった。また、実験データに基づいて不安定振動解析モデルを作成する方法についても、多数のパラメータについて臨界流速を予測することは困難であった。

【0007】

この出願に係る発明は、上述の実情に鑑みて提案されるものであって、三角配列で配置された管群に対する直交流について、多くの時間とコストを要することなく、多数のパラメータについて臨界流速を予測することができるような管群振動予測方法を提供することを目的とする。

【課題を解決するための手段】

【0008】

上述の課題を解決するために、本発明に係る管群振動予測方法の発明は、コンピュータを用いて、回転三角配列に配置された所定方向に延びる複数の所定径の管を含む管群の軸方向に直交して流れる流体の流速及び圧力を予測する方法であって、前記管群の軸方向に直交する面内において、各管をそれぞれ含む複数の領域に分割するステップと、各領域の外周において流体の流速及び圧力をモード展開した外周のモード座標を用いて外周の少なくとも一部の流体の流速及び圧力を表すステップと、管の中心を始点として前記領域に含まれる管との境界をなす内周と外周を結ぶ動径について、前記外周のモード座標と、流体の流速及び圧力を周方向にモード展開した周方向のモード座標とを用いて、各領域内の流体の流速及び圧力を表すステップと、を含み、各領域は、隣接する管の中心を結ぶ直線の少なくとも１つを用いて分割したものである。

【0009】

前記周方向のモード座標は、前記面内において管の中心を軸とした極座標

【数1】

【0010】

について、管の半径をａとして、径方向の

【数2】

【0011】

及び周方向の

【数3】

【0012】

において、時刻ｔにおける流速の径方向及び周方向の成分、圧力

【数4】

【0013】

を各領域の外周における流速の径方向及び周方向の成分、圧力

【数5】

【0014】

が寄与する第１項

【数6】

【0015】

と、整数ｅを領域の番号、整数ｍ及びｋを周方向及び径方向のモード次数、

【数7】

【0016】

として、各領域の内部における流速の径方向及び周方向の成分、圧力を直交基底で展開した第２項

【数8】

【0017】

との和として表したときに、第２項の展開係数

【数9】

【0018】

として与えられることが好ましい。

【0019】

各領域は、前記管群を含む領域について各管を含むように分割した正六角形の少なくとも一部を含み、前記外周のモード座標は、正六角形の外周において流体の流速及び圧力をモード展開したものであることが好ましい。

【0020】

前記外周のモード座標は、前記面内における直交座標

【数10】

【0021】

について、整数ｉ，ｊを正六角形の辺の頂点又は中点の番号、

【数11】

【0022】

をこれら頂点又は中点における座標として、時刻ｔにおける辺上での流速のｘ，ｙ成分、圧力

【数12】

【0023】

を頂点又は中点における流速のｘ，ｙ成分、圧力

【数13】

【0024】

が辺上の位置

【数14】

【0025】

に寄与する第１項及び第２項

【数15】

【0026】

と、整数ｃを辺の番号、整数ｎをモード次数として、辺上における流速のｘ，ｙ成分及び圧力を辺上の直交基底で展開した第３項

【数16】

【0027】

との和として表したときに、第３項の展開係数

【数17】

【0028】

として与えられることが好ましい。

【0029】

時刻ｔにおける流速の径方向及び周方向の成分、圧力

【数18】

【0030】

の外周における値

【数19】

【0031】

時刻ｔにおける外周での流速のｘ，ｙ成分、圧力

【数20】

【0032】

は、方向余弦を用いて

【数21】

【0033】

なる関係を有することが好ましい。

【0034】

ナビエ・ストークスの式及び連続条件に基づいて、前記モード座標を一般化座標として、前記流速及び圧力の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記流体の定常流による流れ場を決定することが好ましい。

【0035】

流体の流速及び圧力にそれぞれ振動外乱を重畳し、これら流速、圧力及び振動外乱についてナビエ・ストークスの式、連続条件及び前記定常流による流れ場に基づいて、前記振動外乱の一般化座標の時間に関する常微分方程式を導き、この常微分方程式に基づいて前記管群の振動を予測することが好ましい。

【0036】

前記振動外乱の一般化座標の時間に関する常微分方程式に基づいて、前記管群の不安定振動が発生する前記流体の臨界流速を決定することが好ましい。

【0037】

前記常微分方程式は、ナビエ・ストークスの式及び連続条件に基づいて、ガレルキン法を用いて求めることが好ましい。

【0038】

前記管群は加圧水型原子炉の蒸気発生器における伝熱管群を構成し、前記流体は前記蒸気発生器における二次冷却水であることが好ましい。

【発明の効果】

【0039】

この発明は、回転三角配列により配置された管群の正六角形のセルについて、モード展開により解析的に得られた常微分方程式を使用しているので、多くの時間とコストを要することなく、多数のパラメータについて臨界流速を予測することができる。

【図面の簡単な説明】

【0040】

【図1】管群振動予測方法の一連の流れを示す図である。

【図2】回転三角配列の管群を示す断面図である。

【図3】１セルにおける座標設定を示す図である。

【図4】三角配列の数値例題を示す図である。

【図5】換算流速の臨界値の質量減衰パラメータ依存性を示す図である。

【図6】正方配列の管群を示す断面図である。

【図7】管群振動を予測する計算装置の構成を示す図である。

【図8】計算装置における一連の動作の流れを示す図である。

【図9】粘性力の導出を示す図である。

【発明を実施するための形態】

【0041】

以下、本発明に係る管群振動予測方法について、図面を参照して詳細に説明する。

【0042】

〔１．本実施の形態の概要〕
図１は、本実施の形態の管群振動予測方法の一連の流れを示す図である。最初のステップＳ１１では、管群の軸の方向に直交する平面を管の属する正六角形領域に分割する。この領域を以下、セルと称する。ステップＳ１２では、セル外周（各辺）での未知量（流速成分、圧力）を、セルの各頂点での値の１次補間と補間誤差を補うモード展開式の和で表す。ステップＳ１３では、セルの外周と内周（管周囲）を結ぶ任意周方向座標の動径上で、任意の動径座標での未知量を、ステップＳ１２で設定した外周での解に適合する１次式とモード展開式の和で表す。ステップＳ１４では、ステップＳ１２〜Ｓ１３によって得た解の表示式を流体の支配方程式系に代入して、変分直接法（ガレルキン法）で時間の常微分方程式系に離散化し、管群の運動方程式と連成させて解く。

【0043】

ステップＳ１４に続いて、次の双方のステップ（Ａ）、（Ｂ）をさらに実行することもできる。ここで、ステップ（Ａ）はステップ（Ｂ）の前段階として必要である。

【0044】

（Ａ）各管静止条件下での定常（平均）流れ場の決定（流れ解析）
（Ｂ）振動外乱に関する安定性解析
本実施の形態では、セルを正六角形の形状にすることにより、回転三角配列に配置した管群に適用するものである。ここで、回転三角配列とは、三角配列の内、管群の軸方向に直交する平面内において隣接する管を結ぶ直線の一つの方向に流体が流れるものをいう。

【0045】

〔２．解析方法〕
図２のように流れを受ける管群１０について、不安定振動が発生する臨界流速を求める。この管群１０は、それぞれ半径ａの第１の管１１_１から第１１の管１１_１１までの１１本の管１１_１〜１１_１１が、軸間の間隔（ピッチ）をＴとして、各管１１_１〜１１_１１の軸が延びるｚ方向について回転三角配列に配置され、Ｘ方向に沿って３列に並んでいる。以下では簡単のため、各管１１_１〜１１_１１を特定する添字１〜１１を省略する。

【0046】

管１１の軸方向に流速、圧力、管の変位の変化しない２次元問題を考え、この軸に直交する２次元平面を想定する。この面内において、最初に各管１１の静止条件下での定常（平均）流れ場を決定し（２．１節）、次に振動外乱を重畳させて安定性解析を行う（２．２節）。

【0047】

〔２．１ 流れ場の決定〕
〔２．１．１ 解の表示〕
図３において、各管１１を含む正六角形領域である各セル１２に関して座標系

【数22】

【0048】

を定義する。

【0049】

また、セル１２の頂点又は中点の番号ｉを定義し、頂点又は中点ｉでのｘ，ｙ方向の流速、圧力をそれぞれ

【数23】

【0050】

とする。さらに、各外周辺に沿って、座標ｓをとり、各外周辺上でのｘ，ｙ方向の流速、圧力をそれぞれ

【数24】

【0051】

とする。

【0052】

これら各外周辺上での流速、圧力は、辺の頂点又は中点ｉ，ｊでの値の直線補間項に精度向上のための外周辺上の直交基底についてのモード展開項を付加して、次のように表すことができる。ここで、ｃは辺の番号を表している。

【数25】

【0053】

ここで

【数26】

【0054】

である。

【0055】

また、

【数27】

【0056】

である。

【0057】

流体運動解析は、各セル１２に関する極座標を用いて行う。このため、式（１）の流速成分を、次のように極座標に変換する。

【数28】

【0058】

ここで、方向余弦は

【数29】

【0059】

である。

【0060】

周方向座標

【数30】

【0061】

での内周と外周を結ぶ動径上での流速、圧力を、次のように表す。

【数31】

【0062】

ここで

【数32】

【0063】

である。

【0064】

式（５）右辺第１項は、内周ｒ＝ａ、外周ｒ＝ｒ_ｍａｘでそれぞれ０，１となり、外周で解（３）に等しくなるための項である。式（５）右辺第２項は、精度向上のための直交関数展開項で、

【数33】

【0065】

である。このような一般化座標は、モード次数で展開したものであり、モード座標と称することもある。内周での流速は、境界層が非常に薄いため、法線方向成分のみが０となる境界条件を満たすように設定している。

【0066】

本実施の形態の解析法は、
（１）式（１）のｉ，ｊに相当する点をセルの頂点だけでなく辺の中点にも設定でき、周方向角

【数34】

【0067】

の区間を各管ごとに変えられる、
（２）周方向モード関数を流速成分、圧力

【数35】

【0068】

について個別に設定できる、という２つの条件を満たすものである。

【0069】

このような条件（１），（２）を満たすことにより、

【数36】

【0070】

について、例えば

【数37】

【0071】

のように、回転三角配列に適用することを可能にしている。

【0072】

式（３）に式（１）を代入した結果を、式（５）右辺第１項に代入する。これによって、各セル１２内の任意位置での流速と圧力が、頂点又は中点でのｘ，ｙ方向の流速と圧力、式（１）右辺第３項のモード展開係数（モード座標）、式（５）右辺第２項のモード展開係数（モード座標）とで表される。すなわち、流体の無限自由度が、これらの時間変数に縮小する。そこで、このようにして得た式を自由度縮小式と称する。

【0073】

〔２．１．２ ガレルキン法〕
自由度縮小式を、下記の連続条件とナビエ・ストークス方程式の変分重み付き表示式に代入して、未知の時間変数に関する常微分方程式を導く。

【数38】

【数39】

【数40】

【0074】

ここで

【数41】

【0075】

である。

【0076】

常微分方程式は、次のマトリックス方程式の形に導く。

【数42】

【0077】

ここで

【数43】

【0078】

である。

【0079】

式（１０）を時間積分し、時間変動の小さい定常近似解に達した解によって流れ場を決定する。このとき数値計算の安定化を行った。多用されている人工粘性は流れ場を変えてしまうので使わず、運動方程式の時間微分を陽に含む非定常慣性項について係数ρ_ｆを

【数44】

【0080】

とする人工慣性を導入した。粘性項と異なり非定常慣性項は定常状態で０に近づくため、求めるべき定常近似解への影響は小さいことに留意した工夫である。

【0081】

〔２．２ 振動解析〕
〔２．２．１ 解の表示〕
解を次のように設定する。

【数45】

【0082】

これは、前節で求めた定常成分

【数46】

【0083】

に振動外乱

【数47】

【0084】

が重畳した形である。

【0085】

２．１．１節の自由度低減手法を、振動成分に関する方程式系の導出にも用いる。まず、セル１２の外周辺上でのｘ，ｙ方向の流速、圧力に関する振動成分を、頂点又は中点ｉ，ｊでの値

【数48】

【0086】

の直線補間と、モード展開項の和で表し、式（１）に対応して次式を得る。

【数49】

【0087】

流体運動解析は、図２のような、各セル１２に関する極座標を用いて行うため、流速成分を極座標に変換し、周方向座標

【数50】

【0088】

での内周と外周を結ぶ動径上での流速、圧力（振動成分）を、式（６）に対応して次のように表す。

【数51】

【0089】

ここで

【数52】

【0090】

である。

【0091】

このようにして、各セル１２内の任意位置での流速と圧力（振動成分）を、頂点又は中点でのｘ，ｙ方向の流速と圧力、式（１２）右辺第３項のモード展開係数、式（１４）右辺第２項のモード展開係数で表す。すなわち、流体の無限自由度を、これらの時間変数に縮小する。

【0092】

〔２．２．２ ガレルキン法〕
ナビエ・ストークス方程式は、変分重み付きの形で次のように表される。

【数53】

【数54】

【0093】

運動学的条件は、管１１表面のｘ，ｙ方向の速度（ｄ／ｄｔ）ｑ_ｘ，（ｄ／ｄｔ）ｑ_ｙの法線方向成分を管１１表面での湧き出しと考え、湧き出し項を流体連続条件式に導入することによって、次のようになる。

【数55】

【0094】

各管の運動方程式は、次式によって与えられる。

【数56】

【0095】

ここで

【数57】

【0096】

である。

【0097】

式（１４）を式（１５），（１６），（１７），（１８）に代入することにより、未知の時間変数に関する常微分方程式系を下記のマトリックス方程式の形に導く。

【数58】

【0098】

ここで

【数59】

【0099】

は一般化座標によって構成される列ベクトルであり、マトリックス

【数60】

【0100】

の固有値を求めることによって安定性を判定する。

【0101】

〔３．計算結果〕
〔３．１ 計算条件〕
次の緒元について計算を行った。

【数61】

【0102】

流れ場の決定では、粘性係数の代わりに乱流粘性係数をμの６０倍の値に設定した（ＣＦＤにおいては２０〜１００倍の値を用いることに基づく）。

【0103】

管群１０の問題では、管１１が多数あるため、図４の太い実線の直線で仕切られ、境界Ａ，Ｂにそれぞれ正、負の、大きさ

【数62】

【0104】

の湧き出しを与えた流体領域を縮小モデルとした。正および負の湧き出しの影響の少ない管Ｄの

【数63】

【0105】

での流れ場を図４の１１本の管の

【数64】

【0106】

の領域に設定し、これらの各管の

【数65】

【0107】

での流れ場は、

【数66】

【0108】

を結ぶ直線に関する対称性より決定した。

【0109】

〔３．２ 臨界流速の従来データとの比較〕
質量減衰パラメータと呼ばれる次の量

【数67】

【0110】

を、減衰比ζ_bを変えることによって変え、換算流速（reduced velocity）と称されている下記の量の臨界値を求めた。

【数68】

【0111】

ただし、

【数69】

【0112】

である。

【0113】

管１１の間の隙間の流速Ｖ_ｃは、流れ場を決めるときの式（９）中の湧き出し、吸込みソースの強さＶ_ｓｒｃのＴ／（Ｔ−２ａ）倍である。臨界値を従来データ（非特許文献１）と比較した結果を図５に示す。図中の“△”は従来データであり、図中の曲線Ａは本解析により得られた値である。本解析により従来データに近い値が得られていることが確かめられる。

【0114】

上記は、流体密度の大きい液体の場合を対象とした、質量減衰パラメータの小さい範囲の結果である。しかし工学上は、蒸気発生器（ＳＧ）のように、ボイド率９０％、すなわち流体密度が１／１０に低下し、式（２０）の質量減衰パラメータがより大きい場合が重要になる。この場合の結果を図中の曲線Ｂで同様に示した。この場合、

【数70】

【0115】

である。

【0116】

〔３．３自由度低減効果〕
本解析では、自由度数がＣＦＤに比べると少なくて済むことを確かめる。解の表示式（１２），（１３）より、流体解析に関する自由度数Ｎは次式によって求められる。

【数71】

【0117】

ここで、Ｎ_１は頂点の個数、Ｎ_２は辺の個数、Ｎ_１３は管の個数（セルの個数）であり、係数３は、ｘ，ｙ方向の流速成分と圧力とで３個の変数があることによる。図４で計算に用いた配列の場合について自由度数を求めると、Ｎ_１＝３６，Ｎ_２＝３６，ｎ_ｍａｘ＝２，Ｎ_３＝１１であり、ｍ_ｍａｘ＝１，ｋ_ｍａｘ＝１で上のように従来データを説明し得る解析結果が得られる。

【0118】

この場合の自由度数はＮ＝３９０である。ＣＦＤの数百万オーダに比べるとかなり少なくて済む。安定判別の固有値を求めるために要する計算時間は１ケースについて２分程度である（ただし臨界流速を求めるためには、固有値計算を１０回程度繰り返す必要がある）。このような低次元の固有値問題に帰着できることは、ＣＦＤ使用の際に流体と管１１の周波数差により流体解析に合わせて時間きざみを非常に小さくして時間積分する必要が生じる問題の回避に有効である。

【0119】

〔４．比較例〕
本実施の形態は回転三角配列に配置された管群に適用されるが、例えば正方配列により配置された管群に適用される場合には、対称性を考慮してモード座標の設定や流れの計算を簡略化することができる。ここでは、回転三角配列の場合との比較のため、正方配列により配置された管群への適用について簡単に説明する。

【0120】

図６は、正方配列により配置された管群１０を示す図である。この図においては、各管１１が正方配列により配置されている他は回転三角配列により配置された図２と同様の構成であるので、対応する符号を付して説明を省略する。

【0121】

このような正方配列の場合、内周と外周を結ぶ動径上での流速、圧力の式（５）は、次のように簡略化することができる。

【数72】

【0122】

ここで

【数73】

【0123】

である。

【0124】

そして、流れ場の計算は、図６において、３個の領域

【数74】

【0125】

での流れ場は同じである。そこで、１個の領域のみについて、流れ場を決定する計算を行うことができる。次に、Ｘ方向についても管１１が多数であるため周期対称な問題である。そこで、出入り口での湧き出し、吸込みの影響の少ない領域

【数75】

【0126】

で求まった流れ場を９本の管１１の周りのセル１２に設定して安定性解析を行うことができる。

【0127】

本実施の形態の回転三角配列の場合は、このような対称性を利用することができないので、式（５）のように周方向の角度に対する依存性を含む流速、圧力を導出する必要がある。

【0128】

〔５．計算装置〕
前述したような臨界流速の予測は、図７に示すような計算装置２０によって実現することができる。この計算装置２０は、ＣＰＵ、ＤＳＰの如き演算部２１、ＲＡＭ、ＲＯＭ、ハードディスクの如き記憶部２２、ＬＣＤ、プリンタの如き出力部２３、キーボード、マウスの如き入力部２４を含み、例えばパーソナルコンピュータを利用することができる。

【0129】

図８に示す計算装置２０の一連の動作は、記憶部２２に格納された臨界流速算定プログラム２２ａを演算部２１が読み出して実行することにより実現される。この臨界流速算定プログラム２２ａは、前述のような手順によって得られた流速及び圧力の一般化座標の時間に関する常微分方程式と振動外乱の一般化座標の時間に関する常微分方程式の表現を含んでいる。

【0130】

最初のステップＳ２１においては、モデルを設定する。ここでは、前述したような管群１０のモデルを設定するものとする。入力部２４は、このモデルについて、管１１の径ａ、管１１を配置した間隔Ｔ、管１１の単位長さ当たりの質量、管１１の固有振動数、流体の密度、流体の粘性係数の少なくとも１つの値を入力値として受け取る。演算部２１は、入力部２４が受け取った入力値を記憶部２２に格納する。

【0131】

ステップＳ２２においては、演算部２１は、記憶部２２に格納された入力値を読み出し、その数値計算部２１ａにおいて、この入力値に基づいて前記常微分方程式の表現を用いてこの管群における臨界流速について数値計算する。演算部２１は、得られた臨界流速の値を記憶部２２に格納する。

【0132】

ステップＳ２３においては、演算部２１は、記憶部２２に格納された臨界流速の値と、同じく記憶部２２に格納された所定の閾値２２ｂとを読み出す。演算部２１は、その判定部２１ｂにおいて、臨界流速値が閾値２２ｂを超えた場合にはＯＫと判定して一連のステップを終了する。一方、臨界流速が閾値を越えない場合にはＮＧとして判定して前のステップＳ２１のモデル設定に手順を戻す。なお、閾値２２ｂは、入力部２４を介して設定することができる。

【0133】

このような一連の工程において、臨界流速の値が所定の閾値を超えるまでモデル設定、数値計算、判定のループを繰り返すことにより、臨界流速が閾値を超えるようなモデル設定を可能としている。また、前述のような常微分方程式の表現を利用することにより、精度の高いモデル設定を可能としている。なお、これに限られず、臨界流速が閾値内に収まるようにモデル設定することもできる。

【0134】

なお、このような臨界流速の算定は、記憶部２２に格納した臨界流速算定プログラム２２ａのような、前述の常微分方程式の表現を含み、モデル設定、数値計算、判定のステップを有するプログラムによっても提供することができる。

【0135】

〔６．付録：式（１８）における粘性力の導出〕
図９に示すように液体とｚ＝０の平面３１で接する物体３２のｘ方向の速度が次式

【数76】

【0136】

のようにある周波数ωで振動する場合の境界層流れを考える。ｘ方向の流速をｖ_ｘ、圧力をｐとし、これらの境界層のすぐ外側の主流での値を、上添字（ｍ）を付けて表わす。

【0137】

主流に関する運動方程式は

【数77】

【0138】

となる。

【0139】

この運動方程式の解を

【数78】

【0140】

とおき、式（Ａ２）に代入すると、次式が得られる。

【数79】

【0141】

境界層内流れを解析するため、ｚ方向の流速は非常に小さいと仮定すると、ｚ方向のナビエ・ストークス方程式は

【数80】

【0142】

に帰着する。

【0143】

従って、条件

【数81】

【0144】

によって次式が要求される。

【数82】

【0145】

境界層は非常に薄いので、近似

【数83】

【0146】

をｘ方向のナビエ・ストークス方程式に適用できる。このようにして次式が得られる。

【数84】

【0147】

式（Ａ７）を下記の境界条件の下で解く。

【数85】

【数86】

【0148】

解を

【数87】

【0149】

と表わし、式（Ａ１０）と式（Ａ３）の第２式を式（Ａ７）〜（Ａ９）に代入すると、次のようになる。

【数88】

【数89】

【数90】

【0150】

式（Ａ４）を式（Ａ１１）に代入して次式を得る。

【数91】

【0151】

式（Ａ１４）に関する同次方程式の特性根αは、

【数92】

【0152】

によって定まり、解が有限である条件より、

【数93】

【0153】

となる。

【0154】

従って、式（Ａ１４）の一般解は次のようになる。

【数94】

【0155】

Ｃ_ｘ（ｘ）は任意関数で、式（Ａ１３）の条件を課すことにより、次のように物体速度と主流速度で表わせる。

【数95】

【0156】

式（Ａ１６）を式（Ａ１５）に代入すると、

【数96】

【0157】

となる。粘性によってｚ＝０で接線方向に作用する力は、

【数97】

【0158】

である。

【0159】

本文の円筒管の問題では、ｘを周方向とみて、周方向に下記の力が作用する。

【数98】

【0160】

虚数単位を陽に含むので

【数99】

【0161】

と変形し、第２項目ではｉωを時間微分作用素とみなす。また、ωには近似的に流体付加質量を考慮しない管の固有振動数を与える。

【0162】

なお、上述の実施の形態は、本発明の一具体例を示すものであり、本発明を限定するものではない。例えば、この管群振動の予測が適用されるのは加圧水型原子炉の蒸気発生器に限らず、例えば一般の熱交換器やボイラに適用することもできる。また、管は中空なものに限らず、一般の構造物に適用することができ、例えば建物を風が通過する際の振動の予測にも適用することができる。

【符号の説明】

【0163】

１０ 管群
１１ 管
１２ セル
２０ 計算装置

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】