流体構造連成解析方法及び流体構造連成解析プログラム
【課題】 キャビテーションを考慮して、実際の挙動により近いシミュレーションを行うことのできる流体構造連成解析方法及び流体構造連成解析プログラムを提供する。
【解決手段】 コンピュータ２０の固体計算手段２１は、運動量保存式を計算して構造物の変位量を算出する。コンピュータ２０の流体計算手段２２は、構造物の変位量に基づいて境界流速とヤコビアンを計算し、圧力をコンピュータ２０の気泡成長計算手段２３に引渡す。気泡成長計算手段２３は、圧力に基づいて気泡半径についての時間発展方程式を計算して、気泡半径、気泡成長速度及びボイド率を流体計算手段２２に引渡す。流体計算手段２２は、ヤコビアン、気泡半径、気泡成長速度及びボイド率を用いて、気泡の質量の項を組み込んだ質量保存式を用いて圧力を算出し、これを固体計算手段２１に引渡す。コンピュータ２０は、流体計算手段から引渡された圧力を用いて変位量を算出する。
【解決手段】 コンピュータ２０の固体計算手段２１は、運動量保存式を計算して構造物の変位量を算出する。コンピュータ２０の流体計算手段２２は、構造物の変位量に基づいて境界流速とヤコビアンを計算し、圧力をコンピュータ２０の気泡成長計算手段２３に引渡す。気泡成長計算手段２３は、圧力に基づいて気泡半径についての時間発展方程式を計算して、気泡半径、気泡成長速度及びボイド率を流体計算手段２２に引渡す。流体計算手段２２は、ヤコビアン、気泡半径、気泡成長速度及びボイド率を用いて、気泡の質量の項を組み込んだ質量保存式を用いて圧力を算出し、これを固体計算手段２１に引渡す。コンピュータ２０は、流体計算手段から引渡された圧力を用いて変位量を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、構造物の中の流体の挙動を解析するための流体構造連成解析方法及び流体構造連成解析プログラムに関する。
【背景技術】
【０００２】
構造物の中の流体の挙動を解析するために、各種シミュレーション技術が検討されている。例えば、核破砕タイプ中性子源の開発においても、計算機シミュレーションの検討が進められている。この開発では、加速器で加速された陽子をターゲットとなる物質に照射し、ターゲットに核破砕を起こすことで中性子を発生させる。ここでは、ターゲットには、ステンレス鋼の容器（構造物）と、その中を循環する液体水銀（流体）から構成された水銀ターゲットを用いる。そして、この水銀ターゲットの中の現象を予め予測するため計算機シミュレーションが大切である。この場合、構造物及び流体が相互に影響し合っているので、構造物及び流体に関する方程式を連立させた解析が必要である。このような構造物及び流体に関するシミュレーションを行なうための解析装置が開示されている（例えば、特許文献１参考。）。
【０００３】
特許文献１においては、流体解析による流体の解及び構造解析による構造全体の解のうち、流体と構造が接触する境界における流体の解と構造の解とが一致する共通解を求める。そして、この求められた共通解を用いて、再度、流体解析と構造解析を実行することにより安定解を算出している。
【特許文献１】特開２０００−３５２５４５号公報（第１頁）
【発明の開示】
【発明が解決しようとする課題】
【０００４】
ところで、上述の例のような水銀ターゲットの場合、内部で核破砕という高密度の熱の発生を伴う現象が生じるため、水銀ターゲットの容器には大きな力がかかる。この場合、高密度の熱の発生により液体水銀が急激に膨張し、その結果、生じる衝撃波によりターゲット外壁に大きな応力が発生する。そして、外壁に歪が生じ、水銀中にキャビテーションが発生する可能性がある。このキャビテーションは、外壁の劣化等の影響を与えるため、キャビテーションが発生した場合の外壁への影響の計算することが大切である。
【０００５】
しかし、このような状況下でのキャビテーションを考慮した流体構造連成解析プログラムはない。
本発明は、上述の課題に鑑みてなされ、その目的は、キャビテーションを考慮して、実際の挙動により近いシミュレーションを行うための流体構造連成解析方法及び流体構造連成解析プログラムを提供することにある。
【課題を解決するための手段】
【０００６】
上記問題点を解決するために、請求項１に記載の発明は、コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析方法であって、前記コンピュータが、運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出段階と、流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出段階と、前記圧力算出段階で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定段階と、前記誤差が所定範囲になった場合、前記算出した圧力を用いて
前記構造物の変位を算出する変位算出段階とを含むことを要旨とする。
【０００７】
請求項２に記載の発明は、請求項１に記載の流体構造連成解析方法において、前記コンピュータは、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算段階を更に含むことを要旨とする。
【０００８】
請求項３に記載の発明は、請求項２に記載の流体構造連成解析方法において、前記コンピュータは、有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断段階を更に含み、気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うことを要旨とする。
【０００９】
請求項４に記載の発明は、請求項２又は３に記載の流体構造連成解析方法において、前記コンピュータは、前記圧力算出段階において収束判定に至るまでの繰り返し回数を計数し、前記繰り返し回数が基準回数になっても、前記収束判定段階において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うことを要旨とする。
【００１０】
請求項５に記載の発明は、コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析プログラムであって、前記コンピュータを、運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出手段と、流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出手段、前記圧力算出手段で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定手段、及び前記誤差が所定範囲になった場合、前記算出した圧力を用いて前記構造物の変位を算出する変位算出手段として機能させることを要旨とする。
【００１１】
請求項６に記載の発明は、請求項５に記載の流体構造連成解析プログラムにおいて、前記コンピュータを、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算手段として更に機能させることを要旨とする。
【００１２】
請求項７に記載の発明は、請求項６に記載の流体構造連成解析プログラムにおいて、前記コンピュータを、有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断手段として更に機能させ、気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うようにさせたことを要旨とする。
【００１３】
請求項８に記載の発明は、請求項６又は７に記載の流体構造連成解析プログラムにおいて、前記コンピュータを、前記圧力算出手段において収束判定に至るまでの繰り返し回数を計数し、前記繰り返し回数が基準回数になっても、前記収束判定手段において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うようにさせたことを要旨とする。
【００１４】
（作用）
請求項１又は５に記載の発明によれば、コンピュータは、流体の気泡成長計算から算出
した気泡半径、気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式から、流体の圧力を算出する。このため、急激な圧力変化によって流体に発生するキャビテーションによる気泡を質量保存式に入れて、シミュレーションを行うことができる。従って、キャビテーションを考慮した流体についての計算と、この流体を内包した構造物についての計算とを連成させたシミュレーションを、実際の現象により近づけて行うことができる。
【００１５】
請求項２又は６に記載の発明によれば、コンピュータは、流体の質量保存式及び運動量保存式から圧力が算出された場合には、その圧力を用いて気泡成長計算を行う。このため、圧力変動に応じた気泡の成長を反映したシミュレーションを行うことができる。
【００１６】
請求項３又は７に記載の発明によれば、コンピュータは、有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断し、大きいメッシュについてのみ気泡成長計算を行う。従って、気泡の成長の影響が少ないメッシュについては、気泡成長計算を行わないので、計算量を少なくすることができる。また、この場合であっても、気泡の成長が大きくないメッシュについて、気泡の質量に関する項を含む質量保存式と運動量保存式とから圧力を算出するので、実際の現象に近いシミュレーションを行うことができる。
【００１７】
請求項４又は８に記載の発明によれば、コンピュータは、密度が所定範囲の誤差にならなかった場合には、算出された圧力を用いて再計算を行う。また、コンピュータは、再計算を所定回数行っても、所定範囲の誤差にならない場合には、時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて圧力を再計算する。このため、気泡について、より正確なシミュレーションができるので、実際の現象に近いシミュレーションを行うことができる。
【発明の効果】
【００１８】
本発明によれば、キャビテーションを考慮して、構造物及び流体の挙動に関するシミュレーションを行うことができる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明を具体化した一実施形態を図１〜図９に基づいて説明する。本実施形態においては、水銀ターゲットにおける圧力波伝播のシミュレーションに適用する。ここでは、水銀ターゲットの容器を構造物、その中の水銀を流体の挙動としたときの式を連立して解析を行う。
【００２０】
このシミュレーションを行うコンピュータ２０は、図１に示すように、シミュレーションの結果を出力するための出力手段３０、各指示を入力するための入力手段４０や初期値や設定値に関するデータを記憶するデータ記憶手段５０に接続されている。
【００２１】
出力手段３０は、シミュレーション結果を出力するために用いられ、具体的には、ディスプレイやデータ記憶部などを用いる。入力手段４０は、シミュレーションを行うために各種指示を入力するために用いられ、具体的には、キーボードやマウスなどである。データ記憶手段５０は、シミュレーションを行うために必要なデータを予め記憶させるために用いられる。このデータ記憶手段５０には、例えば、解析を行う固体及び流体の形状やマトリクスの形状、解析に用いられる物性値などに関するデータが記憶されている。
【００２２】
コンピュータ２０は、固体計算手段２１、流体計算手段２２及び気泡成長計算手段２３を備える。コンピュータ２０は、これらを用いて境界流速算出段階、圧力算出段階、収束判定段階、変位算出段階、気泡成長計算段階及び判断段階を行う。コンピュータ２０は、
流体構造連成解析プログラムを実行することにより、境界流速算出手段、圧力算出手段、収束判定手段、変位算出手段、気泡成長計算手段及び判断手段として機能する。ここで、各手段の機能を図２〜図４を用いて説明する。
【００２３】
図２に示すように、固体計算手段２１は、構造物（固体）に関して、式（１−１）、式（１−２）及び式（１−３）の計算処理を実行し、構造物の変位量ｕを算出する。
式（１−１）は、力学的境界条件式であり、構造物が、その周囲の流体から受けた圧力ｐを節点荷重ｔw に変換するための計算式である。なお、この式（１−１）において、ｖaは定数である。
【００２４】
式（１−２）は、応力テンソルとひずみテンソルの構造物についての関係式であり、節点荷重ｔw を応力テンソルσに変換するための計算式である。この式（１−２）において、Ｖv は法線ベクトルである。
【００２５】
式（１−３）は、構造物に関する運動量保存式であり、式（１−２）で算出した応力テンソルσを用いて、構造物の変位量ｕを算出するための計算式である。この式（１−３）において、ρsは固体の密度、∇は微分演算子、ｂfは外力である。なお、本実施形態では、外力ｂfを「０」と設定する。
【００２６】
図３に示すように、流体計算手段２２は、流体に関して、式（２−１）〜式（２−８）の計算処理を実行して、流体の圧力ｐを算出する。
式（２−１）は、固体と流体の境界面（カップリングサーフェース）の座標値Ｘを算出するための計算式である。この式（２−１）において、Ｘnew は新たに算出する境界面の座標値であり、Ｘold は座標値Ｘnew を算出する直前の境界面の座標値である。
【００２７】
式（２−２）は、固体と流体の境界面の境界流速ｖw を算出するための計算式である。ここで、ｄｔはクーラン条件から決まる時間間隔である。
式（２−３）は、直交座標系（ｘ1、ｘ2、ｘ3、ｔ）を一般座標系（ξ1、ξ2、ξ3、Ｔ）に変換するヤコビアンＪT を算出するための計算式である。
【００２８】
式（２−４）は、流体における質量保存式である。この式（２−４）においては、後述するバブルモデルを組み込んでいる。この式（２−４）において、ρは流体の密度、Ｕi
は流速ベクトル（の反変成分）、πは円周率、ｍは単位体積あたりの気泡数、Ｒは気泡
半径、Ｒ'（＝dR/dt）は気泡の半径の時間微分値であり、いわゆる気泡成長速度である。αは、流体におけるボイド率であり、後述する気泡成長計算手段２３によって算出される。
【００２９】
式（２−５）は、流体の運動量保存式である。この式（２−５）において、ｖは速度、∇は微分演算子、Ｆは運動量保存式に現れるその他の項を示している。この式（２−５）及び式（２−４）の２式は、連立することにより流体に関する圧力に関するポアソン方程式となり、これらから流体計算手段２２は圧力ｐを算出する。
【００３０】
式（２−６）は、流体における状態方程式である。この式（２−６）は、本実施形態ではＴａｉｔの方程式であり、圧力ｐに基づいて単位質量あたりの体積Ｖを算出する。この式（２−６）においてＴtは流体の温度であり、初期条件として与えられる値である。ま
た、式（２−６）において、Ｖ（Ｔt、０）＝Ａ・Ｔt ＋Ｂ、Ｂ（Ｔt）＝Ｂ０・Ｔt＋Ｂ1、Ｃ（Ｔt）＝Ｃ０・Ｔt＋Ｃ1である。ここで、Ａ、Ｂ、Ｃはパラメータであり、Ａ、Ｂ
、Ｃ、Ｂ1、Ｃ1の値は、初期値として与えられる。
【００３１】
式（２−７）は、式（２−６）から算出された単位質量あたりの体積Ｖから、状態方程
式に基づく密度ρeos を算出するための式である。
式（２−８）は、収束判定を行うためのパラメータの算出式である。この式（２−８）の値が収束判断値εより小さくなった場合に、流体計算手段２２は、この式（２−４）及び式（２−５）から算出した値を圧力ｐとして決定する。
【００３２】
図４に示すように、気泡成長計算手段２３は、式（３−１）〜式（３−５）の計算処理を実行して、流体に含まれる気泡の成長を計算する。
式（３−１）は、気泡に関する圧力の関係式である。この式（３−１）において、ｐt は気泡内部圧力である。ｐt'（＝ｄｐt／ｄｔb）は内部圧力変化であり、気泡内部圧力ｐtの時間の１階微分である。γは比熱比、Ｋｈは熱伝達係数、Ｔtbは気泡内温度である。
τは、熱伝達係数Ｋｈを温度θで気泡内温度Ｔtbまで積分した値であり、便宜上、ここでは熱伝達積分係数と称する。
【００３３】
式（３−２）は、気泡に関する温度の式である。この式（３−２）と式（３−１）とを連立して解くことにより、気泡内部圧力ｐtを算出する。この式（３−２）において、Ｄ
は温度拡散係数であり、ｙは気泡半径Ｒによって正規化された気泡の外壁から中心点までの距離である。
【００３４】
式（３−３）は、気泡壁の外側の圧力ｐBを算出するための式である。この式（３−３
）には、上記式（３−１）及び式（３−２）を連立して算出した気泡内部圧力ｐtが代入
される。この式（３−３）において、σstは表面張力、μは流体の粘性係数である。
【００３５】
式（３−４）は、気泡半径Ｒ及び気泡成長速度Ｒ'を算出するための時間発展方程式で
ある。この式（３−４）において、Ｒ”（＝d2Ｒ/dt2 ）は気泡半径の時間の２階微分、
ｃは流体の音速である。ｐ∞は流体の静水圧であり、定数パラメータとして設定される値である。
式（３−５）は、ボイド率αを算出するための式である。
【００３６】
従って、コンピュータ２０は、図１に示すように、固体計算手段２１は流体計算手段２２に、算出した構造物の変位量ｕを引渡し、流体計算手段２２は固体計算手段２１に、流体計算手段２２で算出した圧力ｐを引渡す。また、流体計算手段２２は気泡成長計算手段２３に圧力ｐを引渡し、気泡成長計算手段２３は流体計算手段２２に、気泡成長計算手段２３で算出した気泡半径Ｒ、気泡成長速度Ｒ’、ボイド率αを引渡す。
【００３７】
次に、コンピュータ２０の計算処理について、図５〜図９を用いて説明する。ここでは、各記号に付された上付き文字「ｎ」は、タイムステップｎにおける値であることを示している。このタイムステップｎとは、シミュレーションにより構造物及び流体の状態を求める時間をクーラン条件から決まる時間間隔ｄｔで分割したときの時間回数である。また、後述する流体計算において、各記号に付された上付き文字「ｋ」は、１つのタイムステップｎにおける処理ループの回数を示す。
【００３８】
まず、計算処理において、初期条件として、節点荷重ｔw 、有限要素のメッシュ分割の個数、固体の材料の物性、流体の物性、設定温度等に関するデータ、終了タイムステップ数などの値が入力手段４０を介して入力されて、データ記憶手段５０に記憶される。本実施形態では、次の初期条件が設定される。まず、初期の壁面は大気圧と内部の流体の圧力がつりあっていると考えて、節点荷重ｔwに「０」の値が設定される。また、節点荷重ｔw
以外の外力はないと仮定して外力ｂf＝０の値が設定される。更に、通常は各メッシュ中には複数の気泡が含まれているが、本実施形態においては、各メッシュ中に含まれる気泡群を平均化して、１つの気泡の成長起動として代表させるため、単位体積あたりの気泡数ｍを指定する。また、流体計算における１つのタイムステップに対して、気泡成長計算手
段２３が気泡成長の計算を行うときの気泡計算タイムステップｎb数（以下、必要ループ
数と称する）が設定される。更に、水銀ターゲットのシミュレーションにおいては、各メッシュの流体に作用する圧力が初期条件として与えられる。また、本実施形態では内壁の変位量ｕを考慮するために、収束判断値εは１．０×１０-7に設定している。
【００３９】
（固体計算）
まず、コンピュータ２０が固体のシミュレーション計算処理を実行する場合を説明する。具体的には、コンピュータ２０の固体計算手段２１が、まず、力学的境界条件の設定を行う（ステップＳ１−１）。タイムステップｎ＝１の場合、処理の開始の指示を受けた固体計算手段２１は、力学的境界条件として、初期条件として入力された節点荷重ｔwn=1
を用いる。
【００４０】
また、タイムステップｎ≧２のときには、固体計算手段２１は、後述するように流体計算手段２２において算出された流体の圧力ｐの引渡しを受ける。この場合、ステップＳ１−１における力学的境界条件として、後述する流体計算から算出される流体の圧力ｐを式（１−１）に代入して、節点荷重ｔwn を求める。
【００４１】
そして、設定された節点荷重ｔwnを式（１−２）に代入して、固体のカップリングサーフェース上の応力テンソルσを求める。
次に、固体計算手段２１は、算出した応力テンソルσnから式（１−３）を用いて、固
体におけるカップリングサーフェースの変位量ｕnを算出する（ステップＳ１−２）。ま
た、式（１−３）の固体の密度ρsとして、初期条件としてデータ記憶手段５０に記憶さ
れている固体の材料の物性値を、固体計算手段２１は取得して用いる。
【００４２】
そして、固体計算手段２１は、タイムステップが終了タイムステップ数になったか否かを判断する（ステップＳ１−３）。具体的には、固体計算手段２１は、タイムステップ数をカウントし、予め定められた所定数に到達したかどうかを判断する。終了タイムステップ数が終了していない場合（ステップＳ１−３において「ＮＯ」の場合）、固体計算手段２１は、構造物（内壁）の変位量ｕn を流体計算手段２２に引渡す（ステップＳ１−４）。
【００４３】
（流体計算の第１プロセス）
次に、流体計算について説明する。ここでは、流体計算は、気泡成長計算を入れ子構造にしている。そこで、気泡成長計算の前のプロセスを「第１プロセス」、気泡成長計算の後のプロセスを「第２プロセス」と呼ぶ。
【００４４】
固体計算手段２１から構造物の変位量ｕnを引渡された流体計算手段２２は、境界流速
とヤコビアンの計算を行う（ステップＳ２−１）。具体的には、流体計算手段２２は、固体計算手段２１において算出された構造物の変位量ｕnを、式（２−１）及び式（２−２
）に代入し、構造物である内壁と流体である水銀の境界面（カップリングサーフェース）の座標値Ｘn 及び境界流速ｖwnを算出する。ここで、タイムステップｎ＝１の場合、壁を変位量ｕn=1を計算していないため、座標値Ｘold は、初期条件による値、すなわち流体
の各メッシュ要素の初期座標値を用いる。
そして、流体計算手段２２は、算出した流体の座標値Ｘn 及び境界流速ｖwnに基づいて、ヤコビアンＪTを更新する。
【００４５】
次に、流体計算手段２２は、気泡成長計算を起動する（ステップＳ２−２）。具体的には、圧力ｐnの値を気泡成長計算手段２３に引渡す。ここで、タイムステップｎ＝１のと
きには、流体計算手段２２は、初期条件で設定された圧力ｐn=1を気泡成長計算手段２３
に引渡す。また、タイムステップｎ≧２のときには、後述するステップにおいて流体計算
手段２２が算出した圧力ｐ、すなわちタイムステップ「ｎ−１」の圧力ｐn-1の値を気泡
成長計算手段２３に引渡す。
【００４６】
（気泡成長計算）
気泡成長計算においては、コンピュータ２０の気泡成長計算手段２３は、気泡の周囲の流体の圧力ｐnの更新を行う（ステップＳ３−１）。具体的には、初期条件としては、気
泡の周囲の圧力ｐn=1には、設定された流体計算における初期圧力を用いる。また、初期
条件として気泡内部圧力ｐtには、静的平衡状態にあると仮定して、気泡壁の外側の圧力
ｐB＝ｐ＋ｐ∞とした式（３−３）を用いて算出する。また、タイムステップｎ≧２のと
きには、タイムステップｎ−１で算出した流体の圧力ｐn-1を用いる。
【００４７】
次に、気泡成長計算手段２３は、気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを算出する。これは、ステップＳ３−２を予め設定された必要ループ回、繰り返し行うことにより算出する。
【００４８】
次に、ステップＳ３−２の繰り返し処理について図６を用いて詳述する。ここでは、各記号に付された上付き文字「ｊ」は、１つの気泡計算タイムステップｎbにおける処理ル
ープの回数を示す。なお、この気泡計算タイムステップｎbは、１つのタイムステップｎ
において必要ループ数、行われる。
【００４９】
そこで、気泡成長計算手段２３は、タイムステップｎの気泡半径Ｒjについての時間発
展方程式の計算を行う（ステップＳ３−２）。具体的には、気泡成長計算手段２３は、まず、式（３−１）及び式（３−２）を連立させて解くことにより、気泡内部圧力ｐtjと、熱伝達積分係数τjが求まり、この熱伝達積分係数τjから気泡内部温度Ｔtbjが求まる。
この算出した気泡内部圧力ｐtjを式（３−３）に代入して、気泡壁の外側の圧力ｐBjを算出する。この算出した気泡壁の外側の圧力ｐBjを用いて式（３−４）を解く。これにより、この式（３−１）の解として、気泡半径Ｒj及び気泡成長速度Ｒ'j（＝ｄＲj/ｄｔb ）
が得られる。
【００５０】
次に、気泡成長計算手段２３は、予め設定された必要ループを終了したか否かを判断する（ステップＳ３−３）。そして、気泡成長計算手段２３は、予め設定された必要ループが終了していない場合（ステップＳ３−３において「ＮＯ」の場合）には、算出した気泡半径Ｒjと、気泡成長速度Ｒ'jとを用いて、ステップＳ３−２の再計算を行う。
【００５１】
一方、気泡成長計算手段２３は、必要ループが終了した（ステップＳ３−３において「ＹＥＳ」）と判断すると、ボイド率αjの設定を行う（ステップＳ３−４）。具体的には
、気泡成長計算手段２３は、算出した気泡半径Ｒjを式（３−５）に代入することにより
、ボイド率αjを算出する。
【００５２】
次に、気泡成長計算手段２３は、算出した気泡半径Ｒj、気泡成長速度Ｒ'jを気泡半径
Ｒn、気泡成長速度Ｒ'nとして、またボイド率αnを流体計算手段２２に引渡す（ステップＳ３−５）。すなわち、気泡成長計算手段２３は、必要ループに応じてステップＳ３−２を繰り返し行うことにより、より正確な気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを算出する。
【００５３】
（流体計算の第２プロセス）
気泡成長計算手段２３から気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを引渡された流体計算手段２２は、運動量保存式、ポアソン方程式、流速補正及び質量保存式の計算を行う（ステップＳ２−３）。この計算は、３次元座標で表現される全メッシュについて行われる。この計算について図７〜図９を用いて詳述する。
【００５４】
図７に示すように、流体計算手段２２は、まず、陽解法を用いて運動量保存式の式（２−５）を解いて、仮流速ｖpn,kを算出する（計算ステップ２３１）。
次に、流体計算手段２２は、式（２−５）及び式（２−４）を連立して得られるポアソン方程式を陰解法により解いて、圧力ｐn,kを算出する（計算ステップ２３２）。ここで
、流体計算手段２２は、ステップＳ２−１において更新されたヤコビアンＪT 及び、気泡成長計算手段２３から引渡された気泡半径Ｒn 、気泡成長速度Ｒ'n、ボイド率αn を用いる。また、流体計算手段２２は、この処理を行う直前に算出された密度ρを用いる。例えば、同一タイムステップｎではステップＳ２−３の処理が１回目の場合（ｋ＝１の場合）には、その前のタイムステップｎ−１で算出した密度ρn-1 を用いる。また、同一タイムステップｎでステップＳ２−３の処理が２回以降の場合（ｋ≧２の場合）には、その直前の処理ループで算出した密度ρn,k-1 を用いる。
【００５５】
次に、流体計算手段２２は、図８に示すように、流速の修正計算を行う（計算ステップ２３３）。具体的には、流体計算手段２２は、計算ステップ２３１で式（２−５）を陽解法を用いて算出した仮流速ｖpn,kと、計算ステップ２３２でポアソン方程式を陰解法を用いて算出した圧力ｐn,k とを用いて、式（２−５）の再計算を行う。
【００５６】
更に、流体計算手段２２は、質量保存式の計算を行う（計算ステップ２３４）。具体的には、ステップＳ２−１において更新されたヤコビアンＪT 、気泡成長計算手段２３から引渡された気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αn及び修正後の流速ｖn,kを用
いて、式（２−４）の再計算を行う。そして、この質量保存式の計算により密度ρn,kを
算出する。
【００５７】
次に、流体計算手段２２は、状態方程式の計算を行う（ステップ２−４）。具体的には、図９に示すように、式（２−４）及び式（２−５）を連立して得られるポアソン方程式を陰解法により解いて算出した圧力ｐn,kを式（２−６）に代入する。そして、式（２−
６）から得たＶ（Ｔｔ，ｐ）の値を式（２−７）を代入して、密度ρeos を算出する。
【００５８】
次に、流体計算手段２２は、算出した密度ρn,kを、このタイムステップｎの密度ρnとして設定するか否かを判断する（ステップＳ２−５）。具体的には、質量保存式から求めた密度ρn,kが状態方程式から算出した密度ρeosに対して収束したか否かの判定を行う（ステップＳ２−５）。すなわち、式（２−８）において算出した全メッシュの最大値が、収束判断値εより小さいか否かを算出する。
【００５９】
そして、流体計算手段２２は、いずれかのメッシュにおける式（２−８）により算出した値が収束判断値εより大きい場合には、ループ回数を１つ増加させて、ループ回数ｋをｋ＋１とする（ステップＳ２−６）。そして、ループ回数ｋ＋１とした状態で、上記ステップＳ２−３以降の処理を再度行う。これにより、密度ρn,k を用いて、ステップＳ２−３が再計算されて、新たな密度ρn,k+1が算出される。このように、収束判定を満たさな
い場合（ステップＳ２−５において「ＮＯ」の場合）には、新たに算出される密度ρn,k+1 についてステップＳ２−３を繰り返し行うことにより、収束条件を満たすような密度ρnを求める。
【００６０】
一方、全メッシュの式（２−８）の最大値が収束判断値εより小さい場合（ステップＳ２−５において「ＹＥＳ」の場合）には、タイムステップが終了タイムステップ数になったか否かを判断する（ステップＳ２−７）。具体的には、流体計算手段２２は、タイムステップ数をカウントし、予め定められた所定数に到達したかどうかを判断する。
【００６１】
そして、終了タイムステップ数に到達していない場合（ステップＳ２−７において「Ｎ
Ｏ」の場合）、流体計算手段２２は、算出した密度ρn,kを、このタイムステップｎの密
度ρnとして決定して記憶する（ステップＳ２−８）。なお、この記憶された密度ρn は
、次のタイムステップｎ＋１についての流体計算で用いられる。
【００６２】
更に、流体計算手段２２は、算出した圧力ｐnを固体計算手段２１に引渡す（ステップ
Ｓ２−９）。そして、固体計算手段２１は、受け取った圧力ｐnに基づいて次のタイムス
テップｎ＋１に関してステップＳ１−１以降の処理を行う。以上により、流体計算が完了する。
【００６３】
このようにして、コンピュータ２０は、各タイムステップｎについて、固体計算手段２１において構造物の内壁の変位を計算し、流体計算手段２２において流体である水銀の挙動の計算を行う。そして、シミュレーションによって状況を算出する時間分に対応するようなタイムステップｎを繰り返して行うことにより、水銀ターゲットの所定時間におけるシミュレーションが求められる。
【００６４】
本実施形態によれば、以下のような効果を得ることができる。
・ 本実施形態では、コンピュータ２０の固体計算手段２１は、初期条件又は流体計算手段２２から引渡された圧力に基づいて変位量ｕnを算出し、流体計算手段２２に引渡す
（ステップＳ１−４）。流体計算手段２２は、変位量ｕnに基づいて境界流速及びヤコビ
アンの計算を行い（ステップＳ２−１）、気泡成長計算を起動する（ステップＳ２−２）。この起動において流体計算手段２２は、流体の圧力ｐnをコンピュータ２０の気泡成長
計算手段２３に引渡す。気泡成長計算手段２３は、この圧力ｐnを用いてを時間発展方程
式を計算し（ステップＳ３−２）、算出した気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド
率αnを流体計算手段２２に引渡す（ステップＳ３−５）。流体計算手段２２は、気泡半
径Ｒ、気泡成長速度Ｒ'n及びボイド率αnを用いて、運動量保存式、ポアソン方程式、流
速補正及び質量保存式を計算し（ステップＳ２−３）、密度ρn,k 及び圧力ｐn,kを算出
する。算出した密度ρn,kが状態方程式から算出した密度ρeosに対して収束判定条件を満たさない場合には、算出された圧力ｐn,kを用いて、ステップＳ２−３以降処理を繰り返
して行う。一方、収束判定条件を満たす場合（ステップＳ２−５において「ＹＥＳ」の場合）には、算出された圧力ｐnを固体計算手段２１に引渡す。固体計算手段２１は、引渡
された圧力ｐnに基づいて変位量ｕn+1を計算する。従って、コンピュータ２０は、流体の気体成長計算から算出した気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを用いて、気泡の質量に関する項を含む質量保存式及び運動量保存式から流体の圧力ｐを算出する。このため、急激な圧力変化によって流体に発生するキャビテーションによる気泡を質量保存式に入れて、シミュレーションを行うことができる。従って、キャビテーションを考慮した流体（ここでは水銀）についての計算と、この流体を内包した構造物（ここでは容器）についての計算とを連成させたシミュレーションを、実際の現象に沿って行うことができる。
【００６５】
・ 本実施形態では、気泡成長計算手段２３は、気泡半径Ｒn、気泡成長速度Ｒ'n及び
ボイド率αnを算出する場合、タイムステップｎ≧２のときには、その前のタイムステッ
プｎ−１において算出した圧力ｐn-1を用いて時間発展方程式の式（３−４）を計算する
。このため、圧力ｐの変動に応じた気泡の成長を反映してシミュレーションを行うことができる。
【００６６】
・ 本実施形態では、気泡成長計算手段２３は、流体計算における１のタイムステップｎにおいて、必要ループの回数分だけ気泡成長計算のステップＳ３−２を繰り返し行う。これにより、流体計算におけるタイムステップの時間間隔よりも、より短い時間間隔で気泡成長の計算を行うことができる。従って、気泡成長に関する気泡半径Ｒn、気泡成長速
度Ｒ'n及びボイド率αnをより正確に得て、シミュレーションを行うことができる。
【００６７】
また、上記実施形態は以下のように変更してもよい。
○ 上記実施形態においては、コンピュータ２０は、すべてのメッシュに対して気泡成長計算を起動した上で、運動量保存式、ポアソン方程式、流速補正及び質量保存式の計算を行った（ステップＳ２−３）。これに代えて、気泡の成長の影響が少ないメッシュについては、気泡計算を省略してもよい。具体的には、図１０に示すように、流体計算手段２２が、境界流速とヤコビアンの計算をした後（ステップＳ２−１）、圧力ｐの時間微分が基準値ｋe以上か否かの判定を行う（ステップＳ５−１）。そして、圧力ｐの時間微分が
基準値ｋe以上の場合（ステップＳ５−１において「ＹＥＳ」の場合）にはステップＳ２
−２を行う。一方、圧力ｐの時間微分が基準値ｋeよりも小さい場合（ステップＳ５−１
において「ＮＯ」の場合）には、気泡成長計算をスキップして、運動量保存式、ポアソン方程式、流速補正及び質量保存式の計算を行う（ステップＳ２−３）。これにより、気泡に変化の少ないメッシュについては、気泡成長計算を省略できるので、コンピュータ２０の計算負荷を軽減することができる。なお、このメッシュについても、先のタイムステップで計算した気泡を考慮した質量保存式を用いて圧力を計算する。このため、このメッシュに関しても気泡を考慮したシミュレーションを行うことができる。
【００６８】
○ 上記実施形態においては、コンピュータは、式（２−８）の式が収束判断値εよりも大きい場合には、ループ回数ｋをｋ＋１としてステップＳ２−３以降の処理を繰り返して行った。これに代えて、図１１に示すように、式（２−８）の式が収束判断値εよりも大きい場合には、タイムステップの変更処理（ステップＳ６−１）を行ってから、ループ回数ｋの変更を行ってもよい。
【００６９】
この場合、流体計算手段２２に、変更を行なうための基準として、変更回数ｋrに関す
るデータを記憶させておく。そして、このタイムステップの変更処理（ステップＳ６−１）においては、流体計算手段２２は、まず、ループ回数ｋが変更回数ｋrであるか否かを
判断する。そして、ループ回数ｋが変更回数ｋrより少ない場合には、そのループ回数ｋ
をｋ＋１として（ステップＳ２−６）、上記ステップＳ２−３以降の処理を行う。
【００７０】
一方、ループ回数ｋが変更回数ｋr以上の場合には、気泡成長計算における気泡計算の
必要タイムステップｎb数を増加させて、気泡計算における時間間隔ｄｔbの値をより短い値に再設定し直す。例えば、タイムステップｎb数を１０倍にして、時間間隔ｄｔbを１／１０に設定し直す。そして、流体計算手段２２は、収束判断値εよりも式（２−８）の値が大きいと判断したときの圧力ｐn,kをタイムステップｎの圧力ｐnとして算出する。その後、流体計算手段２２は、次のタイムステップｎ＋１における気泡成長計算を行う場合には、気泡計算タイムステップｎbの設定し直した時間間隔ｄｔbの値を用いて算出する。従って、流体計算の収束状況に応じて気泡成長計算のパラメータを変更し、より的確な圧力ｐnを算出することができる。
【００７１】
また、気泡計算タイムステップｎbの時間間隔ｄｔbの値をより短い値に再設定し直した後、ステップＳ２−３〜Ｓ２−５のループ回数が変更回数ｋrよりも少ない状態で収束判
断値εを満たすような場合には、必要タイムステップｎb数を、元の値に戻してもよい。
具体的には、流体計算手段２２は、収束判定条件を満たすと判断した場合（ステップＳ２−５において「ＹＥＳ」の場合）のループ回数ｋが所定値以下であったときに、気泡成長計算手段２３が行う処理の必要ループ数の設定値を取得する。そして、この必要ループ数の設定値が初期値よりも多い場合には、流体計算手段２２は、その必要ループ数を元の初期値に設定し直す。
【００７２】
また、気泡計算の必要タイムステップｎb数の値を設定し直した場合に、この気泡計算
タイムステップｎbの時間間隔ｄｔbの値を用いて、そのときのタイムステップｎの気泡半
径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを算出し、そのタイムステップｎの圧力ｐnを
算出してもよい。すなわち、この場合には、図１１のステップＳ６−１の処理を行った後、ステップＳ２−２の処理を行うようにジャンプしてもよい。これにより、そのタイムステップｎの圧力ｐnを、より正確に算出して、より実際の現象に近いシミュレーションを
行うことができる。
【００７３】
更に、流体計算手段２２が、式（２−８）の値が収束判断値εを満たすために行ったループ回数ｋが所定の回数ｋｑ以上である場合に、次の気泡成長計算における必要タイムステップｎb数を増加させて、気泡計算における時間間隔ｄｔbの値をより短い値に再設定し直すようにしてもよい。これにより、収束判断値εにより収束しにくい状態には、ステップＳ２−３〜Ｓ２−５のループ回数を減少させて、これらの計算量を減少することができる。
【００７４】
○ 上記実施形態においては、固体計算手段２１と流体計算手段２２とは、１つのタイムステップを行う度に、構造物の変位量ｕの引渡しと、流体の圧力ｐの引渡しを行った。固体計算手段２１及び流体計算手段２２は、変位量ｕ及び圧力ｐの引渡しを同じ時間間隔のタイミングであれば、１つのタイムステップ毎に行わなくてもよい。例えば、固体計算におけるタイムステップを２回行うときの時間間隔と、流体計算におけるタイムステップを３回行うときの時間間隔とが同じである場合には、その最小公倍数である６回のタイムステップで、固体計算手段２１及び流体計算手段２２が、変位量ｕ及び圧力ｐの引渡しを行うようにしてもよい。
【００７５】
○ 上記実施形態においては、構造として水銀ターゲットの容器の内壁、流体として容器内部の水銀を対象とした。本発明は、これに限られることなく、急激な圧力変化によりキャビテーションが発生する流体と構造とが接している状態でのシミュレーションに用いることができる。
【図面の簡単な説明】
【００７６】
【図１】実施形態におけるシステムの概略図。
【図２】固体計算手段の内部構成を説明するための説明図。
【図３】流体計算手段の内部構成を説明するための説明図。
【図４】気泡成長計算手段の内部構成を説明するための説明図。
【図５】実施形態における処理手順を説明するための流れ図。
【図６】気泡成長計算の算出手順を説明するための流れ図。
【図７】運動量保存式及びポアソン方程式の算出手順を説明するための流れ図。
【図８】流速の修正計算及び質量保存式の算出手順を説明するための流れ図。
【図９】収束条件の処理手順を説明するための流れ図。
【図１０】第１変更例における処理手順を説明するための流れ図。
【図１１】第２変更例における処理手順を説明するための流れ図。
【符号の説明】
【００７７】
α…ボイド率、ρ…密度、ρeos…状態方程式から算出される密度、ｐ…圧力、Ｒ…気
泡半径、Ｒ’ …気泡成長速度、ｄｔb…気泡成長速度を算出するための時間間隔、ｖw…
境界流速、２０…コンピュータ、２１…固体計算手段、２２…流体計算手段、２３…気泡成長計算手段。
【技術分野】
【０００１】
本発明は、構造物の中の流体の挙動を解析するための流体構造連成解析方法及び流体構造連成解析プログラムに関する。
【背景技術】
【０００２】
構造物の中の流体の挙動を解析するために、各種シミュレーション技術が検討されている。例えば、核破砕タイプ中性子源の開発においても、計算機シミュレーションの検討が進められている。この開発では、加速器で加速された陽子をターゲットとなる物質に照射し、ターゲットに核破砕を起こすことで中性子を発生させる。ここでは、ターゲットには、ステンレス鋼の容器（構造物）と、その中を循環する液体水銀（流体）から構成された水銀ターゲットを用いる。そして、この水銀ターゲットの中の現象を予め予測するため計算機シミュレーションが大切である。この場合、構造物及び流体が相互に影響し合っているので、構造物及び流体に関する方程式を連立させた解析が必要である。このような構造物及び流体に関するシミュレーションを行なうための解析装置が開示されている（例えば、特許文献１参考。）。
【０００３】
特許文献１においては、流体解析による流体の解及び構造解析による構造全体の解のうち、流体と構造が接触する境界における流体の解と構造の解とが一致する共通解を求める。そして、この求められた共通解を用いて、再度、流体解析と構造解析を実行することにより安定解を算出している。
【特許文献１】特開２０００−３５２５４５号公報（第１頁）
【発明の開示】
【発明が解決しようとする課題】
【０００４】
ところで、上述の例のような水銀ターゲットの場合、内部で核破砕という高密度の熱の発生を伴う現象が生じるため、水銀ターゲットの容器には大きな力がかかる。この場合、高密度の熱の発生により液体水銀が急激に膨張し、その結果、生じる衝撃波によりターゲット外壁に大きな応力が発生する。そして、外壁に歪が生じ、水銀中にキャビテーションが発生する可能性がある。このキャビテーションは、外壁の劣化等の影響を与えるため、キャビテーションが発生した場合の外壁への影響の計算することが大切である。
【０００５】
しかし、このような状況下でのキャビテーションを考慮した流体構造連成解析プログラムはない。
本発明は、上述の課題に鑑みてなされ、その目的は、キャビテーションを考慮して、実際の挙動により近いシミュレーションを行うための流体構造連成解析方法及び流体構造連成解析プログラムを提供することにある。
【課題を解決するための手段】
【０００６】
上記問題点を解決するために、請求項１に記載の発明は、コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析方法であって、前記コンピュータが、運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出段階と、流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出段階と、前記圧力算出段階で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定段階と、前記誤差が所定範囲になった場合、前記算出した圧力を用いて
前記構造物の変位を算出する変位算出段階とを含むことを要旨とする。
【０００７】
請求項２に記載の発明は、請求項１に記載の流体構造連成解析方法において、前記コンピュータは、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算段階を更に含むことを要旨とする。
【０００８】
請求項３に記載の発明は、請求項２に記載の流体構造連成解析方法において、前記コンピュータは、有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断段階を更に含み、気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うことを要旨とする。
【０００９】
請求項４に記載の発明は、請求項２又は３に記載の流体構造連成解析方法において、前記コンピュータは、前記圧力算出段階において収束判定に至るまでの繰り返し回数を計数し、前記繰り返し回数が基準回数になっても、前記収束判定段階において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うことを要旨とする。
【００１０】
請求項５に記載の発明は、コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析プログラムであって、前記コンピュータを、運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出手段と、流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出手段、前記圧力算出手段で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定手段、及び前記誤差が所定範囲になった場合、前記算出した圧力を用いて前記構造物の変位を算出する変位算出手段として機能させることを要旨とする。
【００１１】
請求項６に記載の発明は、請求項５に記載の流体構造連成解析プログラムにおいて、前記コンピュータを、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算手段として更に機能させることを要旨とする。
【００１２】
請求項７に記載の発明は、請求項６に記載の流体構造連成解析プログラムにおいて、前記コンピュータを、有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断手段として更に機能させ、気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うようにさせたことを要旨とする。
【００１３】
請求項８に記載の発明は、請求項６又は７に記載の流体構造連成解析プログラムにおいて、前記コンピュータを、前記圧力算出手段において収束判定に至るまでの繰り返し回数を計数し、前記繰り返し回数が基準回数になっても、前記収束判定手段において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うようにさせたことを要旨とする。
【００１４】
（作用）
請求項１又は５に記載の発明によれば、コンピュータは、流体の気泡成長計算から算出
した気泡半径、気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式から、流体の圧力を算出する。このため、急激な圧力変化によって流体に発生するキャビテーションによる気泡を質量保存式に入れて、シミュレーションを行うことができる。従って、キャビテーションを考慮した流体についての計算と、この流体を内包した構造物についての計算とを連成させたシミュレーションを、実際の現象により近づけて行うことができる。
【００１５】
請求項２又は６に記載の発明によれば、コンピュータは、流体の質量保存式及び運動量保存式から圧力が算出された場合には、その圧力を用いて気泡成長計算を行う。このため、圧力変動に応じた気泡の成長を反映したシミュレーションを行うことができる。
【００１６】
請求項３又は７に記載の発明によれば、コンピュータは、有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断し、大きいメッシュについてのみ気泡成長計算を行う。従って、気泡の成長の影響が少ないメッシュについては、気泡成長計算を行わないので、計算量を少なくすることができる。また、この場合であっても、気泡の成長が大きくないメッシュについて、気泡の質量に関する項を含む質量保存式と運動量保存式とから圧力を算出するので、実際の現象に近いシミュレーションを行うことができる。
【００１７】
請求項４又は８に記載の発明によれば、コンピュータは、密度が所定範囲の誤差にならなかった場合には、算出された圧力を用いて再計算を行う。また、コンピュータは、再計算を所定回数行っても、所定範囲の誤差にならない場合には、時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて圧力を再計算する。このため、気泡について、より正確なシミュレーションができるので、実際の現象に近いシミュレーションを行うことができる。
【発明の効果】
【００１８】
本発明によれば、キャビテーションを考慮して、構造物及び流体の挙動に関するシミュレーションを行うことができる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明を具体化した一実施形態を図１〜図９に基づいて説明する。本実施形態においては、水銀ターゲットにおける圧力波伝播のシミュレーションに適用する。ここでは、水銀ターゲットの容器を構造物、その中の水銀を流体の挙動としたときの式を連立して解析を行う。
【００２０】
このシミュレーションを行うコンピュータ２０は、図１に示すように、シミュレーションの結果を出力するための出力手段３０、各指示を入力するための入力手段４０や初期値や設定値に関するデータを記憶するデータ記憶手段５０に接続されている。
【００２１】
出力手段３０は、シミュレーション結果を出力するために用いられ、具体的には、ディスプレイやデータ記憶部などを用いる。入力手段４０は、シミュレーションを行うために各種指示を入力するために用いられ、具体的には、キーボードやマウスなどである。データ記憶手段５０は、シミュレーションを行うために必要なデータを予め記憶させるために用いられる。このデータ記憶手段５０には、例えば、解析を行う固体及び流体の形状やマトリクスの形状、解析に用いられる物性値などに関するデータが記憶されている。
【００２２】
コンピュータ２０は、固体計算手段２１、流体計算手段２２及び気泡成長計算手段２３を備える。コンピュータ２０は、これらを用いて境界流速算出段階、圧力算出段階、収束判定段階、変位算出段階、気泡成長計算段階及び判断段階を行う。コンピュータ２０は、
流体構造連成解析プログラムを実行することにより、境界流速算出手段、圧力算出手段、収束判定手段、変位算出手段、気泡成長計算手段及び判断手段として機能する。ここで、各手段の機能を図２〜図４を用いて説明する。
【００２３】
図２に示すように、固体計算手段２１は、構造物（固体）に関して、式（１−１）、式（１−２）及び式（１−３）の計算処理を実行し、構造物の変位量ｕを算出する。
式（１−１）は、力学的境界条件式であり、構造物が、その周囲の流体から受けた圧力ｐを節点荷重ｔw に変換するための計算式である。なお、この式（１−１）において、ｖaは定数である。
【００２４】
式（１−２）は、応力テンソルとひずみテンソルの構造物についての関係式であり、節点荷重ｔw を応力テンソルσに変換するための計算式である。この式（１−２）において、Ｖv は法線ベクトルである。
【００２５】
式（１−３）は、構造物に関する運動量保存式であり、式（１−２）で算出した応力テンソルσを用いて、構造物の変位量ｕを算出するための計算式である。この式（１−３）において、ρsは固体の密度、∇は微分演算子、ｂfは外力である。なお、本実施形態では、外力ｂfを「０」と設定する。
【００２６】
図３に示すように、流体計算手段２２は、流体に関して、式（２−１）〜式（２−８）の計算処理を実行して、流体の圧力ｐを算出する。
式（２−１）は、固体と流体の境界面（カップリングサーフェース）の座標値Ｘを算出するための計算式である。この式（２−１）において、Ｘnew は新たに算出する境界面の座標値であり、Ｘold は座標値Ｘnew を算出する直前の境界面の座標値である。
【００２７】
式（２−２）は、固体と流体の境界面の境界流速ｖw を算出するための計算式である。ここで、ｄｔはクーラン条件から決まる時間間隔である。
式（２−３）は、直交座標系（ｘ1、ｘ2、ｘ3、ｔ）を一般座標系（ξ1、ξ2、ξ3、Ｔ）に変換するヤコビアンＪT を算出するための計算式である。
【００２８】
式（２−４）は、流体における質量保存式である。この式（２−４）においては、後述するバブルモデルを組み込んでいる。この式（２−４）において、ρは流体の密度、Ｕi
は流速ベクトル（の反変成分）、πは円周率、ｍは単位体積あたりの気泡数、Ｒは気泡
半径、Ｒ'（＝dR/dt）は気泡の半径の時間微分値であり、いわゆる気泡成長速度である。αは、流体におけるボイド率であり、後述する気泡成長計算手段２３によって算出される。
【００２９】
式（２−５）は、流体の運動量保存式である。この式（２−５）において、ｖは速度、∇は微分演算子、Ｆは運動量保存式に現れるその他の項を示している。この式（２−５）及び式（２−４）の２式は、連立することにより流体に関する圧力に関するポアソン方程式となり、これらから流体計算手段２２は圧力ｐを算出する。
【００３０】
式（２−６）は、流体における状態方程式である。この式（２−６）は、本実施形態ではＴａｉｔの方程式であり、圧力ｐに基づいて単位質量あたりの体積Ｖを算出する。この式（２−６）においてＴtは流体の温度であり、初期条件として与えられる値である。ま
た、式（２−６）において、Ｖ（Ｔt、０）＝Ａ・Ｔt ＋Ｂ、Ｂ（Ｔt）＝Ｂ０・Ｔt＋Ｂ1、Ｃ（Ｔt）＝Ｃ０・Ｔt＋Ｃ1である。ここで、Ａ、Ｂ、Ｃはパラメータであり、Ａ、Ｂ
、Ｃ、Ｂ1、Ｃ1の値は、初期値として与えられる。
【００３１】
式（２−７）は、式（２−６）から算出された単位質量あたりの体積Ｖから、状態方程
式に基づく密度ρeos を算出するための式である。
式（２−８）は、収束判定を行うためのパラメータの算出式である。この式（２−８）の値が収束判断値εより小さくなった場合に、流体計算手段２２は、この式（２−４）及び式（２−５）から算出した値を圧力ｐとして決定する。
【００３２】
図４に示すように、気泡成長計算手段２３は、式（３−１）〜式（３−５）の計算処理を実行して、流体に含まれる気泡の成長を計算する。
式（３−１）は、気泡に関する圧力の関係式である。この式（３−１）において、ｐt は気泡内部圧力である。ｐt'（＝ｄｐt／ｄｔb）は内部圧力変化であり、気泡内部圧力ｐtの時間の１階微分である。γは比熱比、Ｋｈは熱伝達係数、Ｔtbは気泡内温度である。
τは、熱伝達係数Ｋｈを温度θで気泡内温度Ｔtbまで積分した値であり、便宜上、ここでは熱伝達積分係数と称する。
【００３３】
式（３−２）は、気泡に関する温度の式である。この式（３−２）と式（３−１）とを連立して解くことにより、気泡内部圧力ｐtを算出する。この式（３−２）において、Ｄ
は温度拡散係数であり、ｙは気泡半径Ｒによって正規化された気泡の外壁から中心点までの距離である。
【００３４】
式（３−３）は、気泡壁の外側の圧力ｐBを算出するための式である。この式（３−３
）には、上記式（３−１）及び式（３−２）を連立して算出した気泡内部圧力ｐtが代入
される。この式（３−３）において、σstは表面張力、μは流体の粘性係数である。
【００３５】
式（３−４）は、気泡半径Ｒ及び気泡成長速度Ｒ'を算出するための時間発展方程式で
ある。この式（３−４）において、Ｒ”（＝d2Ｒ/dt2 ）は気泡半径の時間の２階微分、
ｃは流体の音速である。ｐ∞は流体の静水圧であり、定数パラメータとして設定される値である。
式（３−５）は、ボイド率αを算出するための式である。
【００３６】
従って、コンピュータ２０は、図１に示すように、固体計算手段２１は流体計算手段２２に、算出した構造物の変位量ｕを引渡し、流体計算手段２２は固体計算手段２１に、流体計算手段２２で算出した圧力ｐを引渡す。また、流体計算手段２２は気泡成長計算手段２３に圧力ｐを引渡し、気泡成長計算手段２３は流体計算手段２２に、気泡成長計算手段２３で算出した気泡半径Ｒ、気泡成長速度Ｒ’、ボイド率αを引渡す。
【００３７】
次に、コンピュータ２０の計算処理について、図５〜図９を用いて説明する。ここでは、各記号に付された上付き文字「ｎ」は、タイムステップｎにおける値であることを示している。このタイムステップｎとは、シミュレーションにより構造物及び流体の状態を求める時間をクーラン条件から決まる時間間隔ｄｔで分割したときの時間回数である。また、後述する流体計算において、各記号に付された上付き文字「ｋ」は、１つのタイムステップｎにおける処理ループの回数を示す。
【００３８】
まず、計算処理において、初期条件として、節点荷重ｔw 、有限要素のメッシュ分割の個数、固体の材料の物性、流体の物性、設定温度等に関するデータ、終了タイムステップ数などの値が入力手段４０を介して入力されて、データ記憶手段５０に記憶される。本実施形態では、次の初期条件が設定される。まず、初期の壁面は大気圧と内部の流体の圧力がつりあっていると考えて、節点荷重ｔwに「０」の値が設定される。また、節点荷重ｔw
以外の外力はないと仮定して外力ｂf＝０の値が設定される。更に、通常は各メッシュ中には複数の気泡が含まれているが、本実施形態においては、各メッシュ中に含まれる気泡群を平均化して、１つの気泡の成長起動として代表させるため、単位体積あたりの気泡数ｍを指定する。また、流体計算における１つのタイムステップに対して、気泡成長計算手
段２３が気泡成長の計算を行うときの気泡計算タイムステップｎb数（以下、必要ループ
数と称する）が設定される。更に、水銀ターゲットのシミュレーションにおいては、各メッシュの流体に作用する圧力が初期条件として与えられる。また、本実施形態では内壁の変位量ｕを考慮するために、収束判断値εは１．０×１０-7に設定している。
【００３９】
（固体計算）
まず、コンピュータ２０が固体のシミュレーション計算処理を実行する場合を説明する。具体的には、コンピュータ２０の固体計算手段２１が、まず、力学的境界条件の設定を行う（ステップＳ１−１）。タイムステップｎ＝１の場合、処理の開始の指示を受けた固体計算手段２１は、力学的境界条件として、初期条件として入力された節点荷重ｔwn=1
を用いる。
【００４０】
また、タイムステップｎ≧２のときには、固体計算手段２１は、後述するように流体計算手段２２において算出された流体の圧力ｐの引渡しを受ける。この場合、ステップＳ１−１における力学的境界条件として、後述する流体計算から算出される流体の圧力ｐを式（１−１）に代入して、節点荷重ｔwn を求める。
【００４１】
そして、設定された節点荷重ｔwnを式（１−２）に代入して、固体のカップリングサーフェース上の応力テンソルσを求める。
次に、固体計算手段２１は、算出した応力テンソルσnから式（１−３）を用いて、固
体におけるカップリングサーフェースの変位量ｕnを算出する（ステップＳ１−２）。ま
た、式（１−３）の固体の密度ρsとして、初期条件としてデータ記憶手段５０に記憶さ
れている固体の材料の物性値を、固体計算手段２１は取得して用いる。
【００４２】
そして、固体計算手段２１は、タイムステップが終了タイムステップ数になったか否かを判断する（ステップＳ１−３）。具体的には、固体計算手段２１は、タイムステップ数をカウントし、予め定められた所定数に到達したかどうかを判断する。終了タイムステップ数が終了していない場合（ステップＳ１−３において「ＮＯ」の場合）、固体計算手段２１は、構造物（内壁）の変位量ｕn を流体計算手段２２に引渡す（ステップＳ１−４）。
【００４３】
（流体計算の第１プロセス）
次に、流体計算について説明する。ここでは、流体計算は、気泡成長計算を入れ子構造にしている。そこで、気泡成長計算の前のプロセスを「第１プロセス」、気泡成長計算の後のプロセスを「第２プロセス」と呼ぶ。
【００４４】
固体計算手段２１から構造物の変位量ｕnを引渡された流体計算手段２２は、境界流速
とヤコビアンの計算を行う（ステップＳ２−１）。具体的には、流体計算手段２２は、固体計算手段２１において算出された構造物の変位量ｕnを、式（２−１）及び式（２−２
）に代入し、構造物である内壁と流体である水銀の境界面（カップリングサーフェース）の座標値Ｘn 及び境界流速ｖwnを算出する。ここで、タイムステップｎ＝１の場合、壁を変位量ｕn=1を計算していないため、座標値Ｘold は、初期条件による値、すなわち流体
の各メッシュ要素の初期座標値を用いる。
そして、流体計算手段２２は、算出した流体の座標値Ｘn 及び境界流速ｖwnに基づいて、ヤコビアンＪTを更新する。
【００４５】
次に、流体計算手段２２は、気泡成長計算を起動する（ステップＳ２−２）。具体的には、圧力ｐnの値を気泡成長計算手段２３に引渡す。ここで、タイムステップｎ＝１のと
きには、流体計算手段２２は、初期条件で設定された圧力ｐn=1を気泡成長計算手段２３
に引渡す。また、タイムステップｎ≧２のときには、後述するステップにおいて流体計算
手段２２が算出した圧力ｐ、すなわちタイムステップ「ｎ−１」の圧力ｐn-1の値を気泡
成長計算手段２３に引渡す。
【００４６】
（気泡成長計算）
気泡成長計算においては、コンピュータ２０の気泡成長計算手段２３は、気泡の周囲の流体の圧力ｐnの更新を行う（ステップＳ３−１）。具体的には、初期条件としては、気
泡の周囲の圧力ｐn=1には、設定された流体計算における初期圧力を用いる。また、初期
条件として気泡内部圧力ｐtには、静的平衡状態にあると仮定して、気泡壁の外側の圧力
ｐB＝ｐ＋ｐ∞とした式（３−３）を用いて算出する。また、タイムステップｎ≧２のと
きには、タイムステップｎ−１で算出した流体の圧力ｐn-1を用いる。
【００４７】
次に、気泡成長計算手段２３は、気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを算出する。これは、ステップＳ３−２を予め設定された必要ループ回、繰り返し行うことにより算出する。
【００４８】
次に、ステップＳ３−２の繰り返し処理について図６を用いて詳述する。ここでは、各記号に付された上付き文字「ｊ」は、１つの気泡計算タイムステップｎbにおける処理ル
ープの回数を示す。なお、この気泡計算タイムステップｎbは、１つのタイムステップｎ
において必要ループ数、行われる。
【００４９】
そこで、気泡成長計算手段２３は、タイムステップｎの気泡半径Ｒjについての時間発
展方程式の計算を行う（ステップＳ３−２）。具体的には、気泡成長計算手段２３は、まず、式（３−１）及び式（３−２）を連立させて解くことにより、気泡内部圧力ｐtjと、熱伝達積分係数τjが求まり、この熱伝達積分係数τjから気泡内部温度Ｔtbjが求まる。
この算出した気泡内部圧力ｐtjを式（３−３）に代入して、気泡壁の外側の圧力ｐBjを算出する。この算出した気泡壁の外側の圧力ｐBjを用いて式（３−４）を解く。これにより、この式（３−１）の解として、気泡半径Ｒj及び気泡成長速度Ｒ'j（＝ｄＲj/ｄｔb ）
が得られる。
【００５０】
次に、気泡成長計算手段２３は、予め設定された必要ループを終了したか否かを判断する（ステップＳ３−３）。そして、気泡成長計算手段２３は、予め設定された必要ループが終了していない場合（ステップＳ３−３において「ＮＯ」の場合）には、算出した気泡半径Ｒjと、気泡成長速度Ｒ'jとを用いて、ステップＳ３−２の再計算を行う。
【００５１】
一方、気泡成長計算手段２３は、必要ループが終了した（ステップＳ３−３において「ＹＥＳ」）と判断すると、ボイド率αjの設定を行う（ステップＳ３−４）。具体的には
、気泡成長計算手段２３は、算出した気泡半径Ｒjを式（３−５）に代入することにより
、ボイド率αjを算出する。
【００５２】
次に、気泡成長計算手段２３は、算出した気泡半径Ｒj、気泡成長速度Ｒ'jを気泡半径
Ｒn、気泡成長速度Ｒ'nとして、またボイド率αnを流体計算手段２２に引渡す（ステップＳ３−５）。すなわち、気泡成長計算手段２３は、必要ループに応じてステップＳ３−２を繰り返し行うことにより、より正確な気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを算出する。
【００５３】
（流体計算の第２プロセス）
気泡成長計算手段２３から気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを引渡された流体計算手段２２は、運動量保存式、ポアソン方程式、流速補正及び質量保存式の計算を行う（ステップＳ２−３）。この計算は、３次元座標で表現される全メッシュについて行われる。この計算について図７〜図９を用いて詳述する。
【００５４】
図７に示すように、流体計算手段２２は、まず、陽解法を用いて運動量保存式の式（２−５）を解いて、仮流速ｖpn,kを算出する（計算ステップ２３１）。
次に、流体計算手段２２は、式（２−５）及び式（２−４）を連立して得られるポアソン方程式を陰解法により解いて、圧力ｐn,kを算出する（計算ステップ２３２）。ここで
、流体計算手段２２は、ステップＳ２−１において更新されたヤコビアンＪT 及び、気泡成長計算手段２３から引渡された気泡半径Ｒn 、気泡成長速度Ｒ'n、ボイド率αn を用いる。また、流体計算手段２２は、この処理を行う直前に算出された密度ρを用いる。例えば、同一タイムステップｎではステップＳ２−３の処理が１回目の場合（ｋ＝１の場合）には、その前のタイムステップｎ−１で算出した密度ρn-1 を用いる。また、同一タイムステップｎでステップＳ２−３の処理が２回以降の場合（ｋ≧２の場合）には、その直前の処理ループで算出した密度ρn,k-1 を用いる。
【００５５】
次に、流体計算手段２２は、図８に示すように、流速の修正計算を行う（計算ステップ２３３）。具体的には、流体計算手段２２は、計算ステップ２３１で式（２−５）を陽解法を用いて算出した仮流速ｖpn,kと、計算ステップ２３２でポアソン方程式を陰解法を用いて算出した圧力ｐn,k とを用いて、式（２−５）の再計算を行う。
【００５６】
更に、流体計算手段２２は、質量保存式の計算を行う（計算ステップ２３４）。具体的には、ステップＳ２−１において更新されたヤコビアンＪT 、気泡成長計算手段２３から引渡された気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αn及び修正後の流速ｖn,kを用
いて、式（２−４）の再計算を行う。そして、この質量保存式の計算により密度ρn,kを
算出する。
【００５７】
次に、流体計算手段２２は、状態方程式の計算を行う（ステップ２−４）。具体的には、図９に示すように、式（２−４）及び式（２−５）を連立して得られるポアソン方程式を陰解法により解いて算出した圧力ｐn,kを式（２−６）に代入する。そして、式（２−
６）から得たＶ（Ｔｔ，ｐ）の値を式（２−７）を代入して、密度ρeos を算出する。
【００５８】
次に、流体計算手段２２は、算出した密度ρn,kを、このタイムステップｎの密度ρnとして設定するか否かを判断する（ステップＳ２−５）。具体的には、質量保存式から求めた密度ρn,kが状態方程式から算出した密度ρeosに対して収束したか否かの判定を行う（ステップＳ２−５）。すなわち、式（２−８）において算出した全メッシュの最大値が、収束判断値εより小さいか否かを算出する。
【００５９】
そして、流体計算手段２２は、いずれかのメッシュにおける式（２−８）により算出した値が収束判断値εより大きい場合には、ループ回数を１つ増加させて、ループ回数ｋをｋ＋１とする（ステップＳ２−６）。そして、ループ回数ｋ＋１とした状態で、上記ステップＳ２−３以降の処理を再度行う。これにより、密度ρn,k を用いて、ステップＳ２−３が再計算されて、新たな密度ρn,k+1が算出される。このように、収束判定を満たさな
い場合（ステップＳ２−５において「ＮＯ」の場合）には、新たに算出される密度ρn,k+1 についてステップＳ２−３を繰り返し行うことにより、収束条件を満たすような密度ρnを求める。
【００６０】
一方、全メッシュの式（２−８）の最大値が収束判断値εより小さい場合（ステップＳ２−５において「ＹＥＳ」の場合）には、タイムステップが終了タイムステップ数になったか否かを判断する（ステップＳ２−７）。具体的には、流体計算手段２２は、タイムステップ数をカウントし、予め定められた所定数に到達したかどうかを判断する。
【００６１】
そして、終了タイムステップ数に到達していない場合（ステップＳ２−７において「Ｎ
Ｏ」の場合）、流体計算手段２２は、算出した密度ρn,kを、このタイムステップｎの密
度ρnとして決定して記憶する（ステップＳ２−８）。なお、この記憶された密度ρn は
、次のタイムステップｎ＋１についての流体計算で用いられる。
【００６２】
更に、流体計算手段２２は、算出した圧力ｐnを固体計算手段２１に引渡す（ステップ
Ｓ２−９）。そして、固体計算手段２１は、受け取った圧力ｐnに基づいて次のタイムス
テップｎ＋１に関してステップＳ１−１以降の処理を行う。以上により、流体計算が完了する。
【００６３】
このようにして、コンピュータ２０は、各タイムステップｎについて、固体計算手段２１において構造物の内壁の変位を計算し、流体計算手段２２において流体である水銀の挙動の計算を行う。そして、シミュレーションによって状況を算出する時間分に対応するようなタイムステップｎを繰り返して行うことにより、水銀ターゲットの所定時間におけるシミュレーションが求められる。
【００６４】
本実施形態によれば、以下のような効果を得ることができる。
・ 本実施形態では、コンピュータ２０の固体計算手段２１は、初期条件又は流体計算手段２２から引渡された圧力に基づいて変位量ｕnを算出し、流体計算手段２２に引渡す
（ステップＳ１−４）。流体計算手段２２は、変位量ｕnに基づいて境界流速及びヤコビ
アンの計算を行い（ステップＳ２−１）、気泡成長計算を起動する（ステップＳ２−２）。この起動において流体計算手段２２は、流体の圧力ｐnをコンピュータ２０の気泡成長
計算手段２３に引渡す。気泡成長計算手段２３は、この圧力ｐnを用いてを時間発展方程
式を計算し（ステップＳ３−２）、算出した気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド
率αnを流体計算手段２２に引渡す（ステップＳ３−５）。流体計算手段２２は、気泡半
径Ｒ、気泡成長速度Ｒ'n及びボイド率αnを用いて、運動量保存式、ポアソン方程式、流
速補正及び質量保存式を計算し（ステップＳ２−３）、密度ρn,k 及び圧力ｐn,kを算出
する。算出した密度ρn,kが状態方程式から算出した密度ρeosに対して収束判定条件を満たさない場合には、算出された圧力ｐn,kを用いて、ステップＳ２−３以降処理を繰り返
して行う。一方、収束判定条件を満たす場合（ステップＳ２−５において「ＹＥＳ」の場合）には、算出された圧力ｐnを固体計算手段２１に引渡す。固体計算手段２１は、引渡
された圧力ｐnに基づいて変位量ｕn+1を計算する。従って、コンピュータ２０は、流体の気体成長計算から算出した気泡半径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを用いて、気泡の質量に関する項を含む質量保存式及び運動量保存式から流体の圧力ｐを算出する。このため、急激な圧力変化によって流体に発生するキャビテーションによる気泡を質量保存式に入れて、シミュレーションを行うことができる。従って、キャビテーションを考慮した流体（ここでは水銀）についての計算と、この流体を内包した構造物（ここでは容器）についての計算とを連成させたシミュレーションを、実際の現象に沿って行うことができる。
【００６５】
・ 本実施形態では、気泡成長計算手段２３は、気泡半径Ｒn、気泡成長速度Ｒ'n及び
ボイド率αnを算出する場合、タイムステップｎ≧２のときには、その前のタイムステッ
プｎ−１において算出した圧力ｐn-1を用いて時間発展方程式の式（３−４）を計算する
。このため、圧力ｐの変動に応じた気泡の成長を反映してシミュレーションを行うことができる。
【００６６】
・ 本実施形態では、気泡成長計算手段２３は、流体計算における１のタイムステップｎにおいて、必要ループの回数分だけ気泡成長計算のステップＳ３−２を繰り返し行う。これにより、流体計算におけるタイムステップの時間間隔よりも、より短い時間間隔で気泡成長の計算を行うことができる。従って、気泡成長に関する気泡半径Ｒn、気泡成長速
度Ｒ'n及びボイド率αnをより正確に得て、シミュレーションを行うことができる。
【００６７】
また、上記実施形態は以下のように変更してもよい。
○ 上記実施形態においては、コンピュータ２０は、すべてのメッシュに対して気泡成長計算を起動した上で、運動量保存式、ポアソン方程式、流速補正及び質量保存式の計算を行った（ステップＳ２−３）。これに代えて、気泡の成長の影響が少ないメッシュについては、気泡計算を省略してもよい。具体的には、図１０に示すように、流体計算手段２２が、境界流速とヤコビアンの計算をした後（ステップＳ２−１）、圧力ｐの時間微分が基準値ｋe以上か否かの判定を行う（ステップＳ５−１）。そして、圧力ｐの時間微分が
基準値ｋe以上の場合（ステップＳ５−１において「ＹＥＳ」の場合）にはステップＳ２
−２を行う。一方、圧力ｐの時間微分が基準値ｋeよりも小さい場合（ステップＳ５−１
において「ＮＯ」の場合）には、気泡成長計算をスキップして、運動量保存式、ポアソン方程式、流速補正及び質量保存式の計算を行う（ステップＳ２−３）。これにより、気泡に変化の少ないメッシュについては、気泡成長計算を省略できるので、コンピュータ２０の計算負荷を軽減することができる。なお、このメッシュについても、先のタイムステップで計算した気泡を考慮した質量保存式を用いて圧力を計算する。このため、このメッシュに関しても気泡を考慮したシミュレーションを行うことができる。
【００６８】
○ 上記実施形態においては、コンピュータは、式（２−８）の式が収束判断値εよりも大きい場合には、ループ回数ｋをｋ＋１としてステップＳ２−３以降の処理を繰り返して行った。これに代えて、図１１に示すように、式（２−８）の式が収束判断値εよりも大きい場合には、タイムステップの変更処理（ステップＳ６−１）を行ってから、ループ回数ｋの変更を行ってもよい。
【００６９】
この場合、流体計算手段２２に、変更を行なうための基準として、変更回数ｋrに関す
るデータを記憶させておく。そして、このタイムステップの変更処理（ステップＳ６−１）においては、流体計算手段２２は、まず、ループ回数ｋが変更回数ｋrであるか否かを
判断する。そして、ループ回数ｋが変更回数ｋrより少ない場合には、そのループ回数ｋ
をｋ＋１として（ステップＳ２−６）、上記ステップＳ２−３以降の処理を行う。
【００７０】
一方、ループ回数ｋが変更回数ｋr以上の場合には、気泡成長計算における気泡計算の
必要タイムステップｎb数を増加させて、気泡計算における時間間隔ｄｔbの値をより短い値に再設定し直す。例えば、タイムステップｎb数を１０倍にして、時間間隔ｄｔbを１／１０に設定し直す。そして、流体計算手段２２は、収束判断値εよりも式（２−８）の値が大きいと判断したときの圧力ｐn,kをタイムステップｎの圧力ｐnとして算出する。その後、流体計算手段２２は、次のタイムステップｎ＋１における気泡成長計算を行う場合には、気泡計算タイムステップｎbの設定し直した時間間隔ｄｔbの値を用いて算出する。従って、流体計算の収束状況に応じて気泡成長計算のパラメータを変更し、より的確な圧力ｐnを算出することができる。
【００７１】
また、気泡計算タイムステップｎbの時間間隔ｄｔbの値をより短い値に再設定し直した後、ステップＳ２−３〜Ｓ２−５のループ回数が変更回数ｋrよりも少ない状態で収束判
断値εを満たすような場合には、必要タイムステップｎb数を、元の値に戻してもよい。
具体的には、流体計算手段２２は、収束判定条件を満たすと判断した場合（ステップＳ２−５において「ＹＥＳ」の場合）のループ回数ｋが所定値以下であったときに、気泡成長計算手段２３が行う処理の必要ループ数の設定値を取得する。そして、この必要ループ数の設定値が初期値よりも多い場合には、流体計算手段２２は、その必要ループ数を元の初期値に設定し直す。
【００７２】
また、気泡計算の必要タイムステップｎb数の値を設定し直した場合に、この気泡計算
タイムステップｎbの時間間隔ｄｔbの値を用いて、そのときのタイムステップｎの気泡半
径Ｒn、気泡成長速度Ｒ'n及びボイド率αnを算出し、そのタイムステップｎの圧力ｐnを
算出してもよい。すなわち、この場合には、図１１のステップＳ６−１の処理を行った後、ステップＳ２−２の処理を行うようにジャンプしてもよい。これにより、そのタイムステップｎの圧力ｐnを、より正確に算出して、より実際の現象に近いシミュレーションを
行うことができる。
【００７３】
更に、流体計算手段２２が、式（２−８）の値が収束判断値εを満たすために行ったループ回数ｋが所定の回数ｋｑ以上である場合に、次の気泡成長計算における必要タイムステップｎb数を増加させて、気泡計算における時間間隔ｄｔbの値をより短い値に再設定し直すようにしてもよい。これにより、収束判断値εにより収束しにくい状態には、ステップＳ２−３〜Ｓ２−５のループ回数を減少させて、これらの計算量を減少することができる。
【００７４】
○ 上記実施形態においては、固体計算手段２１と流体計算手段２２とは、１つのタイムステップを行う度に、構造物の変位量ｕの引渡しと、流体の圧力ｐの引渡しを行った。固体計算手段２１及び流体計算手段２２は、変位量ｕ及び圧力ｐの引渡しを同じ時間間隔のタイミングであれば、１つのタイムステップ毎に行わなくてもよい。例えば、固体計算におけるタイムステップを２回行うときの時間間隔と、流体計算におけるタイムステップを３回行うときの時間間隔とが同じである場合には、その最小公倍数である６回のタイムステップで、固体計算手段２１及び流体計算手段２２が、変位量ｕ及び圧力ｐの引渡しを行うようにしてもよい。
【００７５】
○ 上記実施形態においては、構造として水銀ターゲットの容器の内壁、流体として容器内部の水銀を対象とした。本発明は、これに限られることなく、急激な圧力変化によりキャビテーションが発生する流体と構造とが接している状態でのシミュレーションに用いることができる。
【図面の簡単な説明】
【００７６】
【図１】実施形態におけるシステムの概略図。
【図２】固体計算手段の内部構成を説明するための説明図。
【図３】流体計算手段の内部構成を説明するための説明図。
【図４】気泡成長計算手段の内部構成を説明するための説明図。
【図５】実施形態における処理手順を説明するための流れ図。
【図６】気泡成長計算の算出手順を説明するための流れ図。
【図７】運動量保存式及びポアソン方程式の算出手順を説明するための流れ図。
【図８】流速の修正計算及び質量保存式の算出手順を説明するための流れ図。
【図９】収束条件の処理手順を説明するための流れ図。
【図１０】第１変更例における処理手順を説明するための流れ図。
【図１１】第２変更例における処理手順を説明するための流れ図。
【符号の説明】
【００７７】
α…ボイド率、ρ…密度、ρeos…状態方程式から算出される密度、ｐ…圧力、Ｒ…気
泡半径、Ｒ’ …気泡成長速度、ｄｔb…気泡成長速度を算出するための時間間隔、ｖw…
境界流速、２０…コンピュータ、２１…固体計算手段、２２…流体計算手段、２３…気泡成長計算手段。
【特許請求の範囲】
【請求項１】
コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析方法であって、
前記コンピュータが、
運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出段階と、
流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出段階と、
前記圧力算出段階で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定段階と、
前記誤差が所定範囲になった場合、前記算出した圧力を用いて前記構造物の変位を算出する変位算出段階と
を含むことを特徴とする流体構造連成解析方法。
【請求項２】
前記コンピュータは、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算段階を更に含むことを特徴とする請求項１に記載の流体構造連成解析方法。
【請求項３】
前記コンピュータは、
有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断段階を更に含み、
気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うことを特徴とする請求項２に記載の流体構造連成解析方法。
【請求項４】
前記コンピュータは、
前記圧力算出段階において収束判定に至るまでの繰り返し回数を計数し、
前記繰り返し回数が基準回数になっても、前記収束判定段階において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うことを特徴とする請求項２又は３に記載の流体構造連成解析方法。
【請求項５】
コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析プログラムであって、
前記コンピュータを、
運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出手段と、
流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出手段、
前記圧力算出手段で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定手段、及び
前記誤差が所定範囲になった場合、前記算出した圧力を用いて前記構造物の変位を算出する変位算出手段
として機能させることを特徴とする流体構造連成解析プログラム。
【請求項６】
前記コンピュータを、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算手段として更に機能させることを特徴とする請求項５に記載の流体構造連成解析プログラム。
【請求項７】
前記コンピュータを、
有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断手段として更に機能させ、
気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うようにさせたことを特徴とする請求項６に記載の流体構造連成解析プログラム。
【請求項８】
前記コンピュータを、
前記圧力算出手段において収束判定に至るまでの繰り返し回数を計数し、
前記繰り返し回数が基準回数になっても、前記収束判定手段において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うようにさせたことを特徴とする請求項６又は７に記載の流体構造連成解析プログラム。
【請求項１】
コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析方法であって、
前記コンピュータが、
運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出段階と、
流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出段階と、
前記圧力算出段階で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定段階と、
前記誤差が所定範囲になった場合、前記算出した圧力を用いて前記構造物の変位を算出する変位算出段階と
を含むことを特徴とする流体構造連成解析方法。
【請求項２】
前記コンピュータは、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算段階を更に含むことを特徴とする請求項１に記載の流体構造連成解析方法。
【請求項３】
前記コンピュータは、
有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断段階を更に含み、
気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うことを特徴とする請求項２に記載の流体構造連成解析方法。
【請求項４】
前記コンピュータは、
前記圧力算出段階において収束判定に至るまでの繰り返し回数を計数し、
前記繰り返し回数が基準回数になっても、前記収束判定段階において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うことを特徴とする請求項２又は３に記載の流体構造連成解析方法。
【請求項５】
コンピュータを用いて、構造解析と流体解析とを連成させて、構造物に接する流体シミュレーションを算出するための流体構造連成解析プログラムであって、
前記コンピュータを、
運動量保存式から前記構造物の変位を算出した場合に、この算出された変位に基づいて、構造物と流体の境界流速を算出する境界流速算出手段と、
流体の気泡成長計算から算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて、気泡の質量に関する項を含む質量保存式と、運動量保存式とから、密度及び圧力を算出する圧力算出手段、
前記圧力算出手段で算出された密度と、状態方程式から算出される密度との誤差を算出し、この誤差が所定範囲になるまで圧力算出を繰り返す収束判定手段、及び
前記誤差が所定範囲になった場合、前記算出した圧力を用いて前記構造物の変位を算出する変位算出手段
として機能させることを特徴とする流体構造連成解析プログラム。
【請求項６】
前記コンピュータを、前記質量保存式及び前記運動量保存式から流体に関する圧力が算出された場合には、その圧力を用いて前記気泡成長計算を行う気泡成長計算手段として更に機能させることを特徴とする請求項５に記載の流体構造連成解析プログラム。
【請求項７】
前記コンピュータを、
有限要素として分割された流体のメッシュにおいて気泡の気泡成長速度が一定値よりも大きいか否かを判断する判断手段として更に機能させ、
気泡の気泡成長速度が一定値よりも大きいメッシュについてのみ前記気泡成長計算を行うようにさせたことを特徴とする請求項６に記載の流体構造連成解析プログラム。
【請求項８】
前記コンピュータを、
前記圧力算出手段において収束判定に至るまでの繰り返し回数を計数し、
前記繰り返し回数が基準回数になっても、前記収束判定手段において算出される誤差が前記所定範囲にならなかった場合、気泡成長速度を算出するための時間間隔を小さくして算出した気泡半径、気泡の気泡成長速度及びボイド率を用いて前記圧力算出を行うようにさせたことを特徴とする請求項６又は７に記載の流体構造連成解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２００６−７２５６６（Ｐ２００６−７２５６６Ａ）
【公開日】平成１８年３月１６日（２００６．３．１６）
【国際特許分類】
【出願番号】特願２００４−２５３４５１（Ｐ２００４−２５３４５１）
【出願日】平成１６年８月３１日（２００４．８．３１）
【出願人】（５９２１３１９０６）みずほ情報総研株式会社 (187)
【出願人】（０００００４０９７）日本原子力研究所 (55)
【Ｆターム（参考）】
【公開日】平成１８年３月１６日（２００６．３．１６）
【国際特許分類】
【出願日】平成１６年８月３１日（２００４．８．３１）
【出願人】（５９２１３１９０６）みずほ情報総研株式会社 (187)
【出願人】（０００００４０９７）日本原子力研究所 (55)
【Ｆターム（参考）】
[ Back to top ]