ホーム > 特許ランキング > 国立大学法人九州大学
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-04-14
(45)【発行日】2023-04-24
(54)【発明の名称】数値計算方法および数値計算プログラム
(51)【国際特許分類】
G06F 30/23 20200101AFI20230417BHJP
G06F 30/28 20200101ALI20230417BHJP
【ＦＩ】
G06F30/23
G06F30/28
【請求項の数】
8
(21)【出願番号】P 2021501394
(86)(22)【出願日】2019-02-25
(86)【国際出願番号】 JP2019007071
(87)【国際公開番号】W WO2020174536
(87)【国際公開日】2020-09-03
【審査請求日】2022-01-14
【新規性喪失の例外の表示】特許法第３０条第２項適用 平成３１年１月２８日、東京大学安田講堂（東京都文京区本郷７丁目３－１）にて開催された、平成３０年度ＳＩＰ「革新的燃焼技術」公開シンポジウム、及び、平成３１年２月１９日に公開された、当該シンポジウムの資料のウェブサイト（ｈｔｔｐｓ：／／ｗｗｗ．ｊｓｔ．ｇｏ．ｊｐ／ｓｉｐ／ｄｌ／ｋ０１／ｓｙｍｐｏｅｎｄ／０４ｓ＿１３．ｐｄｆ）にて発表。
【新規性喪失の例外の表示】特許法第３０条第２項適用 平成３０年７月１２日、ＴＨＥ ＧＲＡＮＤ ＨＡＬＬ（品川）（東京都港区港南２丁目１６－４）にて開催された、ＪＨＰＣＮ：学際大規模情報基盤共同利用・共同研究拠点（第１０回シンポジウム）、及び、当該シンポジウムの資料のウェブサイト（ｈｔｔｐｓ：／／ｊｈｐｃｎ－ｋｙｏｔｅｎ．ｉｔｃ．ｕ－ｔｏｋｙｏ．ａｃ．ｊｐ／ｊａ／ｓｙｍｐｏ／１０ｔｈ／Ｉｎｔｒｏ／ＥＸ１７３２０＿Ｐｏｓｔｅｒ．ｐｄｆ）にて発表。
【新規性喪失の例外の表示】特許法第３０条第２項適用 平成３０年８月２７日発行の、東京大学情報基盤センター・スーパーコンピューティングニュース Ｖｏｌ．２０，Ｎｏ．Ｓｐｅｃｉａｌ Ｉｓｓｕｅ ２、及び、そのオンライン版（ｈｔｔｐｓ：／／ｗｗｗ．ｃｃ．ｕ－ｔｏｋｙｏ．ａｃ．ｊｐ／ｐｕｂｌｉｃ／ＶＯＬ２０／ｓｐｅｃｉａｌ２／０２．２０１８０８ＳＰ－１＿３．ｐｄｆ）にて発表。
【新規性喪失の例外の表示】特許法第３０条第２項適用 平成３０年５月１１日、九州大学情報基盤研究開発センター（福岡県福岡市西区元岡７４４）にて開催された、先駆的科学計算に関するフォーラム２０１８、及び、当該フォーラムの資料のウェブサイト（ｈｔｔｐｓ：／／ｗｗｗ．ｃｃ．ｋｙｕｓｈｕ－ｕ．ａｃ．ｊｐ／ｓｃｐ／ｄｏｃ／ｕｓｅｒｓ／ｆｏｒｕｍ／ｆｏｒｕｍ２０１８０５１１／ｆｕｋａｇａｗａ．ｐｄｆ）にて発表。
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成３０年度、国立研究開発法人科学技術振興機構、戦略的イノベーション創造プログラム事業「自動車用内燃機関摺動面潤滑モデルの確立および設計支援ソフトウェアへの展開」委託研究、産業技術力強化法第１９条の適用を受ける特許出願
(73)【特許権者】
【識別番号】504145342
【氏名又は名称】国立大学法人九州大学
(74)【代理人】
【識別番号】100193378
【弁理士】
【氏名又は名称】野口 明生
(72)【発明者】
【氏名】深川 宏樹
【審査官】松浦 功
(56)【参考文献】
【文献】特開２００５－３３９４９２（ＪＰ，Ａ）
【文献】向井信彦 外４名，粒子法による大動脈弁の開閉シミュレーション，日本バーチャルリアリティ学会論文誌，特定非営利活動法人日本バーチャルリアリティ学会，2013年，第18巻，第4号，pp. 447-453
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｆ ３０／００ －３０／２８
Ｇ１６Ｚ ９９／００
Ｇｏｏｇｌｅ Ｓｃｈｏｌａｒ
(57)【特許請求の範囲】
【請求項1】
複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータを用いて、前記階層型の記憶装置に記憶されている第１時刻における流体および構造体の物理量から微小時間後の第２時刻における前記流体および前記構造体の物理量を連成解析する数値計算方法であって、前記第１時刻のみにおける流体の物理量から前記第２時刻における前記流体の物理量を流体の支配方程式を用いて計算する処理と、
前記第１時刻のみにおける構造体の物理量から前記第２時刻における前記構造体の物理量を構造体の支配方程式を用いて計算する処理と、
前記第１時刻のみにおける流体の圧力と構造体の応力を対応させて、前記第２時刻における流体と構造体との境界面の位置を境界面の支配方程式を用いて計算する処理と、
を含むことを特徴とする数値計算方法。
【請求項2】
前記流体の支配方程式は、ナビエ・ストークス方程式および連続の方程式であることを特徴とする請求項１に記載の数値計算方法。【請求項3】
前記構造体の支配方程式は、応力と歪の関係を表す構成方程式であることを特徴とする請求項１または請求項２に記載の数値計算方法。【請求項4】
前記境界面の支配方程式は、速度場の時間微分と前記流体の圧力および前記構造体の応力の釣り合いとの関係を表す方程式であることを特徴とする請求項１から請求項３のいずれか１項に記載の数値計算方法。【請求項5】
複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータに対して、前記階層型の記憶装置に記憶されている第１時刻における流体および構造体の物理量から微小時間後の第２時刻における前記流体および前記構造体の物理量を連成解析する数値計算を指令する数値計算プログラムであって、前記第１時刻のみにおける流体の物理量から前記第２時刻における前記流体の物理量を流体の支配方程式を用いて計算する処理と、
前記第１時刻のみにおける構造体の物理量から前記第２時刻における前記構造体の物理量を構造体の支配方程式を用いて計算する処理と、
前記第１時刻のみにおける流体の圧力と構造体の応力を対応させて、前記第２時刻における流体と構造体との境界面の位置を境界面の支配方程式を用いて計算する処理と、
を含むことを特徴とする数値計算プログラム。
【請求項6】
前記支配方程式の全てが、場の時間微分を、前記場および前記場の一階空間微分ならびに二階空間微分で表現したものであり、前記支配方程式の全てを時刻および位置で離散化した方程式を用いて処理をすることを特徴とする請求項５に記載の数値計算プログラム。【請求項7】
前記離散化した方程式を前記複数の演算器に割り当てて並列処理することを特徴とする請求項６に記載の数値計算プログラム。【請求項8】
前記離散化した方程式は、適切な関数Ｆを用いて下記式の形で表現されることを特徴とする請求項７に記載の数値計算プログラム。ただし、ｍは時刻の離散化のインデックスであり、ｌは位置の離散化のインデックスであり、ｊ（ｌ）はインデックスｌの位置の周辺の点に割り振られたインデックスとしたときに、φ（ｌ，ｍ）は、インデックスｌの位置におけるインデックスｍの時刻の場を表す。【数1】
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、数値計算方法および数値計算プログラムに関する。
【背景技術】
【0002】
現在の製品設計ではＣＡＤ（Ｃｏｍｐｕｔｅｒ－Ａｉｄｅｄ Ｄｅｓｉｇｎ）が一般的に用いられている。さらに、製品設計で作成されるＣＡＤデータを研究・開発工程におけるエンジニアリングにも用いる、いわゆるＣＡＥ（Ｃｏｍｐｕｔｅｒ－Ａｉｄｅｄ ｅｎｇｉｎｅｅｒｉｎｇ）も実用化されてきた。ＣＡＤデータを用いてシミュレーションを行うことで製品性能の事前検証ができれば研究・開発工程を高速化できるからである。
【0003】
工業製品の検証では流体構造連成解析が必要とされる。流体構造連成解析とは、流体の流動と構造体の変形の相互作用を考慮した解析であり、例えば軸受の事前検証では、潤滑油の流れで生じる力で軸受が変形する様子を解析する必要がある。この相互作用に関し、従来の流体構造連成解析では、流体と構造体で別個に方程式を解き、境界面を修正して解を収束させる反復法が用いられてきた（例えば特許文献１参照）。
【先行技術文献】
【特許文献】
【0004】
【文献】特開２００６－７２５６６号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
しかしながら、反復法を用いた流体構造連成解析で解を収束させるためには、流体と構造体の方程式を反復して計算する必要がある。一方、現在のコンピュータの主流は、複数の演算器（コア）と階層構造のメモリを採用しており、メモリバンド幅が演算性能に比べて低い（Ｂｙｔｅ／Ｆｌｏｐが低い）。このような構成のコンピュータの演算性能を引き出すにはレジスタやキャッシュメモリを有効に使い、メモリにおけるデータ通信量を抑えることができる参照の局所性が高い計算方法を用いる必要がある。この観点からは、反復法を用いた流体構造連成解析は、別個に解いた方程式の解における誤差を比較する必要があり、現在の主流となっている計算機のアーキテクチャに適するものではない。
【0006】
本発明は、上記に鑑みてなされたものであり、その目的は、参照の局所性が高い流体構造連成解析の数値計算方法および数値計算プログラムを提供することである。
【課題を解決するための手段】
【0007】
上記課題を解決するために、本発明の一態様にかかる数値計算方法は、複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータを用いて、前記階層型の記憶装置に記憶されている第１時刻における流体および構造体の物理量から微小時間後の第２時刻における前記流体および前記構造体の物理量を連成解析する数値計算方法であって、前記第１時刻のみにおける流体の物理量から前記第２時刻における前記流体の物理量を流体の支配方程式を用いて計算する処理と、前記第１時刻のみにおける構造体の物理量から前記第２時刻における前記構造体の物理量を構造体の支配方程式を用いて計算する処理と、前記第１時刻のみにおける流体の圧力と構造体の応力を対応させて、前記第２時刻における流体と構造体との境界面の位置を境界面の支配方程式を用いて計算する処理と、を含むことを特徴とする。
【0008】
また、本発明の一態様にかかる数値計算プログラムは、複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータに対して、前記階層型の記憶装置に記憶されている第１時刻における流体および構造体の物理量から微小時間後の第２時刻における前記流体および前記構造体の物理量を連成解析する数値計算を指令する数値計算プログラムであって、前記第１時刻のみにおける流体の物理量から前記第２時刻における前記流体の物理量を流体の支配方程式を用いて計算する処理と、前記第１時刻のみにおける構造体の物理量から前記第２時刻における前記構造体の物理量を構造体の支配方程式を用いて計算する処理と、前記第１時刻のみにおける流体の圧力と構造体の応力を対応させて、前記第２時刻における流体と構造体との境界面の位置を境界面の支配方程式を用いて計算する処理と、を含むことを特徴とする。
【発明の効果】
【0009】
本発明によれば、参照の局所性が高い流体構造連成解析の数値計算方法および数値計算プログラムを提供することができる。
【図面の簡単な説明】
【0010】
【発明を実施するための形態】
【0011】
以下、図面を参照しながら、本発明の実施形態にかかる数値計算方法および数値計算プログラムを詳細に説明する。ただし、以下の説明で参照される図面は模式的なものであり、寸法またはその比率が実際のものとは異なる場合がある。
【0012】
〔装置構成〕
図１は、共有メモリ型の階層構造の概念を示す図であり、図２は、分散メモリ型の階層構造の概念を示す図である。本発明の実施形態にかかる数値計算方法および数値計算プログラムは、複数の演算器とこれらが演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータにおいて好適に実施される。そこで、図１および図２を参照しながら、複数の演算器と階層構造の記憶装置を備えるコンピュータにおける計算処理について説明する。なお、複数の演算器を備えるコンピュータは並列計算機や並列コンピュータとも呼ばれ、この並列計算機上の計算処理は並列処理や並列コンピューティングとも呼ばれる。
図１は、共有メモリ型の階層構造の概念を示す図であり、図２は、分散メモリ型の階層構造の概念を示す図である。本発明の実施形態にかかる数値計算方法および数値計算プログラムは、複数の演算器とこれらが演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータにおいて好適に実施される。そこで、図１および図２を参照しながら、複数の演算器と階層構造の記憶装置を備えるコンピュータにおける計算処理について説明する。なお、複数の演算器を備えるコンピュータは並列計算機や並列コンピュータとも呼ばれ、この並列計算機上の計算処理は並列処理や並列コンピューティングとも呼ばれる。
【0013】
図１に示される概念図には、複数の演算器１０１～１０ｎとこれらが演算に用いるデータを記憶するための階層型の記憶装置１１とを備える並列計算機１００が記載されている。階層型の記憶装置１１とは、各演算器１０１～１０ｎと主記憶装置１２との間に階層化された記憶装置を配したものである。図１に示される例では、演算器１０１～１０ｎから順に、レジスタ１３１～１３ｎ、Ｌ１キャッシュメモリ１４１～１４ｎ、Ｌ２キャッシュメモリ１５１～１５ｍ、主記憶装置１２の構成を有している。なお、図１に示される構成は、一つの例であり、階層型の記憶装置１１の構成はこれに限定されるものではなく、階層の程度としても、Ｌ２キャッシュメモリ１５１～１５ｍと主記憶装置１２との間にＬ３キャッシュメモリを備えるなどの変形例がある。
【0014】
これらの記憶装置の処理速度は、主記憶装置１２、Ｌ２キャッシュメモリ１５１～１５ｍ、Ｌ１キャッシュメモリ１４１～１４ｎ、レジスタ１３１～１３ｎの順で高速であり、メモリのデータ転送能力であるメモリバンド幅は、演算器１０１～１０ｎの演算性能に比べて低い（Ｂｙｔｅ／Ｆｌｏｐが低い）。このような特性から、参照の局所性を高めることが演算性能を発揮させるために重要となる。
【0015】
なお、参照の局所性が高い状態とは、一度参照されたデータが再度参照されることや当該データに近いアドレスのデータが参照される可能性が高いという状態である。言い換えると、演算器１０１～１０ｎがレジスタ１３１～１３ｎまたはＬ１キャッシュメモリ１４１～１４ｎもしくはＬ２キャッシュメモリ１５１～１５ｍなど上位の記憶装置に記憶されているデータを用いて計算を行うことができる状態である。
【0016】
逆に言えば、演算器１０１～１０ｎが主記憶装置１２に記憶されているデータを用いなければいけない状況を少なくする戦略が必要である。例えば、そのような状況は、演算器１０１が計算した結果を用いて演算器１０ｎが計算をしなければいけないときに起きる。また、演算器１０１が計算した結果を用いて演算器１０ｎが計算をしなければいけない場合には、演算器１０１の計算が終わらなければ演算器１０ｎの計算を開始することができず、待ち時間が発生してしまうという問題も発生してしまう。
【0017】
図２は、分散メモリ型の階層構造の概念を示す図であり、大型計算機システムなどで採用されている。図２に示される並列計算機２００では、各ノード１６１～１６ｎが複数の演算器とこれらが演算に用いるデータを記憶するための階層型の記憶装置とを備えており、さらに、各ノード１６１～１６ｎにおける主記憶装置１２１～１２ｎがネットワーク１７を介して互いに共有されている。ここで、主記憶装置１２１～１２ｎの共有とは、主記憶装置１２１～１２ｎに単一のアドレスを付与することで、各ノード１６１～１６ｎの演算器が、主記憶装置１２１～１２ｎに記憶されているデータを相互に参照することができることをいう。
【0018】
上記のような構成のコンピュータでも、参照の局所性を高めることが演算性能を発揮させるために重要となる。特に、あるノード１６１～１６ｍの演算器が、他のノード１６１～１６ｍの主記憶装置１２１～１２ｎに記憶されているデータを参照するような状況を可能な限り少なくすることが重要である。
【0019】
上記観点から、反復法を用いた流体構造連成解析は、現在の主流となっている計算機のアーキテクチャに適するものではない。反復法を用いた流体構造連成解析は、別個に解いた方程式の解における誤差を比較する必要があり、参照の局所性が低いからである。そこで、本発明の実施形態にかかる数値計算方法および数値計算プログラムは、陽解法を用いた流体構造連成解析を行う。以下、陽解法に関する説明を行う。
【0020】
〔陽解法〕
陽解法とは、第１時刻tにおける状態量のみから微小時間後の第2時刻t+Δtにおける位置xの状態量を計算する方法である。一方、陰解法は、第２時刻t+Δtにおける位置xの状態量を計算するために、第２時刻t+Δtにおける他の位置xiの状態量も用いる方法である。端的には、陰解法は、第２時刻t+Δtにおける位置xの状態量を計算するために、第２時刻t+Δtにおける他の位置xiの状態量も含んだ連立方程式を解く必要があり、陽解法では、このような連立方程式を解く必要がない。
陽解法とは、第１時刻tにおける状態量のみから微小時間後の第2時刻t+Δtにおける位置xの状態量を計算する方法である。一方、陰解法は、第２時刻t+Δtにおける位置xの状態量を計算するために、第２時刻t+Δtにおける他の位置xiの状態量も用いる方法である。端的には、陰解法は、第２時刻t+Δtにおける位置xの状態量を計算するために、第２時刻t+Δtにおける他の位置xiの状態量も含んだ連立方程式を解く必要があり、陽解法では、このような連立方程式を解く必要がない。
【0021】
このように、陽解法では、連立方程式を解く必要がないので、時間ステップ当たりの計算量およびメモリ転送量が少ない。さらに、第２時刻t+Δtにおける位置xの状態量を計算するために、第１時刻tにおける位置xの近傍のみの状態量から計算するようにすれば、より一層、参照の局所性が高まる。そこで、以下で参照の局所性が高い陽解法を適用し得る条件について説明する。
【0022】
結論としては、ある領域の場φ(n,t)の時間発展が、以下の式（１）に示すように、場φ(n,t)と当該場の空間に関する一階微分φ′(n,t)と二階微分φ″(n,t)の関数fで表現することができればよい。なお、位置をn=(n1,n2,n3)とし、時刻をtとしている。
【0023】
【数1】
【0024】
上記式（１）のように表現される場φ(n,t)を時刻および位置に関して離散化する。時刻tの離散化は、自然数mを用いてt(m+1)=t(m)+Δtとする。位置ｎの離散化は、自然数lから位置n=(n1,n2,n3)への関数n(l)とする。そして、自然数lに対応した位置n(l)の周辺の点を自然数の組(l,j)から自然数への関数k(l,j)を用いてn(k(l,j))と表す。
【0025】
すると、離散化された場の一階微分φ′(n(l),t(m))と二階微分φ″(n(l),t(m))は、位置n(l)の周辺の点における場φ(n(k(l,j)),t(m))の関数で近似できる。したがって、f(φ(n,t),φ′(n,t),φ″(n,t))を離散化したものは、F(φ(n(k(l,j)),t(m)))と近似できる。
【0026】
結果、上記式（１）は、以下のように離散化される。
【0027】
【数2】
【0028】
また、上記式（２）における位置n(l)と時間t(m)を離散化に用いたインデックス(l,m)で同一すれば、以下のように表現される。なお、j(l)は、これはインデックスlの位置n(l)の周辺の点に割り振られたインデックスである。
【0029】
【数3】
【0030】
上記式（３）から解るように、φ(l,m+1)がφ(l,m)とφ(j(l),m)で定まる。ここで重要なことは、式（３）の右辺には、m+1が含まれないことである。このことは、第１時刻t(m)における状態量のみから微小時間後の第2時刻t(m+1)における位置n(l)の状態量を計算することができることを意味している。
【0031】
また、式（３）の右辺には、j(l)が含まれているが、これはインデックスlの位置n(l)の周辺の点に割り振られたインデックスとしたのだから、φ(j(l),m)は、第１時刻t(m)における位置n(l)の周辺の点の場φの値である。図３は、参照の局所性が高い陽解法の概念を示す図である。なお、図３では、位置および時間をインデックスで同一視して表記している。上記式（３）に示される方程式は、図３に示されるように、時刻m+1における位置lのφ(l,m+1)を求める際に、時刻mにおける位置lのφ(l,m)とその近傍の位置j-1(l),j+1(l)等におけるφ(j(l),m)を用いるだけでよいことを意味している。インデックスlは、φの値をメモリに記憶する際のアドレスの近さに対応するので、φ(l,m+1)の計算に際し、φ(l,m)とφ(j(l),m)という参照の局所性が高いデータを用いていることになる。
【0032】
以上の議論から解るように、上記式（１）のように表現できる場合、参照の局所性が高い陽解法を用いることができる。したがって、陽解法を用いた流体構造連成解析を行うためには、流体の支配方程式、構造体の支配方程式、および流体と構造体の境界面の支配方程式を上記式（１）のように表現することが次の目標となる。
【0033】
〔流体の支配方程式〕
流体の支配方程式には、連続の式とナビエ・ストークス方程式を用いる。連続の式は、下記式（４）のように表現され、ナビエ・ストークス方程式は下記式（５）のように表現される。
流体の支配方程式には、連続の式とナビエ・ストークス方程式を用いる。連続の式は、下記式（４）のように表現され、ナビエ・ストークス方程式は下記式（５）のように表現される。
【0034】
【数4】
【数5】
【0035】
ここで、ρは流体内の質量密度であり、*ρは、その双対ホッジであり、要素内の総質量を表す。速度場u#に関するＬｉｅ微分をLu#と書く。計量テンソルをg( , )とし、u:=g(u#, )と定める。圧力Pは質量密度ρの関数とする。なお、Δはホッジラプラシアンとし、すなわちΔ=d*d*-*d*dである。
【0036】
上記式（４）（５）と式（１）を比較すると解るように、式（４）（５）は、局所性が高い陽解法を適用し得る式（１）の形を有している。したがって、式（４）（５）に陽解法を適用し、*ρとｕの時間発展を計算することができる。本発明の実施形態にかかる数値計算方法および数値計算プログラムは、第１時刻ｔにおける流体の物理量*ρおよびｕから微小時間後の第2時刻t+Δtにおける流体の物理量*ρおよびｕを上記（４）（５）を用いて計算する。
【0037】
〔構造体の支配方程式〕
構造体の時間発展を記述するためには、構造体の変形に伴う体積素片の位置を追跡する。このとき、初期時刻における体積素片の位置を示す座標値を用いて各体積素片をラベル付けすることができる。このように構造体の体積素片に固定された座標(t,n1,n2,n3)をラグランジュ座標もしくは自然座標という。一方、空間に固定された通常の空間座標(t,x1,x2,x3)はオイラー座標といわれる。ラグランジュ座標の各座標を時刻ごとにオイラー座標に変換する関数Xi(t,n)は、構造体の時間発展を一意に定める。
構造体の時間発展を記述するためには、構造体の変形に伴う体積素片の位置を追跡する。このとき、初期時刻における体積素片の位置を示す座標値を用いて各体積素片をラベル付けすることができる。このように構造体の体積素片に固定された座標(t,n1,n2,n3)をラグランジュ座標もしくは自然座標という。一方、空間に固定された通常の空間座標(t,x1,x2,x3)はオイラー座標といわれる。ラグランジュ座標の各座標を時刻ごとにオイラー座標に変換する関数Xi(t,n)は、構造体の時間発展を一意に定める。
【0038】
ラグランジュ座標における計量をg=gijdni×dnjとすると、構造体の歪みの時間変化はラグランジュ座標における計量テンソルgijの時間変化を用いて下記式（６）で表すことができる。なお、構造体の応力テンソルki
jは、ei
j:=gjkekjの関数となる。
【0039】
【数6】
【0040】
オイラー座標におけるＭｕｓｉｃａｌ ｉｓｏｍｏｒｐｈｉｓｍをオイラー座標における計量δijdxi×dxjを用いて行う。すなわち、Xi:=δijXjである。
【0041】
構造体の支配方程式では、ナビエ・ストークス方程式に対応する速度場として下記式（７）を用いる。
【0042】
【数7】
【0043】
この速度場に関する構造体の支配方程式は、下記式（８）のように表現される。ただし、下式におけるsは変位ベクトルsi=δij(Xj-nj)である。
【0044】
【数8】
【0045】
ここで、上記式（８）における圧力Ｐは、構造体の歪みeを用いて以下のように計算できる。
【0046】
【数9】
【0047】
ここで、ラメ定数λとμは、ヤング率Ｅとポアソン比νを用いて、それぞれ以下のように表すことができる。
【0048】
【数10】
【0049】
なお、上記式（７）を座標表示した場合は下記式（１１）のように表現される。
【0050】
【数11】
【数12】
【0051】
上記式（８）（１１）と式（１）を比較すると解るように、式（８）（１１）は、局所性が高い陽解法を適用し得る式（１）の形を有している。したがって、式（８）（１１）に陽解法を適用し、速度場ｕの時間発展を計算することができる。本発明の実施形態にかかる数値計算方法および数値計算プログラムは、第１時刻ｔにおける構造体の物理量ｕから微小時間後の第2時刻t+Δtにおける構造体の物理量ｕを上記式（８）（１１）を用いて計算する。
【0052】
〔境界面の支配方程式〕
流体と構造体の境界面では、境界面の位置の速度場uが下記式（１３）に従うとする。
流体と構造体の境界面では、境界面の位置の速度場uが下記式（１３）に従うとする。
【0053】
【数13】
【0054】
ただし、流体側のＰは、流体の密度の関数としての流体の圧力を用い、構造体側のＰは、先述の式（９）で表される歪みeの関数としての圧力を用いる。
【0055】
上記式（１３）と式（１）を比較すると解るように、式（１３）は、局所性が高い陽解法を適用し得る式の形（１）を有している。したがって、式（１３）に陽解法を適用し、境界面の位置の速度場ｕの時間発展を計算することができる。本発明の実施形態にかかる数値計算方法および数値計算プログラムは、第１時刻ｔにおける流体の圧力と構造体の応力を対応させて、微小時間後の第2時刻t+Δtにおける流体と構造体との境界面の位置を境界面の位置を、上記式（１３）を用いて計算する。
【0056】
以上説明したように、本発明の実施形態にかかる数値計算方法および数値計算プログラムが用いる、流体の支配方程式、構造体の支配方程式、および境界面の支配方程式は、局所性が高い陽解法を適用し得る式（１）の形を有している。したがって、複数の演算器と演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータを用いて、階層型の記憶装置に記憶されている第１時刻における流体および構造体の物理量から微小時間後の第２時刻における流体および構造体の物理量を連成解析する際に、図３に示したような局所性が高い陽解法を適用し得る。
【0057】
また、本発明の実施形態にかかる数値計算方法および数値計算プログラムでは、流体の支配方程式、構造体の支配方程式、および境界面の支配方程式の全てが式（１）の形を有しているので、これらを時刻および位置で離散化した方程式は、式（３）の形をしている。したがって、並列計算機における演算器にとっては、流体の支配方程式と構造体の支配方程式と境界面の支配方程式とを区別することなく均等に処理することができる。例えば、ＭＰＩを使って並列計算する場合には、解析モデルを等分に領域分割することで、ロードバランスをよくできる。
【0058】
〔流体軸受の流体構造連成解析における実施例〕
以下、上記説明した数値計算方法を用いた流体軸受の流体構造連成解析の実施例について説明する。
以下、上記説明した数値計算方法を用いた流体軸受の流体構造連成解析の実施例について説明する。
【0059】
図４は、実施例に用いる流体軸受の模式図である。図４に示されるように、流体軸受２０は、中心に軸２１を配し、軸２１の外側に軸受２２を配し、軸２１と軸受２２との隙間に潤滑油２３を充填した構成を有する。当該構成を有ることにより流体軸受２０は、潤滑油２３の流動によって軸２１と軸受２２との関係が自在に回転可能となっている。
【0060】
本実施例に用いる流体軸受２０の具体的な寸法は、軸２１の半径が３０ｍｍであり、軸受２２の厚さが３０ｍｍであり、軸方向の長さが１０ｍｍであり、軸２１と軸受２２との隙間が３０μｍである。また、軸２１と軸受２２との軸心の偏心量は、１．５μｍであり、軸２１と軸受２２との隙間と偏心量との比である偏心率は０．９５である。なお、図４では、各構成の配置が見やすいように軸２１と軸受２２との隙間における潤滑油２３の膜厚を実際の比率よりも大きく表示してある。
【0061】
軸２１および軸受２２における構造格子の個数は、円周方向に１００個であり、軸方向に５個であり、径方向に５個である。潤滑油２３における構造行使の個数は、径方向に１層であることの他、軸２１および軸受２２の場合と同じである。なお、計算精度を向上させるために、圧力変化が激しい隙間の最狭部では構造格子の大きさを小さくしている。
【0062】
潤滑油２３の特性は、大気圧下での密度ρが１０００ｋｇ／ｍ３であり、粘度ηが１０ｍＰａ・ｓである。軸２１の回転数は１０００ｒｐｍであり、軸２１と潤滑油２３との境界面の移動速度ｕは３．１４ｍ／ｓである。
【0063】
【0064】
図５のように自然座標系(n1,n2,n3)から空間座標系(x1,x2,x3)への座標変換を考えると、各格子の自然座標ｎから空間座標ｘへの写像は、形状関数Nα(n)を用いてxi(n)=xi
αNα(n)で与えられる。ここでxi
αは、節点番号αの節点の空間座標におけるｉ成分の値である。
【0065】
６面体の２次要素では節点が２０個ある。表１は要素の中心点を原点にした場合の節点番号と自然座標との対応を表している。
【0066】
【表1】
【0067】
【表2】
【0068】
ただし、上記表２に記載の関数l(r),m(r)は以下で与えられる。
【0069】
【数14】
【数15】
【0070】
構造格子の各要素の形状の情報は、自然座標上の計量テンソルｇに集約される。この計量テンソルの成分をｇｉｊとすると、体積要素は以下のように計算できる。
【0071】
【数16】
【0072】
一方、計量テンソルｇｉｊは、空間座標上の計量テンソルδｋｌから以下のように計算できる。なお、空間座標がカルテシアン座標の場合にはδｋｌがクロネッカーデルタになる。
【0073】
【数17】
【0074】
上記計量テンソルｇｉｊを用いて、流体の支配方程式（式（４）（５））を表現すると以下のようになる。
【0075】
【数18】
【数19】
【0076】
〔第１の検証〕
以上のような設定の下、本発明の実施形態にかかる数値計算方法の検証を行う。第１の検証では、軸２１および軸受２２が剛体であると仮定し、潤滑油２３の圧力分布と密度分布の計算を行う。
以上のような設定の下、本発明の実施形態にかかる数値計算方法の検証を行う。第１の検証では、軸２１および軸受２２が剛体であると仮定し、潤滑油２３の圧力分布と密度分布の計算を行う。
【0077】
時間刻みを０．１μｓとし、自然座標系での空間刻みを２．０として、流体の支配方程式（式（４）（５）すなわち式（１８）（１９））を離散化して陽解法を適用して潤滑油２３の圧力分布と密度分布の計算を行った。図６は、圧度分布の計算結果を示すグラフであり、図７は、密度分布の計算結果を示すグラフである。
【0078】
図６および図７のいずれも、潤滑油２３の軸方向中央部の値を記載している。また、横軸のθは、周方向の角度であり、隙間が最も狭くなる位置を０とし、最も広くなる位置を±πとしている。なお、潤滑油２３は、θの正方向に流れているとして、周方向の境界条件を周期的境界条件としている。
【0079】
図６および図７に示されるグラフから解るように、大気圧下では１０００ｋｇ／ｍ３であった質量密度が低圧領域では３００ｋｇ／ｍ３まで下がっている。これは、低圧領域では気泡が発生し、平均質量密度が大きく下がっていることを示している。
【0080】
この計算結果の精度を検証するために、本発明の実施形態にかかる数値計算方法の計算結果とレイノルズ方程式をガウス＝ザイデル法の計算結果とを比較する。レイノルズ方程式とは、質量保存の式とナビエ・ストークス方程式から導出される式であり、次式で与えられる。
【0081】
【数20】
【0082】
ここでｈは潤滑油２３の膜厚であり、ｕｍは潤滑油２３を挟む境界面の移動速度の平均値である。潤滑油２３内の圧力Ｐの分布はρ，ｈ，ｕｍの分布から計算できる。ガウス＝ザイデル法では、式（２０）の質量密度ρを定数として圧力Ｐを計算する。そして、計算された圧力Ｐから質量密度ρを計算する。その後、計算された質量密度ρを式（２０）に再代入して圧力Ｐを再計算する。この操作を解が収束するまで繰り返す。
【0083】
キャビテーションを含めるには圧力Ｐが負圧になるとゼロに置き換える操作をする。この方法ではキャビテーションが発生する箇所で連続の式が満たされないので、一般には正しい圧力分布が得られない。しかし、θ＝±πでＰ＝０となる境界条件下では、正しい圧力分布を得ることができることが知られている。
【0084】
このように圧力分布を計算した結果と本発明の実施形態にかかる数値計算方法で圧力分布を計算した結果とを比較する。図８は、圧力分布の計算結果の比較を示すグラフである。図８に示されるグラフから解るように、レイノルズ方程式にガウス＝ザイデル法を用いて計算した圧力分布（－印）とナビエ・ストークス方程式に本発明の実施形態にかかる数値計算方法を用いて計算した圧力分布（＋印）は一致している。
【0085】
〔第２の検証〕
次に、軸受２２が弾性体であるとし、潤滑油２３の圧力分布の計算を行う。上記第１の検証の結果でも解るように、流体軸受では、軸の回転時に潤滑油２３の膜厚の狭小部では圧力が高くなる。結果、弾性体である軸受２２も変形を受ける。しかも、軸受２２の変形は潤滑油２３の膜厚にも影響を与えるので、膜厚の狭小部における潤滑油２３の圧力にも影響を与える。ここでは、軸受２２が剛体であると仮定した場合と弾性体であるとした場合の計算結果を比較する。
次に、軸受２２が弾性体であるとし、潤滑油２３の圧力分布の計算を行う。上記第１の検証の結果でも解るように、流体軸受では、軸の回転時に潤滑油２３の膜厚の狭小部では圧力が高くなる。結果、弾性体である軸受２２も変形を受ける。しかも、軸受２２の変形は潤滑油２３の膜厚にも影響を与えるので、膜厚の狭小部における潤滑油２３の圧力にも影響を与える。ここでは、軸受２２が剛体であると仮定した場合と弾性体であるとした場合の計算結果を比較する。
【0086】
図９は、計算結果の比較を示すグラフである。図９に示されるグラフには、軸受２２が剛体であると仮定した場合の圧力分布が＋印で記載され、軸受２２が弾性体であるとした場合の圧力分布が－印で記載されている。また、図９に示されるグラフには、軸受２２の変形量が●印で併記されている。
【0087】
図９から解るように、軸受２２が弾性体であるとした場合、潤滑油２３の高圧領域では、軸受２２が押し広げられるので、軸受２２が剛体であると仮定した計算結果よりも潤滑油２３の圧力が低くなっている。このことは、本発明の実施形態にかかる数値計算方法が流体と構造体の相互作用を正しく反映させて連成解析を行っていることを示している。
【0088】
〔並列計算機における性能検証〕
次に、本発明の実施形態にかかる数値計算方法を並列計算機に対する指令として実施した場合、すなわち本発明の実施形態にかかる数値計算プログラムの性能を検証する。
次に、本発明の実施形態にかかる数値計算方法を並列計算機に対する指令として実施した場合、すなわち本発明の実施形態にかかる数値計算プログラムの性能を検証する。
【0089】
既述したように、現在主流のコンピュータは、複数の演算器と階層構造のメモリを採用しており、メモリバンド幅が演算性能に比べて低い。このような構成のコンピュータの演算性能を引き出すにはレジスタやキャッシュメモリを有効に使い、メモリにおけるデータ通信量を抑えることができる参照の局所性が高い計算方法を用いる必要がある。
【0090】
本発明の実施形態にかかる数値計算プログラムは、複数の演算器と階層構造のメモリを採用している並列計算機に対して第１時刻における流体および構造体の物理量から微小時間後の第２時刻における流体および構造体の物理量を連成解析する数値計算の指令であって、第１時刻における流体および構造体の物理量のみから第２時刻における流体および構造体の物理量を計算する陽解法を利用している。したがって、他の演算器が第２時刻における流体および構造体の物理量を計算し終えることを待たずに複数の演算器の各々が第２時刻における流体および構造体の物理量を計算することができる。つまり、いわゆる待ち時間が発生しない。
【0091】
また、本発明の実施形態にかかる数値計算プログラムは、用いる支配方程式の各々は、場の時間微分を、前記場および前記場の一階空間微分ならびに二階空間微分で表現したものであるので、これらを時刻および位置で離散化した場合、先述の式（３）のような形で表現されるので、第１時刻における流体および構造体の物理量のみから第２時刻における流体および構造体の物理量を計算することができるのみならず、離散化した一致のインデックスがメモリのアドレスに対応するので、参照の局所性が高い。つまり、階層構造のメモリにおいて、レジスタやキャッシュを有効に活用することができる。
【0092】
さらに、本発明の実施形態にかかる数値計算プログラムは、離散化した支配方程式がすべて先述の式（３）のような形で表現されるので、並列度が高い。すなわち、離散化した支配方程式を複数の演算器に割り当てて並列処理することで、演算器の処理能力を余すことなく活用することができる。
【0093】
本発明の実施形態にかかる数値計算方法および数値計算プログラムの並列性能の検証では、上記特性が確認される。検証に用いた並列計算機は、九州大学情報基盤研究開発センターのＰＩＭＥＲＧＹ ＣＸ４００であり、先述の図２に示される分散メモリ型の並列計算機である。
【0094】
このような分散メモリ型の並列計算機を用いて、本発明の実施形態にかかる数値計算プログラムで処理を行う数値計算問題を固定して、処理に用いる演算器の個数を増やしながら処理能力の測定を行う、いわゆる強スケーリングでの並列性能を測定する。図１０は、並列性能の検証結果を示すグラフである。
【0095】
図１０に示されるグラフから解るように、本発明の実施形態にかかる数値計算方法および数値計算プログラムでは、演算器の個数に比例して処理能力が理想的に増加している。なお、本発明の実施形態にかかる数値計算方法および数値計算プログラムでは、２３６１６並列を用いた場合に１６００億自由度の計算をすることができる。
【0096】
したがって、本発明の実施形態にかかる数値計算プログラムは、図２に示されるような分散メモリ型の並列計算機に対する指令として好適であり、本発明の実施形態にかかる数値計算方法は、図２に示されるような分散メモリ型の並列計算機における処理として好適である。
【0097】
さらに、本発明の実施形態にかかる数値計算方法および数値計算プログラムのＧＰＧＰＵ（Ｇｅｎｅｒａｌ－ｐｕｒｐｏｓｅ ｃｏｍｐｕｔｉｎｇ ｏｎ ｇｒａｐｈｉｃｓ ｐｒｏｃｅｓｓｉｎｇ ｕｎｉｔｓ）への適合性を検証する。
【0098】
ＧＰＧＰＵとは、ＧＰＵは画像処理のための演算装置であるが、これを画像処理以外の目的にも応用する技術である。ＧＰＵは、単体では計算処理を行うことはできないが、別途のＣＰＵからの指令によって数値計算をすることも可能である。このとき、ＧＰＵは、実質的に図１に示されるような共有メモリ型の並列計算機として機能する。しかも、ＧＰＵは、ＣＰＵよりも演算器の個数が多く、メモリバンド幅も大きい。そこで、ＧＰＧＰＵが並列計算機の代替として活用されることも多い。
【0099】
本検証で使用するＧＰＵはＮＶＩＤＩＡ ＴＥＳＬＡ Ｐ１００ ｆｏｒ ＮＶＬｉｎｋ―Ｏｐｔｉｍｉｚｅｄ Ｓｅｒｖｅｒｓであり、Ｂｙｔｅ／Ｆｌｏｐは０．１３８（Ｂ／Ｆ）＝０．７３２（ＴＢｙｔｅ／ｓ）／５．３（ＴＦｌｏｐ／ｓ）である。したがって、疎行列・ベクトル積（ＳｐＭＶ）が１０（Ｂ／Ｆ）であることを考慮すると、ＧＰＵ一枚あたりのＳｐＭＶの計算速度は、高々７３．２（ＧＦｌｏｐ／ｓ）＝０．７３２（ＴＢｙｔｅ）／１０（Ｂｙｔｅ／Ｆｌｏｐ）にしかならない。
【0100】
一方、本発明の実施形態にかかる数値計算方法および数値計算プログラムを適用したＧＰＧＰＵでは、計量テンソルｇのＧＰＵ一枚あたりの計算速度は３３７（ＧＦｌｏｐ／ｓ）である。このことは、以下のことを意味する。
【0101】
先述したように、本発明の実施形態にかかる数値計算方法および数値計算プログラムでは、連立方程式を解く必要がない。一方、比較的計算負荷が高い計量テンソルｇの計算をする必要がある。連立方程式を共役勾配法（反復法）で解くと、ＳｐＭＶが全体の計算における殆どを占める。ＳｐＭＶは、上述したように、Ｂ／Ｆ値が高く、ＧＰＵの演算性能を活かせない。したがって、本発明の実施形態にかかる数値計算方法および数値計算プログラムを適用したＧＰＧＰＵは、連立方程式を解く必要がある従来の処理よりもＧＰＵの性能を有効活用することができ、計算速度が向上している。
【0102】
以上のように、本発明の実施形態にかかる数値計算プログラムは、図１に示されるような共有メモリ型の並列計算機に対する指令として好適であり、本発明の実施形態にかかる数値計算方法は、図１に示されるような共有メモリ型の並列計算機における処理として好適である。また、図２に示されるような分散メモリ型の並列計算機における各ノードが共有メモリ型の並列計算機として構成されている、共有分散メモリ型の並列計算機に対しても、本発明の実施形態にかかる数値計算方法および数値計算プログラムが有効であることは言うまでもない。
【0103】
以上、図面を参照しながら本発明を実施形態に基づいて説明してきたが、本発明は上記の実施形態よって限定されるものではない。上記説明した数値計算方法は、適切に数値計算プログラムとして実装することができ、逆に上記説明した数値計算プログラムの実行は、数値計算方法の実施とみなし得る。また、上記説明した数値計算プログラムは、コンピュータが読み取り可能な媒体に記録された状態で生産等の実施を行うことができる。
【符号の説明】
【0104】
１００，２００ 並列計算機
１０１～１０ｎ 演算器
１１ 階層型の記憶装置
１２ 主記憶装置
１３１～１３ｎ レジスタ
１４１～１４ｎ Ｌ１キャッシュメモリ
１５１～１５ｍ Ｌ２キャッシュメモリ
１６１～１６ｎ ノード
１７ ネットワーク
２０ 流体軸受
２１ 軸
２２ 軸受
２３ 潤滑油
１０１～１０ｎ 演算器
１１ 階層型の記憶装置
１２ 主記憶装置
１３１～１３ｎ レジスタ
１４１～１４ｎ Ｌ１キャッシュメモリ
１５１～１５ｍ Ｌ２キャッシュメモリ
１６１～１６ｎ ノード
１７ ネットワーク
２０ 流体軸受
２１ 軸
２２ 軸受
２３ 潤滑油
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】