数値流体計算方法、数値流体計算装置、プログラム、及び、記録媒体
【課題】ボリュームＣＡＤシステムと一体化した数値流体計算装置を実現する。
【解決手段】数値流体計算装置１は、Ｎ≦Ｎmaxとなるように対象物Ωを６面体の内部セル及び（６＋Ｎ）面体の境界セルに分割するモデリング部１１と、境界セルの三角切断面ｋｔの面積Ａ_ktの面積を算出する境界セル情報生成部１３と、解くべき代数方程式に含まれる、対流流束ρＵφの三角切断面ｋｔ上での積分値、及び、拡散流束Γ∇φの各三角切断面上での積分値を、三角切断面の面積Ａ_ktと境界条件φ_kt及び（∇φ）_ktとから算出する境界セル係数設定部１５とを備える。
【解決手段】数値流体計算装置１は、Ｎ≦Ｎmaxとなるように対象物Ωを６面体の内部セル及び（６＋Ｎ）面体の境界セルに分割するモデリング部１１と、境界セルの三角切断面ｋｔの面積Ａ_ktの面積を算出する境界セル情報生成部１３と、解くべき代数方程式に含まれる、対流流束ρＵφの三角切断面ｋｔ上での積分値、及び、拡散流束Γ∇φの各三角切断面上での積分値を、三角切断面の面積Ａ_ktと境界条件φ_kt及び（∇φ）_ktとから算出する境界セル係数設定部１５とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、有限体積法を用いて対流拡散方程式を解く数値流体計算方法、及び、数値流体計算装置に関する。また、そのような数値流体計算装置としてコンピュータを動作させるためのプログラム、及び、そのようなプログラムが記録された記録媒体に関する。
【背景技術】
【０００２】
コンピュータの進歩により、ものつくりの技術は飛躍的な発展を遂げた。例えば、設計を支援するＣＡＤ（Computer Aided Design）システム、解析を支援するＣＡＥ（Computer Aided Engineering）システム、加工を支援するＣＡＭ（Computer Aided Manufacturing）システム、検査を支援するＣＡＴ（Computer Aided Testing）システムは、ものつくりを行う上でもはや不可欠なツールとなっている。
【０００３】
しかしながら、これらのシステムは、別々のソフトウェアにより実現されており、データ構造にも汎用性がない。例えば、ＣＡＤシステムにおいて広く用いられているワイヤフレームモデルは、物体の表面形状を表現するものであり、物体の内部構造を表現するものではない。したがって、ワイヤフレームモデルをシミュレーションにそのまま使用することはできない。
【０００４】
設計にもシミュレーションにも利用可能な汎用モデル化手法として、ボクセルモデルが知られている。ボクセルモデルは、物体をボクセルと呼ばれる立方体の集合として表現するモデル化手法である。各ボクセルに物理量を担わせることによって、ボクセルモデルを設計にもシミュレーションにも利用することができる。ただし、物体をボクセルの集合として表現した場合、その物体の表面の情報（例えば法線の方向）が欠落してしまう。したがって、ボクセルモデルを用いた数値流体計算では境界値問題を精度良く解くことができない。
【０００５】
ボクセルモデルの改良として、物体の表面と交わるボクセルを切断立方体（「Kitta Cube」と呼ばれる）に置き換えた切断立方体モデルが知られている。切断立方体とは、各辺上に高々１つの切断点が設けられた立方体を、３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体のことである。切断点は、基礎となる立方体の各辺と、物体の表面との交点に設けられる。切断立方体モデルは、物体の表面を良く近似する、切断三角形からなる表面を有するので境界値問題に適している。
【０００６】
しかしながら、切断立方体における三角切断面のパターンは１０64通りに達し、切断点の配置パターンは１４４通りに及ぶ。したがって、切断立方体モデルを用いた数値流体計算アルゴリズムは多数の場合分けを含み、これをコーディングすることは現実的ではない。また、仮にコーディングできたとしても計算量が厖大になり、要求される処理速度を達成することができない。
【０００７】
このような問題を解決するための試みが、非特許文献１に記載されている。具体的には、切断立方体モデルにおける切断立方体をカットセルに置き換える（例えば１辺を共有する２つの三角切断面を１つの切断面に置き換える）ことにより、３次元数値流体計算が可能になるとの予想が示されている。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】雷、外１名、「ＶＣＡＤにより非圧縮性粘性流体のシミュレーション」、第１９回数値流体力学シンポジウム講演要旨集、２００５年１２月、ｐ．５１、ＣＤ−ＲＯＭ ＮＯ．Ａ１−３
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、非特許文献１に記載されているように切断立方体モデルにおける切断立方体をカットセルに置き換えた場合、全てのカットセルの切断面を隣接するカットセルの切断面と矛盾なく整合させることができない。すなわち、各カットセルの切断面により構成されるカットセルモデルの表面に穴があいてしまい、解の精度が大幅に低下するという問題を生じる。
【００１０】
本発明は、上記の問題に鑑みてなされたものであり、有限体積法を用いて対流拡散方程式を解く数値流体計算方法において、切断立方体をカットセル等に置き換えることなく、必要な場合分けの数を低減させた数値流体計算方法を実現することにある。
【課題を解決するための手段】
【００１１】
上記課題を解決するために、本発明に係る数値流体方法は、対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を、コンピュータによって解く数値流体計算方法において、計算対象領域を複数のセルに分割するセル分割ステップであって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割ステップと、上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出ステップと、上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出ステップにて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定ステップとを含んでいる、ことを特徴としている。
【００１２】
また、上記課題を解決するために、本発明に係る数値流体計算装置は、対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を解く数値流体計算装置において、計算対象領域を複数のセルに分割するセル分割手段であって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割手段と、上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出手段と、上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出手段にて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定手段と、を備えている、ことを特徴としている。
【００１３】
上記の構成によれば、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体が用いられる。したがって、計算対象領域の境界を、これらの三角切断面によって精度良く近似することができる。そして、解くべき代数方程式に含まれる対流項及び拡散項を、これらの三角切断面を用いて算出するので、解くべき代数方程式を精度良く算出することができる。したがって、上記物理量の各セル上での値を、精度良く算出することができる。
【００１４】
しかも、各境界セルにおける切断面数（三角切断面の数）Ｎは、Ｎmax以下に制限されている。したがって、解くべき代数方程式（の係数）を算出するために必要な場合分けを、各境界セルにおける切断面数に制限を加えない場合よりも少なくすることができる。
【００１５】
本発明に係る数値流体計算方法は、上記境界セルの各々について、該境界セルに隣接する隣接セルとの間の境界面の面積を算出する境界面積算出ステップと、上記境界セルの各々について、上記代数方程式に含まれる係数であって、該境界セルに隣接する隣接セル上の上記物理量に係る係数を、上記境界面積算出ステップにて算出された各境界面の面積を用いて算出する係数算出ステップとを更に含んでいる、ことが好ましい。
【００１６】
上記の構成によれば、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体が用いられる。そして、解くべき代数方程式に含まれる係数を、これらの境界セルの境界面の面積を用いて算出するので、解くべき代数方程式を精度良く算出することができる。したがって、上記物理量の各セル上での値を精度良く算出することができる。
【００１７】
本発明に係る数値流体計算方法において、上記該境界セルの各々に含まれる三角切断面数Ｎの上限値Ｎmaxが４である、ことが好ましい。
【００１８】
上記の構成によれば、境界セルの境界面の面積を算出するために必要な場合分けを、著しく少なくすることができる（図５〜図８参照）。
【００１９】
なお、本発明に係る数値流体計算装置は、コンピュータによって実現してもよい。この場合、コンピュータを上記各手段として動作させることにより、本発明に係る数値流体計算装置をコンピュータにおいて実現するプログラム、及び、そのプログラムを記録したコンピュータ読み取り可能な記録媒体も、本発明の範疇に入る。
【発明の効果】
【００２０】
以上のように、本発明によれば、有限体積法を用いて対流拡散方程式を解く数値流体計算方法において、切断立方体をカットセル等に置き換えることなく、必要な場合分けの数を低減させた数値流体計算方法を実現することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明の実施形態を示すものであり、数値流体計算装置の構成を示すブロック図である。
【図２】本発明の実施形態を示すものであり、コンピュータを用いて実現された数値流体計算装置のハードウェア構成を示すブロック図である。
【図３】本発明の実施形態を示すものであり、切断立方体モデルの構成を示す説明図である。（ａ）は、対象物の断面図であり、（ｂ）は、内部セルの斜視図であり、（ｃ）は、境界セルの斜視図である。
【図４】本発明の実施形態を示すものであり、境界セルの分類を示す説明図である。
【図５】本発明の実施形態を示すものであり、切断点をもたない境界面の構成を示す平面図である。
【図６】本発明の実施形態を示すものであり、２つの切断点を含む境界面の構成を示す平面図である。
【図７】本発明の実施形態を示すものであり、３つの切断点を含む境界面の構成を示す平面図である。
【図８】本発明の実施形態を示すものであり、４つの切断点を含む境界面の構成を示す平面図である。
【図９】本発明の実施例を示すものであり、計算対象とする系の構造を示す模式図である。
【図１０】図９に示す系における、平均Nusseltの評価結果を示すグラフである。（ａ）は、小球面上での平均Nusselt数を示し、（ｂ）は、大球面上での平均Nusselt数を示す。
【発明を実施するための形態】
【００２２】
〔本発明の基礎となる事項〕
本発明の実施形態について説明する前に、本発明の基礎となる事項について説明する。なお、本明細書においては、下付き文字を前に_を付して表し、上付き文字を前に^を付して表す。例えば、ＡeをＡ_eと表し、Ｌ2をＬ^2と表す。
【００２３】
（切断立方体モデル）
本発明の数値流体計算方法に用いる切断立方体（Kitta Cube)モデルについて図３〜図４を参照して説明する。
【００２４】
切断立方体モデルは、発明者らが開発したボリュームＣＡＤシステムにおけるモデリング手法であり、図３（ａ）に示すように、対象物Ωの内部に含まれる６面体の内部セルＩＣと、対象物Ωの境界面∂Ωを含む６＋Ｎ面体の境界セルＢＣとを用いて対象物Ωを表現するモデリング手法である。内部セルは、「非境界セル」と呼ばれることもある。
【００２５】
切断立方体モデルにおいては、内部セルＩＣとして、図３（ｂ）に示すような６面体（例えば立方体）が用いられる。内部セルＩＣにおいては、６つの境界面が、それぞれ、内部セルＩＣに隣接する６つの隣接セル（内部セルの場合もあれば、境界セルの場合もある）との境界を成す。また、内部セルＩＣは、１つの格子点を含む。内部セルＩＣに含まれる格子点の位置は、例えば、内部セルＩＣの重心位置である。
【００２６】
本明細書においては、注目するセル（以下「注目セル」と呼称する）に含まれる格子点をＰと記し、注目セルの６つの境界面をｅ（ｘ軸正方向），ｗ（ｘ軸負方向），ｎ（ｙ軸正方向），ｓ（ｙ軸負方向），ｔ（ｚ軸正方向），ｂ（ｚ軸負方向）と記す。また、境界面ｅ，ｗ，ｓ，ｎ，ｔ，ｂを介して注目セルに隣接するセル（以下「隣接セル」と呼称する）に含まれる格子点を、それぞれ、Ｅ，Ｗ，Ｓ，Ｎ，Ｔ，Ｂと記す。また、境界面ｅ，ｗ，ｓ，ｎ，ｔ，ｂの面積を、それぞれ、Ａ_e，Ａ_w，Ａ_n，Ａ_s，Ａ_t，Ａ_bと記す（図３（ｂ）参照）。
【００２７】
また、切断立方体モデルにおいては、境界セルＢＣとして、図３（ｃ）に示すような（６＋Ｎ）面体が用いられる。すなわち、各辺上に高々１つの切断点が設けられた基礎となる６面体（例えば内部セルと合同な立方体）を、３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体が用いられる。切断点は、基礎となる６面体の各辺と、対象物Ωの境界面∂Ωとの交点に設けられる（各セルのサイズを十分に小さくとれば、各辺上に設けられる切断点の数を１以下にすることができる）。図３（ｃ）に例示した境界セルＢＣは、切断点Ｋ_1，Ｋ_2，Ｋ_3を頂点とする三角切断面ｋｔ_１と、切断点Ｋ_1，Ｋ_2，Ｋ_４を頂点とする三角切断面ｋｔ_２で切断して得られたものである。
【００２８】
境界セルＢＣにおいては、三角切断面ＫＴを除く６つの境界面が、それぞれ、境界セルＢＣに隣接する６つの隣接セル（内部セルの場合もあれば、境界セルの場合もある）との境界を成す。一方、三角切断面ＫＴは、対象物Ωと外界（固体）との間の境界を構成する。また、境界セルＢＣは、内部セルＩＣと同様、１つの格子点を含む。格子点の位置は、例えば、境界セルＢＣの重心位置である。
【００２９】
なお、三角切断面の数Ｎは境界セル毎に異なり得ることに留意されたい。ただし、本発明においては、境界セルの切断面数Ｎに上限を設ける。換言すれば、切断面数Ｎが予め定められた上限値を上回る境界セルが生じないように、セルサイズを十分に小さく設定する。本実施形態においては、境界セルの切断面数Ｎの上限を４とする。換言すれば、切断面数Ｎが４を上回る境界セルが生じないように、セルサイズを十分に小さく設定する。これにより、各境界セルは、図４に示す４つのタイプのうちの何れかに分類される。
【００３０】
（対流拡散方程式）
本発明が対象とする方程式は、対象物Ω上で定義された物理量φ＝φ（ｔ，ｘ，y，ｚ）が従う対流拡散方程式（１）である。対流拡散方程式（１）は、「移流拡散方程式」、あるいは、「輸送方程式」と呼ばれることもある。物理量φとしては、単位質量あたりの運動量の各成分、あるいは、単位質量あたりのエネルギーが想定されるが、これに限定されるものではない。
【００３１】
【数１】
【００３２】
ここで、Ｕ＝（ｕ，ｖ，ｗ）は流速（ベクトル場）、ρは密度（スカラ場）、Γは拡散係数（スカラ場）、Ｓは生成項（スカラ場）である。ρＵφは、対流流束と呼ばれ、Γ∇φは、拡散流束と呼ばれる。
【００３３】
なお、以下の説明では、簡単のために、密度ρ及び拡散係数Γが時間及び位置に依存せずに一定であるものと仮定するが、本発明はこれに限定されるものではない。すなわち、密度ρ及び拡散係数Γは、時間又は位置に依存するものであってもよい。また、以下の説明では、簡単のために、生成項Ｓが物理量φに比例する項Ｓ_P×φと定数項Ｓ_Cとの和で表せるものと仮定するが、本発明はこれに限定されるものではない。すなわち、生成項Ｓが物理量φの２次以上の項を含んでいてもよい。
【００３４】
本発明の数値流体計算方法においては、有限体積法を用いて対流拡散方程式（１）を解く。すなわち、対流拡散方程式（１）を各セル上で積分することによって得られる代数方程式からなる方程式系（連立方程式）をコンピュータによって解く。各代数方程式は、注目セルに含まれる格子点Ｐにおける物理量φ_Pと、隣接セルに含まれる格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝における物理量φ_NBとの間に成り立つ代数方程式である。これらの代数方程式は、対流拡散方程式（１）から以下のようにして導出することができる。
【００３５】
すなわち、まず、対流拡散方程式（１）を各内部セル上で積分することにより式（２）を得る。また、対流拡散方程式（１）を各境界セル上で積分することにより式（３）を得る。
【００３６】
【数２】
【００３７】
【数３】
【００３８】
ここで、Δｔは、時間ステップであり、ΔＶは、注目セルの体積である。また、Ａ_nbは、注目セルの各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の面積であり、Ａ_ktは、注目セルの各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の面積である。なお、式（２）及び式（３）においては、面積Ａ_nbを、境界面ｎｂの外向法線方向を向いた面積ベクトルとして扱っている。また、式（３）においては、面積Ａ_KTを、三角切断面ｋｔの外向法線方向を向いた面積ベクトルとして面積Ａ_KTを扱っている。
【００３９】
また、（ρφ）_Pは、注目セルに含まれる格子点Ｐにおけるρφの値である。（ρＵφ）_nb、及び、（Γ∇φ）_nbは、それぞれ、注目セルの各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝における対流流束ρＵφ、及び、拡散流束Γ∇φの値である。和Σ_nbは、注目セルの６つの境界面の全てを渡る。
【００４０】
また、（ρＵφ）_kt、及び、（Γ∇φ）_ktは、注目セルの各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝（の代表点）における対流流束ρＵφ、及び、拡散流束Γ∇φの値である。和Σ_ktは、注目セルのＮ個の三角切断面の全てを渡る。
【００４１】
次に、適当な差分法を用いて、物理量φの境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝における値φ_nbを、物理量φの格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝における値φ_NBに置き換えることにより代数方程式（４）を得る。
【００４２】
【数４】
【００４３】
・ハイブリッド法
差分法としてハイブリッド法を用いた場合、内部セルに対する代数方程式（４）における係数ａ_NBは、式（５ａ）〜式（５ｈ）により与えられる。また、内部セルに対する代数方程式（４）における定数ｂは、式（５ｉ）により与えられる。
【００４４】
【数５】
【００４５】
ここで、括弧||・||は、括弧内の要素のうちで最大の要素を指す。また、Ｆ_nbは、式（６）により定義される量であり、Ｄ_nbは、式（７）により定義される量である。式（７）において、ｒ_（P→NB）は、注目セルの格子点Ｐを始点とし、隣接セルの格子点ＮＢを終点とするベクトルを表す。ｅ_（P→NB）は、ｒ_（P→NB）と同じ向きをむいた単位ベクトルである。
【００４６】
【数６】
【００４７】
【数７】
【００４８】
一方、差分法としてハイブリッド法を用いた場合、境界セルに対する代数方程式（４）における係数ａ_NBは、内部セルに対する代数方程式（４）と同様、式（５ａ）〜式（５ｈ）により与えられる。また、境界セルに対する代数方程式（４）における定数ｂは、対流流束ρＵφ及び拡散流速Γ∇φの三角切断面上での積分を含む式（８）により与えられる。
【００４９】
【数８】
【００５０】
対流流束ρＵφ及び拡散流速Γ∇φの三角切断面上での積分は、境界条件として与えられたφ及び∇φの境界面∂Ω上での値から算出することができる。
【００５１】
・風上差分法
差分法として風上差分法を用いた場合、内部セルに対する代数方程式（４）における係数ａ_NBは、式（９ａ）〜式（９ｂ）によって与えられる。また、内部セルに対する代数方程式（４）における定数ｂは、式（９ｃ）によって与えられる。Ｆ_nbは、上述した式（６）により定義される量であり、Ｄ_nbは、式（７）により定義される量である。
【００５２】
【数９】
【００５３】
ここで、風上差分法とは、物理量φの境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝における値φ_nbを、式（１０）によって物理量φの格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝における値φ_NBに置き換えることを指す。式（１０）において、Δｒ_（P→NB）は、注目セルの格子点Ｐを始点とし、隣接セルの格子点ＮＢを終点とするベクトルを表し、Δｒ_（NB→P）は、隣接セルの格子点ＮＢを始点とし、注目セルの格子点Ｐを終点とするベクトルを表す。また、Ｍ_nb^+及びＭ_nb^-は式（１１）により定義される。
【００５４】
【数１０】
【００５５】
【数１１】
【００５６】
（１）式における時間項の注目セル上での積分を、式（１２）のように離散化した場合、式（９ｂ）における係数ａ_t、及び、式（９ｃ）における定数ｂ_temporalは、式（１３）のように与えられる。
【００５７】
【数１２】
【００５８】
【数１３】
【００５９】
また、式（９ｃ）における定数ｂ_CT及び定数ｂ_DTは、式（１４ａ）〜式（１４ｂ）ように与えられる。
【００６０】
【数１４】
【００６１】
一方、差分法として風上差分法を用いた場合、境界セルに対する代数方程式（４）における係数ａ_NBは、内部セルに対する代数方程式（４）と同様、式（９ａ）〜式（９ｂ）によって与えられる。また、境界セルに対する代数方程式（４）における定数ｂは、対流流束ρＵφ及び拡散流速Γ∇φの三角切断面上での積分を含む式（１５）によって与えられる。
【００６２】
【数１５】
【００６３】
なお、ここではハイブリッド法、及び、風上差分法から導かれる代数方程式（４）を用いた有限体積法について説明したが、本発明の適用範囲はこれに限定されるものではない。すなわち、中心差分法、べき乗法、指数法（厳密法）など、他の離散化手法から導かれる代数方程式（４）を用いた有限体積法も本発明の適用範囲に含まれる。
【００６４】
また、ここでは流速Ｕを既知として、流速Ｕにより規定される流れ場における物理量φの輸送について説明したが、本発明はこれに限定されるものではない。すなわち、代数方程式（４）と圧力修正方程式とを連立させることにより、非定常非圧縮粘性流体の流速Ｕ*を算出する数値流体計算方法も本発明の範疇に含まれる。
【００６５】
〔実施形態〕
以下、本発明の一実施形態について、図面に基づいて説明する。
【００６６】
（数値流体計算装置の構成）
本実施形態に係る数値流体計算装置１の構成について、図１を参照して説明する。図１は、数値流体計算装置１の構成を示すブロック図である。数値流体計算装置１は、モデリング部１１と、切断面数カウント部１２（同図において「カウント部」と略記）と、境界セル情報生成部１３と、内部セル情報生成部１４と、境界セル係数設定部１５と、内部セル係数設定部１６と、行列演算部１７と、収束判定部１８と、ＵＩ部１９と、パラメータ記憶部２１と、メッシュデータ記憶部２２とを備えている。
【００６７】
モデリング部１１は、対象物Ωを複数のセルに分割するための手段である。具体的には、切断立方体モデルに従って対象物Ωを複数のセルに分割し、各セルの頂点座標を表すメッシュデータを生成する。上述したとおり、モデリング部１１は、内部セルとして６面体を用い、境界セルとして（６＋Ｎ）面体を用いる。モデリング部１１によって生成されたメッシュデータは、切断面カウント部１２に提供される。
【００６８】
切断面数カウント部１２は、モデリング部１１から提供されたメッシュデータを参照して各境界セルにおける切断面数Ｎをカウントし、最大切断面数が４以下であるか否かを判定する。
【００６９】
最大切断面数が４以下である場合、切断面数カウント部１２は、モデリング部１１に「ＯＫ」を返す。切断面数カウント部１２から「ＯＫ」が返されると、モデリング部１１は、各セルの頂点座標を表すメッシュデータをメッシュデータ記憶部２２に格納する。
【００７０】
一方、最大切断面数が４よりも大きい場合、切断面数カウント部１２は、モデリング部１１に「ＮＧ」を返す。切断面数カウント部１２から「ＮＧ」が返されると、モデリング部１１は、対象物Ωをより小さいサイズの内部セルと境界セルとに再分割する。モデリング部１１は、最大切断面数が４以下になるまで（切断面数カウント部１２から「ＯＫ」が返されるまで）再分割を繰り返し、最終的には最大切断面数が４以下となるメッシュデータをメッシュデータ記憶部２２に格納する。
【００７１】
境界セル情報生成部１３は、メッシュデータ記憶部２２に格納されたメッシュデータから、各境界セルについて、（１）体積ΔＶ、（２）格子点Ｐの座標、（３）格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝と格子点Ｐとの距離｜ＮＢ−Ｐ｜、（４）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の代表点の座標、（５）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の面積Ａ_nb、（６）各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の代表点の座標、（７）各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の面積Ａ_kt、（８）各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の法線ベクトルｎ_ktを算出する。境界セル情報生成部１３により算出された諸量は、境界セル情報として境界セル係数設定部１５に提供される。
【００７２】
なお、境界セルの格子点Ｐの座標は、例えば、その境界セルの頂点座標の和を頂点数で除算することによって容易に算出することができる。また、各境界面ｎｂの代表点の座標は、例えば、その境界面ｎｂの頂点座標の和を頂点数で除算することによって容易に算出することができる。また、各境界面ｎｂの面積Ａ_nbは、境界面の構成毎に導出した求積公式（後述）を用いて算出することができる。境界セルの体積ΔＶについても同様である。
【００７３】
また、各三角切断面ｋｔの代表点の座標は、例えば、その三角切断面ｋｔの頂点座標の和を頂点数３で割ることによって容易に算出することができる。また、各三角切断面ｋｔの面積Ａ_ktは、例えば、その三角切断面ｋｔの３つの頂点（切断点）の位置ベクトルをα，β，γとして、Ａ_kt＝｜（β−α）×（γ−α）｜／２により算出することができる（「×」は外積）。また、三角切断面ｋｔの法線ベクトルｎ_ktは、例えば、n_kt＝（β−α）×（γ−α）／｜（β−α）×（γ−α）｜により算出することができる（「×」は外積）。
【００７４】
内部セル情報生成部１４は、メッシュデータ記憶部２２に格納されたメッシュデータから、各内部セルについて、（１）体積ΔＶ，（２）格子点Ｐの座標、（３）各格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝と格子点Ｐとの距離｜ＮＢーＰ｜、（４）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の代表点の座標、（５）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の面積Ａ_nbを算出する。内部セル情報生成部１４により算出された諸量は、内部セル情報として内部セル係数設定部１６に提供される。
【００７５】
なお、内部セルの格子点Ｐの座標は、例えば、その内部セルの頂点座標の和を頂点数８で除算することによって容易に算出することができる。また、各境界面ｎｂの代表点の座標は、例えば、その境界面ｎｂの頂点座標の和を頂点数４で除算することによって容易に算出することができる。また、各境界面ｎｂの面積Ａ_nbは、例えば、互いに隣接する２辺の長さを乗算することによって容易に算出することができる。例えば、内部セルが各辺の長さがＬの立方体である場合、Ａｎｂ＝Ｌ^2（一定）である。また、内部セルの体積ΔＶは、例えば、互いに隣接する３辺の長さを乗算することによって、容易に算出することができる。例えば、内部セルが各辺の長さがＬの立方体である場合、ΔＶ＝Ｌ^3（一定）である。
【００７６】
代数方程式（４）からなる方程式系を構成するために必要なパラメータ{Δｔ，ρ，Γ，Ｓ，Ｕ＝（ｕ，ｖ，ｗ），φ_kt，（∇φ）_kt｝の値は、ＵＩ部１９を介してユーザにより設され、パラメータ記憶部２１に格納される。境界セル係数設定部１５、及び、内部セル係数設定部１６は、パラメータ記憶部２１に記憶されたこれらのパラメータを自由に参照することができる。
【００７７】
境界セル係数設定部１５は、対流流束ρＵφの各三角切断面ｋｔ上での積分の和Σ（ρＵφ）_kt・Ａ_kt、及び、拡散流速Γ∇φの各三角切断面ｋｔ上での積分の和を、パラメータ記憶部２１に格納された境界値φ_kt及び（Δφ）_ktなどを参照して算出する。
【００７８】
具体的には、（１）密度ρと流速Ｕ_ktと境界条件φ_ktとをパラメータ記憶部２１から読み出し、（２）法線ベクトルｎ_ktと面積Ａ_ktとを境界セル情報生成部１３から取得し、（３）流速Ｕ_ktと法線ベクトルｎ_ktとの内積Ｕ_kt・ｎ_ktを算出し、（４）内積Ｕ_kt・ｎ_ktに密度ρ、境界条件φ_kt、及び、面積Ａ_ktを乗じることによって、（ρＵφ）_kt・Ａ_ktを算出する。そして、得られた（ρＵφ）_kt・Ａ_ktを合算することによりΣ（ρＵφ）_kt・Ａ_ktを得る。
【００７９】
同様に、（１）拡散係数Γと境界条件（∇φ）_ktとをパラメータ記憶部２１から読み出し、（２）法線ベクトルｎ_ktと面積Ａ_ktとを境界セル情報生成部１３から取得し、（３）境界条件（∇φ）_ktと法線ベクトルｎ_ktとの内積（∇φ）_kt・ｎ_ktを算出し、（４）内積（∇φ）_kt・ｎ_ktに拡散係数Γ及び面積Ａ_ktを乗じることによって、（Γ∇φ）_kt・Ａ_ktを算出する。そして、得られた（Γ∇φ）_kt・Ａ_ktを合算することによって、Σ（Γ∇φ）_kt・Ａ_ktを得る。
【００８０】
境界セル係数設定部１５は、更に、パラメータ記憶部２１から読み出した｛ρ，Ｕ，Γ，Ｓ_P｝と、境界セル情報生成部１３から取得した{Ａ_nb，｜ＮＢーＰ｜，ΔＶ｝とを参照し、式（５ａ）〜（５ｈ）に従って係数ａ_P，ａ_E，ａ_W，ａ_N，ａ_S，ａ_T，ａ_Bを算出する。また、先に算出した｛Σ（Γ∇φ）_kt・Ａ_kt，Σ（ρＵφ）_kt・Ａ_kt｝と、パラメータ記憶部２１から読み出した｛Δｔ，Ｓ_C，ρ，φ_P^T｝と、境界セル情報生成部１３から取得したΔＶとを参照し、式（８）に従って定数ｂを算出する。
【００８１】
内部セル係数設定部１６は、パラメータ記憶部２１から読み出した｛ρ，Ｕ，Γ，Ｓ_P｝と、境界セル情報生成部１３から取得した{Ａ_nb，｜ＮＢーＰ｜，ΔＶ｝とを参照し、式（５ａ）〜（５ｈ）に従って係数ａ_P，ａ_E，ａ_W，ａ_N，ａ_S，ａ_T，ａ_Bを算出する。また、パラメータ記憶部２１から読み出した｛Δｔ，Ｓ_C，ρ，φ_P^T｝と、境界セル情報生成部１３から取得したΔＶとを参照し、式（５ｉ）に従って定数ｂを算出する。
【００８２】
なお、定数ｂを算出するために境界セル係数設定部１５及び内部セル係数設定部１６が参照するφ_P^tは、ＵＩ部１９を介してユーザによって設定されたものではなく、行列演算部１７によって前回算出されたφ_Pである。ただし、初回の演算において参照するφ_P^tは、ＵＩ部１９を介してユーザによって設定された初期値φ_P^0である。
【００８３】
行列演算部１７は、境界セル係数設定部１５及び内部セル係数設定部１６によって設定された係数ａ_ＮＢ及び定数ｂにより規定される、代数方程式（４）を解くための手段である。代数方程式（４）を解くためのアルゴリズムとしては、ヤコビの反復法、ガウス・ザイデルの反復法、逐次過緩和法（ＳＯＲ：Successive Over Relaxation）などが知られているが、ここでは、収束の速い逐次過緩和法を用いる。
【００８４】
行列演算部１７は、代数方程式（４）の解（各格子点における物理量φの値）を更新する度に、更新後の解を収束判定部１８に提供する。収束判定部１８は、更新前の解と更新後の解との差（のノルム）が予め定められた閾値以下になるか否かを判定する。
【００８５】
更新前の解と更新後の解との差が予め定められた閾値以下である場合、収束判定部１８は、行列演算部１７に「ＯＫ」を返す。収束判定部１８から「ＯＫ」が返されると、行列演算部１７は、更新後の解を外部に出力すると共に、更新後の解をパラメータ記憶部２１に格納する。パラメータ記憶部２１に格納された解は、時刻ｔにおける解φ_P^tとして、時刻ｔ＋Δｔにおける物理量φ_P^t+Δtを算出するために参照される。
【００８６】
更新前の解と更新後の解との差が予め定められた閾値よりも大きい場合、収束判定部１８は、行列演算部１７に「ＮＧ」を返す。収束判定部１８から「ＮＧ」が返されると、行列演算部１７は、再度、解を更新する。行列演算部１７は、更新前の解と更新後の解との差が予め定められた閾値以下になるまで（収束判定部１８から「ＯＫ」が返されるまで）解の更新を繰り返し、最終的には更新前の解との差が予め定められた閾値以下となる収束解を得る。
【００８７】
（コンピュータを用いた構成例）
数値流体計算装置１は、例えば、コンピュータ（電子計算機）を用いて実現することができる。図２は、コンピュータを用いて実現された数値流体計算装置１のハードウェア構成を例示したブロック図である。
【００８８】
数値流体計算装置１は、図２に示したように、バス１１０を介して互いに接続された演算装置１２０と、主記憶装置１３０と、補助記憶装置１４０と、入出力インタフェース１５０とを備えている。演算装置１２０として利用可能なデバイスとしては、ＣＰＵ（Central Processing Unit）を挙げることができる。また、主記憶装置１３０として利用可能なデバイスとしては、例えば、半導体ＲＡＭ（random access memory）を挙げることができる。また、補助記憶装置１４０として利用可能なデバイスとしては、例えば、ハードディスクドライブを挙げることができる。
【００８９】
入出力インタフェース１５０には、図２に示したように、入力装置２００及び出力装置３００が接続される。入力装置２００は、パラメータ{Δｔ，ρ，Γ，Ｓ，Ｕ＝（ｕ，ｖ，ｗ），φ_kt，（∇φ）_kt}の値をユーザが数値流体計算装置１に入力するための手段であり、例えば、キーボードである。入力装置２００を介して入力されたパラメータ{Δｔ，ρ，Γ，Ｓ，Ｕ＝（ｕ，ｖ，ｗ），φ_kt，（∇φ）_kt｝は、演算装置１２０がこれを参照することができるよう主記憶装置１３０に格納される。すなわち、主記憶装置１３０は、パラメータ記憶手段２１として利用される。
【００９０】
一方、出力装置３００は、代数方程式（４）の解（各格子点における物理量φの値）を出力するための手段であり、例えば、モニタである。なお、代数方程式（４）の解を出力装置３００を介して出力する代わりに、代数方程式（４）の解を補助記憶装置１４０に格納するようにしてもよい。
【００９１】
補助記憶装置１４０には、コンピュータを数値流体計算装置１として動作させるための数値流体計算プログラムが格納されている。数値流体計算プログラムは、モデリングプログラムと、切断面数カウントプログラムと、境界セル情報設定プログラムと、境界セル係数設定プログラムと、内部セル係数設定プログラムと、行列演算プログラムと、収束判定プログラムと、ＵＩプログラムとを含んでいる。また、補助記憶装置１４０は、メッシュデータを格納するメッシュデータ記憶部２２としても利用される。
【００９２】
演算装置１２０は、主記憶装置１３０上に展開されたモデリングプログラムに含まれる命令を実行することによって、コンピュータをモデリング部１１として機能させる。同様に、主記憶装置１３０上に展開された切断面数カウントプログラム、境界セル情報設定プログラム、境界セル係数設定プログラム、内部セル係数設定プログラム、行列演算プログラム、収束判定プログラム、及び、ＵＩプログラムに含まれる命令を実行することによって、コンピュータ１を切断面数カウント部１２、境界セル情報生成部１３、内部セル情報生成部１４、境界セル係数設定部１５、内部セル係数設定部１６、行列演算部１７、及び、収束判定部１８として機能させる。
【００９３】
本発明の目的は、上記各プログラムのプログラムコード（実行形式プログラム、中間コードプログラム、ソースプログラム）をコンピュータで読み取り可能に記録した記録媒体を数値流体計算装置１に供給し、数値流体計算装置１がその記録媒体に記録されているプログラムコードを読み出し実行することによっても達成可能である。
【００９４】
上記記録媒体としては、例えば、磁気テープやカセットテープ等のテープ系、フロッピー（登録商標）ディスク／ハードディスク等の磁気ディスクやＣＤ−ＲＯＭ／ＭＯ／ＭＤ／ＤＶＤ／ＣＤ−Ｒ等の光ディスクを含むディスク系、ＩＣカード（メモリカードを含む）／光カード等のカード系、あるいはマスクＲＯＭ／ＥＰＲＯＭ／ＥＥＰＲＯＭ／フラッシュＲＯＭ等の半導体メモリ系などを用いることができる。
【００９５】
また、数値流体計算装置１を通信ネットワークと接続可能に構成し、上記プログラムコードを通信ネットワークを介して数値流体計算装置１に供給するようにしてもよい。この通信ネットワークとしては、とくに限定されず、例えば、インターネット、イントラネット、エキストラネット、ＬＡＮ、ＩＳＤＮ、ＶＡＮ、ＣＡＴＶ通信網、仮想専用網（virtual private network）、電話回線網、移動体通信網、衛星通信網等が利用可能である。また、通信ネットワークを構成する伝送媒体としては、とくに限定されず、例えば、ＩＥＥＥ１３９４、ＵＳＢ、電力線搬送、ケーブルＴＶ回線、電話線、ＡＤＳＬ回線等の有線でも、ＩｒＤＡやリモコンのような赤外線、Ｂｌｕｅｔｏｏｔｈ（登録商標）、８０２．１１無線、ＨＤＲ、携帯電話網、衛星回線、地上波デジタル網等の無線でも利用可能である。なお、本発明は、上記プログラムコードが電子的な伝送で具現化された、搬送波に埋め込まれたコンピュータデータ信号の形態でも実現され得る。
【００９６】
（求積公式）
境界セルにおける境界面ｎｂの面積Ａ_nbは、上述したように、境界面の構成毎に導出された求積公式を用いて算出することができる。以下、この求積公式について、図５〜図８を参照して説明する。なお、以下の説明においては、基礎となる立方体の面に占める境界面ｎｂの割合を開口率Ａ_fluidと記す。
【００９７】
「基礎となる立方体の各辺に高々１つの切断点を設ける」という条件から、考慮すべき組み合わせは、（１）切断点がない場合、（２）切断点が２つある場合、（３）切断点が３つある場合、（４）切断点が４つある場合の何れかに分類される。以下（１）〜（４）の場合について順に説明する。
【００９８】
（１）切断点がない場合
境界面ｎｂ上に切断点がない場合、境界面ｎｂ全体が対象物Ωの内部（流体）に含まれるか（図５（ａ）参照）、又は、境界面ｎｂ全体が対象物Ωの外部（固体）に含まれるか（図５（ｂ）参照）の何れかである。前者の場合、開口率Ａ_fluid＝１であり、後者の場合、開口率Ａ_fluid＝０である。
【００９９】
（２）切断点が２つある場合
境界面ｎｂが２つの切断面を含む場合、境界面ｎｂの構成は、図６（ａ）〜図６（ｃ）に示す３通りのうちの何れかになる。この場合、開口率Ａ_fluidは、式（１６）により算出される。
【０１００】
【数１６】
【０１０１】
（３）切断点が３つある場合
境界面ｎｂ上に切断点が３つある場合、境界面ｎｂの構成は、図７（ａ）〜図７（ｅ）に示す５通りのうちの何れかになる。図７（ａ）〜図７（ｄ）に示す構成の場合、開口率Ａ_fluidは、式（１７）により算出される。
【０１０２】
【数１７】
【０１０３】
一方、図７（ｅ）の構成は、問題が２次元に退化してしまうため望ましくない。このため、図７（ｅ）に示す切断点の配置が生じた場合には、メッシュを切り直す（対象物Ωを再分割する）ことが望ましい。
【０１０４】
（４）切断点が４つある場合
境界面ｎｂ上に切断点が４つある場合、境界面ｎｂの構成は、図８（ａ）〜図８（ｄ）に示す４通りのうちの何れかになる。図８（ａ）〜図８（ｂ）に示す構成の場合、開口率Ａ_fluidは、式（１８）により算出される。
【０１０５】
【数１８】
【０１０６】
なお、図８（ｃ）〜図８（ｄ）に示す構成は、問題が２次元に退化してしまうため望ましくない。このため、図８（ｃ）〜図８（ｄ）に示す切断点の配置が生じた場合には、メッシュを切り直す（対象物Ωを再分割する）ことが望ましい。
【０１０７】
〔実施例〕
最後に、温度Ｔが従う対流拡散方程式に本発明を適用した実施例について説明する。
【０１０８】
本実施例において計算対象とする系を図９に示す。本実施例において計算対象とする系は、図９に示すように、共通の中心を持つ小球面（半径ｄ_i）と大球面（半径ｄ_o）との間の空間である（ｄ_i＜ｄ_o）。境界条件としては、小球面（内側境界）の温度＝Ｔ_i、大球面（外側境界）の温度＝Ｔ_oを課す（Ｔ_i＞Ｔ_o）。
【０１０９】
レーリー数Ｒａを下表の各値に設定したうえで、本発明の数値流体計算方法に従って各セルにおける温度Ｔの値を算出した。なお、レーリー数Ｒａとは、式（１９）により定義される量である。
【０１１０】
【数１９】
【０１１１】
【表１】
【０１１２】
計算精度をチェックするために、各球面上での平均Nusselt数（以下「Ｎｕ数」」と略記する）を評価した。なお、小球面上での平均Nusselt数は、式（２０ａ）のように定義される量であり、大球面上での平均Nusselt数は、式（２０ｂ）のように定義される量である。
【０１１３】
【数２０】
【０１１４】
評価結果を図１０に示す。図１０（ａ）は、小球面上での平均Nusselt数を示し、図１０（ｂ）は、大球面上での平均Nusselt数を示している。同図においては、下記参考文献に記載された平均Nusselt数を比較例として示している。本発明の数値流体計算方法によれば、表１に示した何れのレーリー数に対しても、精度の高い数値流体計算が実現されていることが分かる。
【０１１５】
〔参考文献〕Greg V K, "Natural convection between concentric spheres", International Journal of Heat and Mass Transfer, 35（8), 1991, pp1935-1945
【産業上の利用可能性】
【０１１６】
本発明は数値流体計算に広く利用することができる。特に、（Ｖ）ＣＡＤシステムと連携する数値流体計算に好適に利用することができる。
【符号の説明】
【０１１７】
１ 数値流体計算装置
１１ モデリング部
１２ 切断面数カウント部
１３ 境界セル情報生成部
１４ 内部セル情報生成部
１５ 境界セル係数設定部
１６ 内部セル係数設定部
１７ 行列演算部
１８ 収束判定部
１９ ＵＩ部
２１ パラメータ記憶部
２２ メッシュデータ記憶部
ＩＣ 内部セル
ＢＣ 境界セル
ｅ，ｗ，ｎ，ｓ，ｔ，ｂ 境界面
ｋｔ 三角切断面
Ｐ 注目セルに含まれる格子点
Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ 隣接セルに含まれる格子点
【技術分野】
【０００１】
本発明は、有限体積法を用いて対流拡散方程式を解く数値流体計算方法、及び、数値流体計算装置に関する。また、そのような数値流体計算装置としてコンピュータを動作させるためのプログラム、及び、そのようなプログラムが記録された記録媒体に関する。
【背景技術】
【０００２】
コンピュータの進歩により、ものつくりの技術は飛躍的な発展を遂げた。例えば、設計を支援するＣＡＤ（Computer Aided Design）システム、解析を支援するＣＡＥ（Computer Aided Engineering）システム、加工を支援するＣＡＭ（Computer Aided Manufacturing）システム、検査を支援するＣＡＴ（Computer Aided Testing）システムは、ものつくりを行う上でもはや不可欠なツールとなっている。
【０００３】
しかしながら、これらのシステムは、別々のソフトウェアにより実現されており、データ構造にも汎用性がない。例えば、ＣＡＤシステムにおいて広く用いられているワイヤフレームモデルは、物体の表面形状を表現するものであり、物体の内部構造を表現するものではない。したがって、ワイヤフレームモデルをシミュレーションにそのまま使用することはできない。
【０００４】
設計にもシミュレーションにも利用可能な汎用モデル化手法として、ボクセルモデルが知られている。ボクセルモデルは、物体をボクセルと呼ばれる立方体の集合として表現するモデル化手法である。各ボクセルに物理量を担わせることによって、ボクセルモデルを設計にもシミュレーションにも利用することができる。ただし、物体をボクセルの集合として表現した場合、その物体の表面の情報（例えば法線の方向）が欠落してしまう。したがって、ボクセルモデルを用いた数値流体計算では境界値問題を精度良く解くことができない。
【０００５】
ボクセルモデルの改良として、物体の表面と交わるボクセルを切断立方体（「Kitta Cube」と呼ばれる）に置き換えた切断立方体モデルが知られている。切断立方体とは、各辺上に高々１つの切断点が設けられた立方体を、３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体のことである。切断点は、基礎となる立方体の各辺と、物体の表面との交点に設けられる。切断立方体モデルは、物体の表面を良く近似する、切断三角形からなる表面を有するので境界値問題に適している。
【０００６】
しかしながら、切断立方体における三角切断面のパターンは１０64通りに達し、切断点の配置パターンは１４４通りに及ぶ。したがって、切断立方体モデルを用いた数値流体計算アルゴリズムは多数の場合分けを含み、これをコーディングすることは現実的ではない。また、仮にコーディングできたとしても計算量が厖大になり、要求される処理速度を達成することができない。
【０００７】
このような問題を解決するための試みが、非特許文献１に記載されている。具体的には、切断立方体モデルにおける切断立方体をカットセルに置き換える（例えば１辺を共有する２つの三角切断面を１つの切断面に置き換える）ことにより、３次元数値流体計算が可能になるとの予想が示されている。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】雷、外１名、「ＶＣＡＤにより非圧縮性粘性流体のシミュレーション」、第１９回数値流体力学シンポジウム講演要旨集、２００５年１２月、ｐ．５１、ＣＤ−ＲＯＭ ＮＯ．Ａ１−３
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、非特許文献１に記載されているように切断立方体モデルにおける切断立方体をカットセルに置き換えた場合、全てのカットセルの切断面を隣接するカットセルの切断面と矛盾なく整合させることができない。すなわち、各カットセルの切断面により構成されるカットセルモデルの表面に穴があいてしまい、解の精度が大幅に低下するという問題を生じる。
【００１０】
本発明は、上記の問題に鑑みてなされたものであり、有限体積法を用いて対流拡散方程式を解く数値流体計算方法において、切断立方体をカットセル等に置き換えることなく、必要な場合分けの数を低減させた数値流体計算方法を実現することにある。
【課題を解決するための手段】
【００１１】
上記課題を解決するために、本発明に係る数値流体方法は、対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を、コンピュータによって解く数値流体計算方法において、計算対象領域を複数のセルに分割するセル分割ステップであって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割ステップと、上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出ステップと、上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出ステップにて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定ステップとを含んでいる、ことを特徴としている。
【００１２】
また、上記課題を解決するために、本発明に係る数値流体計算装置は、対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を解く数値流体計算装置において、計算対象領域を複数のセルに分割するセル分割手段であって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割手段と、上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出手段と、上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出手段にて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定手段と、を備えている、ことを特徴としている。
【００１３】
上記の構成によれば、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体が用いられる。したがって、計算対象領域の境界を、これらの三角切断面によって精度良く近似することができる。そして、解くべき代数方程式に含まれる対流項及び拡散項を、これらの三角切断面を用いて算出するので、解くべき代数方程式を精度良く算出することができる。したがって、上記物理量の各セル上での値を、精度良く算出することができる。
【００１４】
しかも、各境界セルにおける切断面数（三角切断面の数）Ｎは、Ｎmax以下に制限されている。したがって、解くべき代数方程式（の係数）を算出するために必要な場合分けを、各境界セルにおける切断面数に制限を加えない場合よりも少なくすることができる。
【００１５】
本発明に係る数値流体計算方法は、上記境界セルの各々について、該境界セルに隣接する隣接セルとの間の境界面の面積を算出する境界面積算出ステップと、上記境界セルの各々について、上記代数方程式に含まれる係数であって、該境界セルに隣接する隣接セル上の上記物理量に係る係数を、上記境界面積算出ステップにて算出された各境界面の面積を用いて算出する係数算出ステップとを更に含んでいる、ことが好ましい。
【００１６】
上記の構成によれば、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体が用いられる。そして、解くべき代数方程式に含まれる係数を、これらの境界セルの境界面の面積を用いて算出するので、解くべき代数方程式を精度良く算出することができる。したがって、上記物理量の各セル上での値を精度良く算出することができる。
【００１７】
本発明に係る数値流体計算方法において、上記該境界セルの各々に含まれる三角切断面数Ｎの上限値Ｎmaxが４である、ことが好ましい。
【００１８】
上記の構成によれば、境界セルの境界面の面積を算出するために必要な場合分けを、著しく少なくすることができる（図５〜図８参照）。
【００１９】
なお、本発明に係る数値流体計算装置は、コンピュータによって実現してもよい。この場合、コンピュータを上記各手段として動作させることにより、本発明に係る数値流体計算装置をコンピュータにおいて実現するプログラム、及び、そのプログラムを記録したコンピュータ読み取り可能な記録媒体も、本発明の範疇に入る。
【発明の効果】
【００２０】
以上のように、本発明によれば、有限体積法を用いて対流拡散方程式を解く数値流体計算方法において、切断立方体をカットセル等に置き換えることなく、必要な場合分けの数を低減させた数値流体計算方法を実現することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明の実施形態を示すものであり、数値流体計算装置の構成を示すブロック図である。
【図２】本発明の実施形態を示すものであり、コンピュータを用いて実現された数値流体計算装置のハードウェア構成を示すブロック図である。
【図３】本発明の実施形態を示すものであり、切断立方体モデルの構成を示す説明図である。（ａ）は、対象物の断面図であり、（ｂ）は、内部セルの斜視図であり、（ｃ）は、境界セルの斜視図である。
【図４】本発明の実施形態を示すものであり、境界セルの分類を示す説明図である。
【図５】本発明の実施形態を示すものであり、切断点をもたない境界面の構成を示す平面図である。
【図６】本発明の実施形態を示すものであり、２つの切断点を含む境界面の構成を示す平面図である。
【図７】本発明の実施形態を示すものであり、３つの切断点を含む境界面の構成を示す平面図である。
【図８】本発明の実施形態を示すものであり、４つの切断点を含む境界面の構成を示す平面図である。
【図９】本発明の実施例を示すものであり、計算対象とする系の構造を示す模式図である。
【図１０】図９に示す系における、平均Nusseltの評価結果を示すグラフである。（ａ）は、小球面上での平均Nusselt数を示し、（ｂ）は、大球面上での平均Nusselt数を示す。
【発明を実施するための形態】
【００２２】
〔本発明の基礎となる事項〕
本発明の実施形態について説明する前に、本発明の基礎となる事項について説明する。なお、本明細書においては、下付き文字を前に_を付して表し、上付き文字を前に^を付して表す。例えば、ＡeをＡ_eと表し、Ｌ2をＬ^2と表す。
【００２３】
（切断立方体モデル）
本発明の数値流体計算方法に用いる切断立方体（Kitta Cube)モデルについて図３〜図４を参照して説明する。
【００２４】
切断立方体モデルは、発明者らが開発したボリュームＣＡＤシステムにおけるモデリング手法であり、図３（ａ）に示すように、対象物Ωの内部に含まれる６面体の内部セルＩＣと、対象物Ωの境界面∂Ωを含む６＋Ｎ面体の境界セルＢＣとを用いて対象物Ωを表現するモデリング手法である。内部セルは、「非境界セル」と呼ばれることもある。
【００２５】
切断立方体モデルにおいては、内部セルＩＣとして、図３（ｂ）に示すような６面体（例えば立方体）が用いられる。内部セルＩＣにおいては、６つの境界面が、それぞれ、内部セルＩＣに隣接する６つの隣接セル（内部セルの場合もあれば、境界セルの場合もある）との境界を成す。また、内部セルＩＣは、１つの格子点を含む。内部セルＩＣに含まれる格子点の位置は、例えば、内部セルＩＣの重心位置である。
【００２６】
本明細書においては、注目するセル（以下「注目セル」と呼称する）に含まれる格子点をＰと記し、注目セルの６つの境界面をｅ（ｘ軸正方向），ｗ（ｘ軸負方向），ｎ（ｙ軸正方向），ｓ（ｙ軸負方向），ｔ（ｚ軸正方向），ｂ（ｚ軸負方向）と記す。また、境界面ｅ，ｗ，ｓ，ｎ，ｔ，ｂを介して注目セルに隣接するセル（以下「隣接セル」と呼称する）に含まれる格子点を、それぞれ、Ｅ，Ｗ，Ｓ，Ｎ，Ｔ，Ｂと記す。また、境界面ｅ，ｗ，ｓ，ｎ，ｔ，ｂの面積を、それぞれ、Ａ_e，Ａ_w，Ａ_n，Ａ_s，Ａ_t，Ａ_bと記す（図３（ｂ）参照）。
【００２７】
また、切断立方体モデルにおいては、境界セルＢＣとして、図３（ｃ）に示すような（６＋Ｎ）面体が用いられる。すなわち、各辺上に高々１つの切断点が設けられた基礎となる６面体（例えば内部セルと合同な立方体）を、３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体が用いられる。切断点は、基礎となる６面体の各辺と、対象物Ωの境界面∂Ωとの交点に設けられる（各セルのサイズを十分に小さくとれば、各辺上に設けられる切断点の数を１以下にすることができる）。図３（ｃ）に例示した境界セルＢＣは、切断点Ｋ_1，Ｋ_2，Ｋ_3を頂点とする三角切断面ｋｔ_１と、切断点Ｋ_1，Ｋ_2，Ｋ_４を頂点とする三角切断面ｋｔ_２で切断して得られたものである。
【００２８】
境界セルＢＣにおいては、三角切断面ＫＴを除く６つの境界面が、それぞれ、境界セルＢＣに隣接する６つの隣接セル（内部セルの場合もあれば、境界セルの場合もある）との境界を成す。一方、三角切断面ＫＴは、対象物Ωと外界（固体）との間の境界を構成する。また、境界セルＢＣは、内部セルＩＣと同様、１つの格子点を含む。格子点の位置は、例えば、境界セルＢＣの重心位置である。
【００２９】
なお、三角切断面の数Ｎは境界セル毎に異なり得ることに留意されたい。ただし、本発明においては、境界セルの切断面数Ｎに上限を設ける。換言すれば、切断面数Ｎが予め定められた上限値を上回る境界セルが生じないように、セルサイズを十分に小さく設定する。本実施形態においては、境界セルの切断面数Ｎの上限を４とする。換言すれば、切断面数Ｎが４を上回る境界セルが生じないように、セルサイズを十分に小さく設定する。これにより、各境界セルは、図４に示す４つのタイプのうちの何れかに分類される。
【００３０】
（対流拡散方程式）
本発明が対象とする方程式は、対象物Ω上で定義された物理量φ＝φ（ｔ，ｘ，y，ｚ）が従う対流拡散方程式（１）である。対流拡散方程式（１）は、「移流拡散方程式」、あるいは、「輸送方程式」と呼ばれることもある。物理量φとしては、単位質量あたりの運動量の各成分、あるいは、単位質量あたりのエネルギーが想定されるが、これに限定されるものではない。
【００３１】
【数１】
【００３２】
ここで、Ｕ＝（ｕ，ｖ，ｗ）は流速（ベクトル場）、ρは密度（スカラ場）、Γは拡散係数（スカラ場）、Ｓは生成項（スカラ場）である。ρＵφは、対流流束と呼ばれ、Γ∇φは、拡散流束と呼ばれる。
【００３３】
なお、以下の説明では、簡単のために、密度ρ及び拡散係数Γが時間及び位置に依存せずに一定であるものと仮定するが、本発明はこれに限定されるものではない。すなわち、密度ρ及び拡散係数Γは、時間又は位置に依存するものであってもよい。また、以下の説明では、簡単のために、生成項Ｓが物理量φに比例する項Ｓ_P×φと定数項Ｓ_Cとの和で表せるものと仮定するが、本発明はこれに限定されるものではない。すなわち、生成項Ｓが物理量φの２次以上の項を含んでいてもよい。
【００３４】
本発明の数値流体計算方法においては、有限体積法を用いて対流拡散方程式（１）を解く。すなわち、対流拡散方程式（１）を各セル上で積分することによって得られる代数方程式からなる方程式系（連立方程式）をコンピュータによって解く。各代数方程式は、注目セルに含まれる格子点Ｐにおける物理量φ_Pと、隣接セルに含まれる格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝における物理量φ_NBとの間に成り立つ代数方程式である。これらの代数方程式は、対流拡散方程式（１）から以下のようにして導出することができる。
【００３５】
すなわち、まず、対流拡散方程式（１）を各内部セル上で積分することにより式（２）を得る。また、対流拡散方程式（１）を各境界セル上で積分することにより式（３）を得る。
【００３６】
【数２】
【００３７】
【数３】
【００３８】
ここで、Δｔは、時間ステップであり、ΔＶは、注目セルの体積である。また、Ａ_nbは、注目セルの各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の面積であり、Ａ_ktは、注目セルの各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の面積である。なお、式（２）及び式（３）においては、面積Ａ_nbを、境界面ｎｂの外向法線方向を向いた面積ベクトルとして扱っている。また、式（３）においては、面積Ａ_KTを、三角切断面ｋｔの外向法線方向を向いた面積ベクトルとして面積Ａ_KTを扱っている。
【００３９】
また、（ρφ）_Pは、注目セルに含まれる格子点Ｐにおけるρφの値である。（ρＵφ）_nb、及び、（Γ∇φ）_nbは、それぞれ、注目セルの各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝における対流流束ρＵφ、及び、拡散流束Γ∇φの値である。和Σ_nbは、注目セルの６つの境界面の全てを渡る。
【００４０】
また、（ρＵφ）_kt、及び、（Γ∇φ）_ktは、注目セルの各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝（の代表点）における対流流束ρＵφ、及び、拡散流束Γ∇φの値である。和Σ_ktは、注目セルのＮ個の三角切断面の全てを渡る。
【００４１】
次に、適当な差分法を用いて、物理量φの境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝における値φ_nbを、物理量φの格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝における値φ_NBに置き換えることにより代数方程式（４）を得る。
【００４２】
【数４】
【００４３】
・ハイブリッド法
差分法としてハイブリッド法を用いた場合、内部セルに対する代数方程式（４）における係数ａ_NBは、式（５ａ）〜式（５ｈ）により与えられる。また、内部セルに対する代数方程式（４）における定数ｂは、式（５ｉ）により与えられる。
【００４４】
【数５】
【００４５】
ここで、括弧||・||は、括弧内の要素のうちで最大の要素を指す。また、Ｆ_nbは、式（６）により定義される量であり、Ｄ_nbは、式（７）により定義される量である。式（７）において、ｒ_（P→NB）は、注目セルの格子点Ｐを始点とし、隣接セルの格子点ＮＢを終点とするベクトルを表す。ｅ_（P→NB）は、ｒ_（P→NB）と同じ向きをむいた単位ベクトルである。
【００４６】
【数６】
【００４７】
【数７】
【００４８】
一方、差分法としてハイブリッド法を用いた場合、境界セルに対する代数方程式（４）における係数ａ_NBは、内部セルに対する代数方程式（４）と同様、式（５ａ）〜式（５ｈ）により与えられる。また、境界セルに対する代数方程式（４）における定数ｂは、対流流束ρＵφ及び拡散流速Γ∇φの三角切断面上での積分を含む式（８）により与えられる。
【００４９】
【数８】
【００５０】
対流流束ρＵφ及び拡散流速Γ∇φの三角切断面上での積分は、境界条件として与えられたφ及び∇φの境界面∂Ω上での値から算出することができる。
【００５１】
・風上差分法
差分法として風上差分法を用いた場合、内部セルに対する代数方程式（４）における係数ａ_NBは、式（９ａ）〜式（９ｂ）によって与えられる。また、内部セルに対する代数方程式（４）における定数ｂは、式（９ｃ）によって与えられる。Ｆ_nbは、上述した式（６）により定義される量であり、Ｄ_nbは、式（７）により定義される量である。
【００５２】
【数９】
【００５３】
ここで、風上差分法とは、物理量φの境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝における値φ_nbを、式（１０）によって物理量φの格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝における値φ_NBに置き換えることを指す。式（１０）において、Δｒ_（P→NB）は、注目セルの格子点Ｐを始点とし、隣接セルの格子点ＮＢを終点とするベクトルを表し、Δｒ_（NB→P）は、隣接セルの格子点ＮＢを始点とし、注目セルの格子点Ｐを終点とするベクトルを表す。また、Ｍ_nb^+及びＭ_nb^-は式（１１）により定義される。
【００５４】
【数１０】
【００５５】
【数１１】
【００５６】
（１）式における時間項の注目セル上での積分を、式（１２）のように離散化した場合、式（９ｂ）における係数ａ_t、及び、式（９ｃ）における定数ｂ_temporalは、式（１３）のように与えられる。
【００５７】
【数１２】
【００５８】
【数１３】
【００５９】
また、式（９ｃ）における定数ｂ_CT及び定数ｂ_DTは、式（１４ａ）〜式（１４ｂ）ように与えられる。
【００６０】
【数１４】
【００６１】
一方、差分法として風上差分法を用いた場合、境界セルに対する代数方程式（４）における係数ａ_NBは、内部セルに対する代数方程式（４）と同様、式（９ａ）〜式（９ｂ）によって与えられる。また、境界セルに対する代数方程式（４）における定数ｂは、対流流束ρＵφ及び拡散流速Γ∇φの三角切断面上での積分を含む式（１５）によって与えられる。
【００６２】
【数１５】
【００６３】
なお、ここではハイブリッド法、及び、風上差分法から導かれる代数方程式（４）を用いた有限体積法について説明したが、本発明の適用範囲はこれに限定されるものではない。すなわち、中心差分法、べき乗法、指数法（厳密法）など、他の離散化手法から導かれる代数方程式（４）を用いた有限体積法も本発明の適用範囲に含まれる。
【００６４】
また、ここでは流速Ｕを既知として、流速Ｕにより規定される流れ場における物理量φの輸送について説明したが、本発明はこれに限定されるものではない。すなわち、代数方程式（４）と圧力修正方程式とを連立させることにより、非定常非圧縮粘性流体の流速Ｕ*を算出する数値流体計算方法も本発明の範疇に含まれる。
【００６５】
〔実施形態〕
以下、本発明の一実施形態について、図面に基づいて説明する。
【００６６】
（数値流体計算装置の構成）
本実施形態に係る数値流体計算装置１の構成について、図１を参照して説明する。図１は、数値流体計算装置１の構成を示すブロック図である。数値流体計算装置１は、モデリング部１１と、切断面数カウント部１２（同図において「カウント部」と略記）と、境界セル情報生成部１３と、内部セル情報生成部１４と、境界セル係数設定部１５と、内部セル係数設定部１６と、行列演算部１７と、収束判定部１８と、ＵＩ部１９と、パラメータ記憶部２１と、メッシュデータ記憶部２２とを備えている。
【００６７】
モデリング部１１は、対象物Ωを複数のセルに分割するための手段である。具体的には、切断立方体モデルに従って対象物Ωを複数のセルに分割し、各セルの頂点座標を表すメッシュデータを生成する。上述したとおり、モデリング部１１は、内部セルとして６面体を用い、境界セルとして（６＋Ｎ）面体を用いる。モデリング部１１によって生成されたメッシュデータは、切断面カウント部１２に提供される。
【００６８】
切断面数カウント部１２は、モデリング部１１から提供されたメッシュデータを参照して各境界セルにおける切断面数Ｎをカウントし、最大切断面数が４以下であるか否かを判定する。
【００６９】
最大切断面数が４以下である場合、切断面数カウント部１２は、モデリング部１１に「ＯＫ」を返す。切断面数カウント部１２から「ＯＫ」が返されると、モデリング部１１は、各セルの頂点座標を表すメッシュデータをメッシュデータ記憶部２２に格納する。
【００７０】
一方、最大切断面数が４よりも大きい場合、切断面数カウント部１２は、モデリング部１１に「ＮＧ」を返す。切断面数カウント部１２から「ＮＧ」が返されると、モデリング部１１は、対象物Ωをより小さいサイズの内部セルと境界セルとに再分割する。モデリング部１１は、最大切断面数が４以下になるまで（切断面数カウント部１２から「ＯＫ」が返されるまで）再分割を繰り返し、最終的には最大切断面数が４以下となるメッシュデータをメッシュデータ記憶部２２に格納する。
【００７１】
境界セル情報生成部１３は、メッシュデータ記憶部２２に格納されたメッシュデータから、各境界セルについて、（１）体積ΔＶ、（２）格子点Ｐの座標、（３）格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝と格子点Ｐとの距離｜ＮＢ−Ｐ｜、（４）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の代表点の座標、（５）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の面積Ａ_nb、（６）各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の代表点の座標、（７）各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の面積Ａ_kt、（８）各三角切断面ｋｔ∈｛ｋｔ_1，ｋｔ_2，…，ｋｔ_N｝の法線ベクトルｎ_ktを算出する。境界セル情報生成部１３により算出された諸量は、境界セル情報として境界セル係数設定部１５に提供される。
【００７２】
なお、境界セルの格子点Ｐの座標は、例えば、その境界セルの頂点座標の和を頂点数で除算することによって容易に算出することができる。また、各境界面ｎｂの代表点の座標は、例えば、その境界面ｎｂの頂点座標の和を頂点数で除算することによって容易に算出することができる。また、各境界面ｎｂの面積Ａ_nbは、境界面の構成毎に導出した求積公式（後述）を用いて算出することができる。境界セルの体積ΔＶについても同様である。
【００７３】
また、各三角切断面ｋｔの代表点の座標は、例えば、その三角切断面ｋｔの頂点座標の和を頂点数３で割ることによって容易に算出することができる。また、各三角切断面ｋｔの面積Ａ_ktは、例えば、その三角切断面ｋｔの３つの頂点（切断点）の位置ベクトルをα，β，γとして、Ａ_kt＝｜（β−α）×（γ−α）｜／２により算出することができる（「×」は外積）。また、三角切断面ｋｔの法線ベクトルｎ_ktは、例えば、n_kt＝（β−α）×（γ−α）／｜（β−α）×（γ−α）｜により算出することができる（「×」は外積）。
【００７４】
内部セル情報生成部１４は、メッシュデータ記憶部２２に格納されたメッシュデータから、各内部セルについて、（１）体積ΔＶ，（２）格子点Ｐの座標、（３）各格子点ＮＢ∈｛Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ｝と格子点Ｐとの距離｜ＮＢーＰ｜、（４）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の代表点の座標、（５）各境界面ｎｂ∈｛ｅ，ｗ，ｎ，ｓ，ｔ，ｂ｝の面積Ａ_nbを算出する。内部セル情報生成部１４により算出された諸量は、内部セル情報として内部セル係数設定部１６に提供される。
【００７５】
なお、内部セルの格子点Ｐの座標は、例えば、その内部セルの頂点座標の和を頂点数８で除算することによって容易に算出することができる。また、各境界面ｎｂの代表点の座標は、例えば、その境界面ｎｂの頂点座標の和を頂点数４で除算することによって容易に算出することができる。また、各境界面ｎｂの面積Ａ_nbは、例えば、互いに隣接する２辺の長さを乗算することによって容易に算出することができる。例えば、内部セルが各辺の長さがＬの立方体である場合、Ａｎｂ＝Ｌ^2（一定）である。また、内部セルの体積ΔＶは、例えば、互いに隣接する３辺の長さを乗算することによって、容易に算出することができる。例えば、内部セルが各辺の長さがＬの立方体である場合、ΔＶ＝Ｌ^3（一定）である。
【００７６】
代数方程式（４）からなる方程式系を構成するために必要なパラメータ{Δｔ，ρ，Γ，Ｓ，Ｕ＝（ｕ，ｖ，ｗ），φ_kt，（∇φ）_kt｝の値は、ＵＩ部１９を介してユーザにより設され、パラメータ記憶部２１に格納される。境界セル係数設定部１５、及び、内部セル係数設定部１６は、パラメータ記憶部２１に記憶されたこれらのパラメータを自由に参照することができる。
【００７７】
境界セル係数設定部１５は、対流流束ρＵφの各三角切断面ｋｔ上での積分の和Σ（ρＵφ）_kt・Ａ_kt、及び、拡散流速Γ∇φの各三角切断面ｋｔ上での積分の和を、パラメータ記憶部２１に格納された境界値φ_kt及び（Δφ）_ktなどを参照して算出する。
【００７８】
具体的には、（１）密度ρと流速Ｕ_ktと境界条件φ_ktとをパラメータ記憶部２１から読み出し、（２）法線ベクトルｎ_ktと面積Ａ_ktとを境界セル情報生成部１３から取得し、（３）流速Ｕ_ktと法線ベクトルｎ_ktとの内積Ｕ_kt・ｎ_ktを算出し、（４）内積Ｕ_kt・ｎ_ktに密度ρ、境界条件φ_kt、及び、面積Ａ_ktを乗じることによって、（ρＵφ）_kt・Ａ_ktを算出する。そして、得られた（ρＵφ）_kt・Ａ_ktを合算することによりΣ（ρＵφ）_kt・Ａ_ktを得る。
【００７９】
同様に、（１）拡散係数Γと境界条件（∇φ）_ktとをパラメータ記憶部２１から読み出し、（２）法線ベクトルｎ_ktと面積Ａ_ktとを境界セル情報生成部１３から取得し、（３）境界条件（∇φ）_ktと法線ベクトルｎ_ktとの内積（∇φ）_kt・ｎ_ktを算出し、（４）内積（∇φ）_kt・ｎ_ktに拡散係数Γ及び面積Ａ_ktを乗じることによって、（Γ∇φ）_kt・Ａ_ktを算出する。そして、得られた（Γ∇φ）_kt・Ａ_ktを合算することによって、Σ（Γ∇φ）_kt・Ａ_ktを得る。
【００８０】
境界セル係数設定部１５は、更に、パラメータ記憶部２１から読み出した｛ρ，Ｕ，Γ，Ｓ_P｝と、境界セル情報生成部１３から取得した{Ａ_nb，｜ＮＢーＰ｜，ΔＶ｝とを参照し、式（５ａ）〜（５ｈ）に従って係数ａ_P，ａ_E，ａ_W，ａ_N，ａ_S，ａ_T，ａ_Bを算出する。また、先に算出した｛Σ（Γ∇φ）_kt・Ａ_kt，Σ（ρＵφ）_kt・Ａ_kt｝と、パラメータ記憶部２１から読み出した｛Δｔ，Ｓ_C，ρ，φ_P^T｝と、境界セル情報生成部１３から取得したΔＶとを参照し、式（８）に従って定数ｂを算出する。
【００８１】
内部セル係数設定部１６は、パラメータ記憶部２１から読み出した｛ρ，Ｕ，Γ，Ｓ_P｝と、境界セル情報生成部１３から取得した{Ａ_nb，｜ＮＢーＰ｜，ΔＶ｝とを参照し、式（５ａ）〜（５ｈ）に従って係数ａ_P，ａ_E，ａ_W，ａ_N，ａ_S，ａ_T，ａ_Bを算出する。また、パラメータ記憶部２１から読み出した｛Δｔ，Ｓ_C，ρ，φ_P^T｝と、境界セル情報生成部１３から取得したΔＶとを参照し、式（５ｉ）に従って定数ｂを算出する。
【００８２】
なお、定数ｂを算出するために境界セル係数設定部１５及び内部セル係数設定部１６が参照するφ_P^tは、ＵＩ部１９を介してユーザによって設定されたものではなく、行列演算部１７によって前回算出されたφ_Pである。ただし、初回の演算において参照するφ_P^tは、ＵＩ部１９を介してユーザによって設定された初期値φ_P^0である。
【００８３】
行列演算部１７は、境界セル係数設定部１５及び内部セル係数設定部１６によって設定された係数ａ_ＮＢ及び定数ｂにより規定される、代数方程式（４）を解くための手段である。代数方程式（４）を解くためのアルゴリズムとしては、ヤコビの反復法、ガウス・ザイデルの反復法、逐次過緩和法（ＳＯＲ：Successive Over Relaxation）などが知られているが、ここでは、収束の速い逐次過緩和法を用いる。
【００８４】
行列演算部１７は、代数方程式（４）の解（各格子点における物理量φの値）を更新する度に、更新後の解を収束判定部１８に提供する。収束判定部１８は、更新前の解と更新後の解との差（のノルム）が予め定められた閾値以下になるか否かを判定する。
【００８５】
更新前の解と更新後の解との差が予め定められた閾値以下である場合、収束判定部１８は、行列演算部１７に「ＯＫ」を返す。収束判定部１８から「ＯＫ」が返されると、行列演算部１７は、更新後の解を外部に出力すると共に、更新後の解をパラメータ記憶部２１に格納する。パラメータ記憶部２１に格納された解は、時刻ｔにおける解φ_P^tとして、時刻ｔ＋Δｔにおける物理量φ_P^t+Δtを算出するために参照される。
【００８６】
更新前の解と更新後の解との差が予め定められた閾値よりも大きい場合、収束判定部１８は、行列演算部１７に「ＮＧ」を返す。収束判定部１８から「ＮＧ」が返されると、行列演算部１７は、再度、解を更新する。行列演算部１７は、更新前の解と更新後の解との差が予め定められた閾値以下になるまで（収束判定部１８から「ＯＫ」が返されるまで）解の更新を繰り返し、最終的には更新前の解との差が予め定められた閾値以下となる収束解を得る。
【００８７】
（コンピュータを用いた構成例）
数値流体計算装置１は、例えば、コンピュータ（電子計算機）を用いて実現することができる。図２は、コンピュータを用いて実現された数値流体計算装置１のハードウェア構成を例示したブロック図である。
【００８８】
数値流体計算装置１は、図２に示したように、バス１１０を介して互いに接続された演算装置１２０と、主記憶装置１３０と、補助記憶装置１４０と、入出力インタフェース１５０とを備えている。演算装置１２０として利用可能なデバイスとしては、ＣＰＵ（Central Processing Unit）を挙げることができる。また、主記憶装置１３０として利用可能なデバイスとしては、例えば、半導体ＲＡＭ（random access memory）を挙げることができる。また、補助記憶装置１４０として利用可能なデバイスとしては、例えば、ハードディスクドライブを挙げることができる。
【００８９】
入出力インタフェース１５０には、図２に示したように、入力装置２００及び出力装置３００が接続される。入力装置２００は、パラメータ{Δｔ，ρ，Γ，Ｓ，Ｕ＝（ｕ，ｖ，ｗ），φ_kt，（∇φ）_kt}の値をユーザが数値流体計算装置１に入力するための手段であり、例えば、キーボードである。入力装置２００を介して入力されたパラメータ{Δｔ，ρ，Γ，Ｓ，Ｕ＝（ｕ，ｖ，ｗ），φ_kt，（∇φ）_kt｝は、演算装置１２０がこれを参照することができるよう主記憶装置１３０に格納される。すなわち、主記憶装置１３０は、パラメータ記憶手段２１として利用される。
【００９０】
一方、出力装置３００は、代数方程式（４）の解（各格子点における物理量φの値）を出力するための手段であり、例えば、モニタである。なお、代数方程式（４）の解を出力装置３００を介して出力する代わりに、代数方程式（４）の解を補助記憶装置１４０に格納するようにしてもよい。
【００９１】
補助記憶装置１４０には、コンピュータを数値流体計算装置１として動作させるための数値流体計算プログラムが格納されている。数値流体計算プログラムは、モデリングプログラムと、切断面数カウントプログラムと、境界セル情報設定プログラムと、境界セル係数設定プログラムと、内部セル係数設定プログラムと、行列演算プログラムと、収束判定プログラムと、ＵＩプログラムとを含んでいる。また、補助記憶装置１４０は、メッシュデータを格納するメッシュデータ記憶部２２としても利用される。
【００９２】
演算装置１２０は、主記憶装置１３０上に展開されたモデリングプログラムに含まれる命令を実行することによって、コンピュータをモデリング部１１として機能させる。同様に、主記憶装置１３０上に展開された切断面数カウントプログラム、境界セル情報設定プログラム、境界セル係数設定プログラム、内部セル係数設定プログラム、行列演算プログラム、収束判定プログラム、及び、ＵＩプログラムに含まれる命令を実行することによって、コンピュータ１を切断面数カウント部１２、境界セル情報生成部１３、内部セル情報生成部１４、境界セル係数設定部１５、内部セル係数設定部１６、行列演算部１７、及び、収束判定部１８として機能させる。
【００９３】
本発明の目的は、上記各プログラムのプログラムコード（実行形式プログラム、中間コードプログラム、ソースプログラム）をコンピュータで読み取り可能に記録した記録媒体を数値流体計算装置１に供給し、数値流体計算装置１がその記録媒体に記録されているプログラムコードを読み出し実行することによっても達成可能である。
【００９４】
上記記録媒体としては、例えば、磁気テープやカセットテープ等のテープ系、フロッピー（登録商標）ディスク／ハードディスク等の磁気ディスクやＣＤ−ＲＯＭ／ＭＯ／ＭＤ／ＤＶＤ／ＣＤ−Ｒ等の光ディスクを含むディスク系、ＩＣカード（メモリカードを含む）／光カード等のカード系、あるいはマスクＲＯＭ／ＥＰＲＯＭ／ＥＥＰＲＯＭ／フラッシュＲＯＭ等の半導体メモリ系などを用いることができる。
【００９５】
また、数値流体計算装置１を通信ネットワークと接続可能に構成し、上記プログラムコードを通信ネットワークを介して数値流体計算装置１に供給するようにしてもよい。この通信ネットワークとしては、とくに限定されず、例えば、インターネット、イントラネット、エキストラネット、ＬＡＮ、ＩＳＤＮ、ＶＡＮ、ＣＡＴＶ通信網、仮想専用網（virtual private network）、電話回線網、移動体通信網、衛星通信網等が利用可能である。また、通信ネットワークを構成する伝送媒体としては、とくに限定されず、例えば、ＩＥＥＥ１３９４、ＵＳＢ、電力線搬送、ケーブルＴＶ回線、電話線、ＡＤＳＬ回線等の有線でも、ＩｒＤＡやリモコンのような赤外線、Ｂｌｕｅｔｏｏｔｈ（登録商標）、８０２．１１無線、ＨＤＲ、携帯電話網、衛星回線、地上波デジタル網等の無線でも利用可能である。なお、本発明は、上記プログラムコードが電子的な伝送で具現化された、搬送波に埋め込まれたコンピュータデータ信号の形態でも実現され得る。
【００９６】
（求積公式）
境界セルにおける境界面ｎｂの面積Ａ_nbは、上述したように、境界面の構成毎に導出された求積公式を用いて算出することができる。以下、この求積公式について、図５〜図８を参照して説明する。なお、以下の説明においては、基礎となる立方体の面に占める境界面ｎｂの割合を開口率Ａ_fluidと記す。
【００９７】
「基礎となる立方体の各辺に高々１つの切断点を設ける」という条件から、考慮すべき組み合わせは、（１）切断点がない場合、（２）切断点が２つある場合、（３）切断点が３つある場合、（４）切断点が４つある場合の何れかに分類される。以下（１）〜（４）の場合について順に説明する。
【００９８】
（１）切断点がない場合
境界面ｎｂ上に切断点がない場合、境界面ｎｂ全体が対象物Ωの内部（流体）に含まれるか（図５（ａ）参照）、又は、境界面ｎｂ全体が対象物Ωの外部（固体）に含まれるか（図５（ｂ）参照）の何れかである。前者の場合、開口率Ａ_fluid＝１であり、後者の場合、開口率Ａ_fluid＝０である。
【００９９】
（２）切断点が２つある場合
境界面ｎｂが２つの切断面を含む場合、境界面ｎｂの構成は、図６（ａ）〜図６（ｃ）に示す３通りのうちの何れかになる。この場合、開口率Ａ_fluidは、式（１６）により算出される。
【０１００】
【数１６】
【０１０１】
（３）切断点が３つある場合
境界面ｎｂ上に切断点が３つある場合、境界面ｎｂの構成は、図７（ａ）〜図７（ｅ）に示す５通りのうちの何れかになる。図７（ａ）〜図７（ｄ）に示す構成の場合、開口率Ａ_fluidは、式（１７）により算出される。
【０１０２】
【数１７】
【０１０３】
一方、図７（ｅ）の構成は、問題が２次元に退化してしまうため望ましくない。このため、図７（ｅ）に示す切断点の配置が生じた場合には、メッシュを切り直す（対象物Ωを再分割する）ことが望ましい。
【０１０４】
（４）切断点が４つある場合
境界面ｎｂ上に切断点が４つある場合、境界面ｎｂの構成は、図８（ａ）〜図８（ｄ）に示す４通りのうちの何れかになる。図８（ａ）〜図８（ｂ）に示す構成の場合、開口率Ａ_fluidは、式（１８）により算出される。
【０１０５】
【数１８】
【０１０６】
なお、図８（ｃ）〜図８（ｄ）に示す構成は、問題が２次元に退化してしまうため望ましくない。このため、図８（ｃ）〜図８（ｄ）に示す切断点の配置が生じた場合には、メッシュを切り直す（対象物Ωを再分割する）ことが望ましい。
【０１０７】
〔実施例〕
最後に、温度Ｔが従う対流拡散方程式に本発明を適用した実施例について説明する。
【０１０８】
本実施例において計算対象とする系を図９に示す。本実施例において計算対象とする系は、図９に示すように、共通の中心を持つ小球面（半径ｄ_i）と大球面（半径ｄ_o）との間の空間である（ｄ_i＜ｄ_o）。境界条件としては、小球面（内側境界）の温度＝Ｔ_i、大球面（外側境界）の温度＝Ｔ_oを課す（Ｔ_i＞Ｔ_o）。
【０１０９】
レーリー数Ｒａを下表の各値に設定したうえで、本発明の数値流体計算方法に従って各セルにおける温度Ｔの値を算出した。なお、レーリー数Ｒａとは、式（１９）により定義される量である。
【０１１０】
【数１９】
【０１１１】
【表１】
【０１１２】
計算精度をチェックするために、各球面上での平均Nusselt数（以下「Ｎｕ数」」と略記する）を評価した。なお、小球面上での平均Nusselt数は、式（２０ａ）のように定義される量であり、大球面上での平均Nusselt数は、式（２０ｂ）のように定義される量である。
【０１１３】
【数２０】
【０１１４】
評価結果を図１０に示す。図１０（ａ）は、小球面上での平均Nusselt数を示し、図１０（ｂ）は、大球面上での平均Nusselt数を示している。同図においては、下記参考文献に記載された平均Nusselt数を比較例として示している。本発明の数値流体計算方法によれば、表１に示した何れのレーリー数に対しても、精度の高い数値流体計算が実現されていることが分かる。
【０１１５】
〔参考文献〕Greg V K, "Natural convection between concentric spheres", International Journal of Heat and Mass Transfer, 35（8), 1991, pp1935-1945
【産業上の利用可能性】
【０１１６】
本発明は数値流体計算に広く利用することができる。特に、（Ｖ）ＣＡＤシステムと連携する数値流体計算に好適に利用することができる。
【符号の説明】
【０１１７】
１ 数値流体計算装置
１１ モデリング部
１２ 切断面数カウント部
１３ 境界セル情報生成部
１４ 内部セル情報生成部
１５ 境界セル係数設定部
１６ 内部セル係数設定部
１７ 行列演算部
１８ 収束判定部
１９ ＵＩ部
２１ パラメータ記憶部
２２ メッシュデータ記憶部
ＩＣ 内部セル
ＢＣ 境界セル
ｅ，ｗ，ｎ，ｓ，ｔ，ｂ 境界面
ｋｔ 三角切断面
Ｐ 注目セルに含まれる格子点
Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂ 隣接セルに含まれる格子点
【特許請求の範囲】
【請求項１】
対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を、コンピュータによって解く数値流体計算方法において、
計算対象領域を複数のセルに分割するセル分割ステップであって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割ステップと、
上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出ステップと、
上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出ステップにて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定ステップと、を含んでいる、
ことを特徴とする数値流体計算方法。
【請求項２】
上記境界セルの各々について、該境界セルに隣接する隣接セルとの間の境界面の面積を算出する境界面積算出ステップと、
上記境界セルの各々について、上記代数方程式に含まれる係数であって、該境界セルに隣接する隣接セル上の上記物理量に係る係数を、上記境界面積算出ステップにて算出された各境界面の面積を用いて算出する係数算出ステップとを更に含んでいる、
ことを特徴とする請求項１に記載の数値流体計算方法。
【請求項３】
上記該境界セルの各々に含まれる三角切断面数Ｎの上限値Ｎmaxが４である、
ことを特徴とする請求項１または２に記載の数値流体計算方法。
【請求項４】
上記対流拡散方程式は、（１）式により表され、
【数１】
上記代数方程式は、各セル上の上記物理量をφP、該セルに隣接する隣接セル上の上記物理量をφE，φW，φN，φS，φT，φBとして、（２）式により表される、
【数２】
ことを特徴とする請求項１から３までの何れか１項に記載の数値流体計算方法。
【請求項５】
対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を解く数値流体計算装置において、
計算対象領域を複数のセルに分割するセル分割手段であって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割手段と、
上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出手段と、
上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出手段にて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定手段と、を備えている、
ことを特徴とする数値流体計算装置。
【請求項６】
コンピュータを請求項５に記載の数値流体計算装置として動作させるためのプログラムであって、上記コンピュータを上記数値流体計算装置が備える各手段として機能させるプログラム。
【請求項７】
請求項６に記載のプログラムが記録された、コンピュータ読み取り可能な記録媒体。
【請求項１】
対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を、コンピュータによって解く数値流体計算方法において、
計算対象領域を複数のセルに分割するセル分割ステップであって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割ステップと、
上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出ステップと、
上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出ステップにて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定ステップと、を含んでいる、
ことを特徴とする数値流体計算方法。
【請求項２】
上記境界セルの各々について、該境界セルに隣接する隣接セルとの間の境界面の面積を算出する境界面積算出ステップと、
上記境界セルの各々について、上記代数方程式に含まれる係数であって、該境界セルに隣接する隣接セル上の上記物理量に係る係数を、上記境界面積算出ステップにて算出された各境界面の面積を用いて算出する係数算出ステップとを更に含んでいる、
ことを特徴とする請求項１に記載の数値流体計算方法。
【請求項３】
上記該境界セルの各々に含まれる三角切断面数Ｎの上限値Ｎmaxが４である、
ことを特徴とする請求項１または２に記載の数値流体計算方法。
【請求項４】
上記対流拡散方程式は、（１）式により表され、
【数１】
上記代数方程式は、各セル上の上記物理量をφP、該セルに隣接する隣接セル上の上記物理量をφE，φW，φN，φS，φT，φBとして、（２）式により表される、
【数２】
ことを特徴とする請求項１から３までの何れか１項に記載の数値流体計算方法。
【請求項５】
対流拡散方程式を各セル上で積分することによって得られる代数方程式であって、各セル上の物理量と、該セルに隣接する隣接セル上の上記物理量との間に成り立つ代数方程式からなる方程式系を解く数値流体計算装置において、
計算対象領域を複数のセルに分割するセル分割手段であって、計算対象領域の境界を含む境界セルとして、各辺上に高々１つの切断点が設けられた６面体を３つの切断点を頂点とするＮ個の三角切断面で切断して得られる（６＋Ｎ）面体（ＮはＮmax以下の自然数）を用いるセル分割手段と、
上記境界セルの各々について、上記Ｎ個の三角切断面の面積を算出する切断面積算出手段と、
上記境界セルの各々について、上記代数方程式に含まれる、対流流束の各三角切断面上での積分値、及び、拡散流束の各三角切断面上での積分値を、上記切断面積算出手段にて算出された各三角切断面の面積と、上記計算対象領域の境界に課された境界条件とから算出する境界条件設定手段と、を備えている、
ことを特徴とする数値流体計算装置。
【請求項６】
コンピュータを請求項５に記載の数値流体計算装置として動作させるためのプログラムであって、上記コンピュータを上記数値流体計算装置が備える各手段として機能させるプログラム。
【請求項７】
請求項６に記載のプログラムが記録された、コンピュータ読み取り可能な記録媒体。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１１−１３８２１０（Ｐ２０１１−１３８２１０Ａ）
【公開日】平成２３年７月１４日（２０１１．７．１４）
【国際特許分類】
【出願番号】特願２００９−２９６２５４（Ｐ２００９−２９６２５４）
【出願日】平成２１年１２月２５日（２００９．１２．２５）
【出願人】（５０３３５９８２１）独立行政法人理化学研究所 (1,056)
【Ｆターム（参考）】
【公開日】平成２３年７月１４日（２０１１．７．１４）
【国際特許分類】
【出願日】平成２１年１２月２５日（２００９．１２．２５）
【出願人】（５０３３５９８２１）独立行政法人理化学研究所 (1,056)
【Ｆターム（参考）】
[ Back to top ]