JP7221062B2 - 流体解析システム、流体解析方法、および流体解析プログラム - Google Patents

Description

本開示の一側面は、流体解析システム、流体解析方法、および流体解析プログラムに関する。

液体の挙動をコンピュータで解析する手法が知られている。例えば、特許文献１には、ＭＰＳ法に基づく計算アルゴリズムを実行して非圧縮性流体に係る流体の挙動解析を行う方法が記載されている。

特開２００９－１２９１９３号公報

ＭＰＳ法は、連続体に関する方程式を数値的に解くための離散化手法である粒子法の一種である。格子法とは異なり、粒子法は流体の挙動を格子に丸め込まないので、飛沫などのような流体の細部の挙動を予測することができ、また、流体の大きな変形を予測することも可能である。その一方で、粒子法では計算の精度を上げようとすると計算時間が増大してしまう。そのため、粒子法を用いた液体の挙動の予測を短時間で実行することが望まれている。

本開示の一側面に係る流体解析システムは少なくとも一つのプロセッサを備える。少なくとも一つのプロセッサは、液体を構成する複数の粒子のそれぞれの第１時点における位置および速度に基づく入力ベクトルを取得し、入力ベクトルを機械学習モデルに入力することで、第１時点から所与の時間ステップが経過した第２時点における、複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成する。

本開示の一側面に係る流体解析方法は、少なくとも一つのプロセッサを備える流体解析システムにより実行される。この流体解析方法は、液体を構成する複数の粒子のそれぞれの第１時点における位置および速度に基づく入力ベクトルを取得するステップと、入力ベクトルを機械学習モデルに入力することで、第１時点から所与の時間ステップが経過した第２時点における、複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成するステップとを含む。

本開示の一側面に係る流体解析プログラムは、液体を構成する複数の粒子のそれぞれの第１時点における位置および速度に基づく入力ベクトルを取得するステップと、入力ベクトルを機械学習モデルに入力することで、第１時点から所与の時間ステップが経過した第２時点における、複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成するステップとをコンピュータに実行させる。

このような側面においては、液体を構成する個々の粒子の位置および速度が機械学習を用いて計算される。粒子法を用いた予測の少なくとも一部に機械学習を導入することで、その予測を短時間で実行することが可能になる。

本開示の一側面によれば、粒子法を用いた液体の挙動の予測を短時間で実行することができる。

ＭＰＳ法の手順の一例を示すフローチャートである。 実施形態に係る流体解析システムで用いられるコンピュータのハードウェア構成の一例を示す図である。 実施形態に係る流体解析システムの機能構成の一例を示す図である。 実施形態に係る流体解析システムによる学習処理の一例を示すフローチャートである。 実施形態に係る流体解析システムによる予測処理の一例を示すフローチャートである。 実施形態に係る流体解析システムによる学習処理の別の例を示すフローチャートである。 実施形態に係る流体解析システムによる予測処理の別の例を示すフローチャートである。 実施形態に係る流体解析システムによる学習処理の別の一例を示すフローチャートである。 実施形態に係る流体解析システムによる予測処理の別の例を示すフローチャートである。

以下、添付図面を参照しながら本開示での実施形態を詳細に説明する。なお、図面の説明において同一または同等の要素には同一の符号を付し、重複する説明を省略する。

［システムの概要］
実施形態に係る流体解析システムは液体の挙動を予測するコンピュータシステムである。より具体的には、流体解析システムは何らかの構造物に関連する液体の挙動を予測する。液体の挙動とは、時間の経過に伴う液体の動きのことをいう。処理対象となる液体の種類、すなわち、モデル化される液体（液体モデル）の種類は何ら限定されず、例えば水、オイルなどでもよい。処理対象となる液体を取り巻く雰囲気（例えば温度）は限定されず、例えば、高温下で液化した物質が処理対象であってもよい。液体に関連する構造物の種類、すなわち計算の前提として設定される構造物の種類も何ら限定されず、人工物でも自然物でもよい。構造物の例として水車、タンク、水洗トイレ、船舶、ギアボックス、堤防、ダムなどが挙げられるが、当然にこれらに限定されない。

流体解析システムは機械学習を利用する。機械学習とは、与えられた情報に基づいて反復的に学習することで法則またはルールを自律的に見つけ出す手法である。入力データから（連続的な）数値を予測する回帰問題を扱うことができる限り、機械学習の具体的な手法は限定されない。例えば、流体解析システムはランダムフォレスト、勾配ブースティング決定木（ＧＢＤＴ）、ニューラルネットワーク、サポートベクター回帰（ＳＶＲ）などの任意の手法が用いられてよい。

流体解析システムは、データを学習することで計算モデル（機械学習モデル）を生成し、この計算モデルを学習済みモデルとして取得することができる。これは学習フェーズに相当する。学習済みモデルは、液体の挙動を予測するために最適であると推定される計算モデル（機械学習モデル）であり、“現実に最適である計算モデル（機械学習モデル）”とは限らないことに留意されたい。流体解析システムは、学習済みモデルを用いて入力データを処理することで、液体の挙動を示す予測データを出力することもできる。これは予測フェーズに相当する。

学習済みモデルはコンピュータシステム間で移植可能である。したがって、或るコンピュータシステムで生成された学習済みモデルを、別のコンピュータシステムで用いることができる。もちろん、一つのコンピュータシステムが学習済みモデルの生成および利用の双方を実行してもよい。すなわち、流体解析システムは、学習フェーズおよび予測フェーズの双方を実行してもよいし、学習フェーズおよび予測フェーズのいずれか一方を実行しなくてもよい。

［粒子法］
流体解析システムは液体の挙動を予測するために粒子法を利用する。粒子法は、連続体に関する方程式を数値的に解くための離散化手法の一つである。格子法とは異なり、粒子法は流体の挙動を格子に丸め込まないので、飛沫などのような流体の細部の挙動を予測することに適しており、また、流体の大きな変形を予測することにも適している。粒子法では、流体を複数の粒子に分割することでモデル化する。この粒子は、分子、液滴などの実体そのものではなく、液体の動きと共に移動する計算点である。個々の粒子には何らかの物理量（例えば、速度、質量、圧力など）が関連付けられる。

粒子法の代表的なものとしてＤＥＭ（Ｄｉｓｃｒｅｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ）法、ＳＰＨ（Ｓｍｏｏｔｈｅｄ Ｐａｒｔｉｃｌｅ Ｈｙｄｒｏｄｙｎａｍｉｃｓ）法、およびＭＰＳ（Ｍｏｖｉｎｇ Ｐａｒｔｉｃｌｅ Ｓｅｍｉ－ｉｍｐｌｉｃｉｔ）法がある。ＤＥＭ法は、粒子間にバネを仮定して粒子の衝突および反発を計算する手法である。ＳＰＨ法は、流体解析の物理方程式であるナビエ・ストークス方程式を計算するために、近傍粒子上の関数値の重み付き平均によって補間関数を定義する手法である。ＭＰＳ法は、そのような補間関数を定義せず、各粒子の位置、速度などについて、重み関数で補正した近傍粒子点の和でナビエ・ストークス方程式を計算する手法である。一般に、ＤＥＭ法は粉体の問題を扱うのに対して、ＳＰＨ法およびＭＰＳ法は液体の問題を扱う。ＳＰＨ法は陽解法（ｆｕｌｌｙ ｅｘｐｌｉｃｉｔ）で計算を行うため、計算ステップを細かく刻む必要があるのに対して、ＭＰＳ法は半陰的解法（ｓｅｍｉ－ｉｍｐｌｉｃｉｔ）なので、計算ステップを大きく刻めるという利点がある。

本実施形態ではＭＰＳ法を応用する。流体解析システムはＭＰＳ法を利用した予測において機械学習モデルを用いることで、計算時間の短縮を可能にする。

流体解析システムの詳細に先立ち、図１を参照しながらＭＰＳ法の手順を処理フローＳ１として説明する。図１はＭＰＳ法の手順の一例を示すフローチャートである。

ステップＳ１１では、粒子数密度などの計算条件が設定される。ステップＳ１２では個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）が設定される。ｒ^０は個々の粒子の初期位置を示し、ｕ^０は個々の粒子の初期速度を示し、Ｐ^０は個々の粒子間の初期圧力を示す。ただしＰ^０は必須ではなく、Ｐ^０が与えられなくても流体解析は可能である。本開示では、粒子の位置および速度は２次元または３次元のベクトルで表され、粒子に作用する圧力はスカラーで表される。上記の通り、粒子とはシミュレーションにおける計算点であり、したがって、粒子の大きさは、解析対象の大きさ、期待する解像度などにより任意に設定されてよい。

ステップＳ１３では個々の粒子の仮速度が計算される。仮速度は、現在の時点から時間ステップΔｔだけ経過した時点における粒子の速度を示す仮の値である。時間ステップは、解析において１回あたりに進める時間間隔である。

粒子の仮速度ｕ^＊は、粘性項および重力・外力項を含む式（１）により得られる。ｇは重力加速度である。したがって、粒子の仮速度ｕ^＊は、該粒子に作用する力に基づいて算出されるということができる。

粘性項は、重み関数ｗ（・）を含む下記式（２）で得られる。

ここで、

である。νは動粘性係数を示し、ｄは計算領域の次元を示す。ｎ^０は一粒子の近傍粒子の粒子数密度の基準値を示す。λ^０は一粒子の影響の範囲内にある近傍粒子との距離の２乗の重み平均値を示す。ｕ_ｉ ^ｋは計算ステップｋにおけるｉ番目の粒子の速度を示す（ｉは粒子番号を示す）。ｒ_ｉ ^ｋは計算ステップｋにおけるｉ番目の粒子の位置を示す。なお、粘性項は式（２）に限定されず、高精度粒子法に基づく別の手法によって計算されてもよい。

重み関数ｗ（・）については一般に式（３）が用いられる。

ここで、｜ｒ_ｊ－ｒ_ｉ｜は二つの粒子ｉ，ｊ間の距離を示す。ｒ_ｅは或る粒子を近傍粒子と見做す距離（近傍粒子の距離）の基準値である。なお、重み関数は式（３）に限定されるものではなく、他の種類の重み関数が用いられてもよい。

ステップＳ１４では個々の粒子の仮位置ｒ^＊が式（４）によって計算される。仮位置は、現在の時点から時間ステップだけ経過した時点における粒子の位置を示す仮の値である。上述したように粒子の仮速度ｕ^＊は該粒子に作用する力に基づいて算出されるので、粒子の仮位置ｒ^＊もその力に基づいて算出されるということができる。

ステップＳ１５では、粒子の仮位置での圧力Ｐ^ｋ＋１が式（５）によって計算される。

ここで、行列Ａの要素ａ_ｉｊは式（６）で示される。

ベクトルｂの要素ｂ_ｉは式（７）で示される。

ρ^０は液体の密度の基準値を示す。ｎ_ｉ ^＊はｉ番目の粒子の仮位置ｒ_ｉ ^＊での粒子密度であり次式で定義される。

なお、行列Ａおよびベクトルｂの計算方法は上記のものに限定されず、他の手法が用いられてもよい。

ステップＳ１６では、圧力勾配∇Ｐ^ｋ＋１による速度変化量ｕ´が式（８）によって計算される。

圧力勾配∇Ｐ^ｋ＋１は式（９）で示される。

ここで、Ｐ＾_ｉ ^ｋ＋１は、ｉ番目の粒子と近傍粒子との圧力差が負圧にならないように補正するための値である。圧力勾配の計算方法は式（９）に限定されず、他の手法が用いられてもよい。

ステップＳ１７では、粒子の速度および位置が更新される。これは、仮に求めた粒子の速度および位置を確定することを意味する。粒子の速度ｕ^ｋ＋１は式（１０）で示され、更新後の粒子の位置ｒ^ｋ＋１は式（１１）で示される。

ステップＳ１８に示すように、ステップＳ１３～Ｓ１７の処理は、予測における所与の最終時刻に到達するまで繰り返し実行される。その最終時刻にまだ到達していない場合には（ステップＳ１８においてＮＯ）、ステップＳ１９において時間ステップΔｔの分だけ時間が進められて、その後、更新後の粒子の速度および位置に対してステップＳ１３～Ｓ１７の処理が実行される。ステップＳ１８において終了条件が満たされるまで個々の粒子の速度および位置が更新され続けることで、所与の初期時刻から所与の最終時刻までの間の液体の挙動が予測される。

ＭＰＳ法による処理（処理フローＳ１）の中で計算時間が特に掛かるのが仮位置での圧力の計算（ステップＳ１５）である。この処理では圧力Ｐ^ｋ＋１を連立一次方程式で解く必要があり、そのためには行列Ａおよびベクトルｂを用いて反復法もしくは直接法で計算する必要がある。そのため、複数の計算ノードによる分散メモリ型並列計算ではノード間のデータ通信量が大きな負荷となり、高速化が困難である。一方、ステップＳ１５以外の処理は、複数の計算ノードによる並列計算に適している。例えば、液体モデルが存在する空間であるシミュレーション空間を複数の区画に分割して個々の計算ノードが個々の区画について計算を実行することで、計算の高速化が可能である。本実施形態では、少なくとも仮位置での圧力の計算（ステップＳ１５）を含む処理を機械学習で実行することで、液体の挙動を予測する処理の高速化を図る。

［システムの構成］
図２は実施形態に係る流体解析システム１０を構成するコンピュータ１００の一般的なハードウェア構成の一例を示す。例えば、コンピュータ１００はプロセッサ１０１、主記憶部１０２、補助記憶部１０３、通信制御部１０４、入力装置１０５、および出力装置１０６を備える。プロセッサ１０１はオペレーティングシステムおよびアプリケーション・プログラムを実行する。主記憶部１０２は例えばＲＯＭおよびＲＡＭで構成される。補助記憶部１０３は例えばハードディスクまたはフラッシュメモリで構成され、一般に主記憶部１０２よりも大量のデータを記憶する。通信制御部１０４は例えばネットワークカードまたは無線通信モジュールで構成される。入力装置１０５は例えばキーボード、マウス、タッチパネルなどで構成される。出力装置１０６は例えばモニタおよびスピーカで構成される。

流体解析システム１０の各機能要素は、補助記憶部１０３に予め記憶される流体解析プログラム１１０により実現される。具体的には、各機能要素は、プロセッサ１０１または主記憶部１０２の上に流体解析プログラム１１０を読み込ませてその流体解析プログラム１１０を実行させることで実現される。プロセッサ１０１はその流体解析プログラム１１０に従って、通信制御部１０４、入力装置１０５、または出力装置１０６を動作させ、主記憶部１０２または補助記憶部１０３におけるデータの読み出しおよび書き込みを行う。処理に必要なデータまたはデータベースは主記憶部１０２または補助記憶部１０３内に格納される。

流体解析プログラム１１０は、例えば、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、半導体メモリなどの有形の記録媒体に固定的に記録された上で提供されてもよい。あるいは、流体解析プログラム１１０は、搬送波に重畳されたデータ信号として通信ネットワークを介して提供されてもよい。

流体解析システム１０は１台のコンピュータ１００で構成されてもよいし、複数台のコンピュータ１００で構成されてもよい。複数台のコンピュータ１００を用いる場合には、これらのコンピュータ１００がインターネットやイントラネットなどの通信ネットワークを介して接続されることで、論理的に一つの流体解析システム１０が構築される。

図３は実施形態に係る流体解析システム１０の機能構成の一例を示す図である。本実施形態では、流体解析システム１０は機能要素として学習部１１および予測部１２を備える。学習部１１は、機械学習を実行することで学習済みモデルを生成および取得する機能要素である。予測部１２は、その学習済みモデルを用いて入力データから予測データを求める機能要素である。

［システムの動作］
流体解析システム１０の動作を説明するとともに本実施形態に係る流体解析方法について説明する。まず、図４および図５を参照しながら、ＭＰＳ法のすべてを機械学習に置き換える手法の一例を説明する。図４は流体解析システム１０の学習処理の一例を処理フローＳ２として示すフローチャートである。図５は流体解析システム１０の予測処理の一例を処理フローＳ３として示すフローチャートである。図４および図５はそれぞれ、学習フェーズおよび予測フェーズに対応する。処理フローＳ２，Ｓ３のいずれも、本開示に係る流体解析方法に対応し得る。処理フローＳ２は学習済みモデルの生成方法の一例であるともいえる。

図４を参照しながら、学習部１１により実行される処理フローＳ２を説明する。学習部１１は、ステップＳ２１では一つの液体モデルを示す教師データに基づいて計算条件を取得し、ステップＳ２２ではその教師データに基づいて個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）を取得する。ステップＳ２１，Ｓ２２はそれぞれステップＳ１１，Ｓ１２に対応する。

教師データは、初期時刻（粒子の初期条件に対応する時刻）から所与の最終時刻までの間の個々の時点（すなわち、時間ステップΔｔずつ進む個々の時点）における、全粒子のそれぞれの位置および速度を示すデータである。機械学習で得られる出力ベクトルが各粒子の圧力をさらに含む場合には、教師データは各時点における全粒子のそれぞれの圧力をさらに含む。教師データの準備方法は限定されず、例えば、教師データは、別途実行されるコンピュータ・シミュレーションによって生成されてもよい。一般に、複数の液体モデルのそれぞれについて教師データが用意される。用意される教師データの個数（液体モデルの個数）は限定されない。

教師データの取得方法は限定されない。例えば、学習部１１は流体解析システム１０のユーザ操作により入力されたデータを受け付けてもよいし、所与のデータベースにアクセスしてデータを読み出してもよいし、他のコンピュータシステムから送られてきたデータを受信してもよい。

ステップＳ２３では、学習部１１は、複数の粒子の位置および速度を示す入力ベクトルを用いて機械学習を実行する。具体的には、学習部１１は、第１計算ステップｋにおける複数の粒子の位置および速度を示す入力ベクトルを教師データから生成する。そして、学習部１１はその入力ベクトルを機械学習モデルに入力し、第２計算ステップｋ＋１における該複数の粒子の位置および速度を示す出力ベクトルを得る。ステップＳ２３はステップＳ１３～Ｓ１９に対応する。第１計算ステップは第１時点に対応し、第２計算ステップは第２時点に対応する。第１時点は或る時刻ｔであり、第２時点は該時刻ｔから時間ステップΔｔだけ進んだ時刻（ｔ＋Δｔ）である。第１時点の初期値は粒子の初期条件に対応する時点であり、その後、計算ステップが１ずつ進むごとに第１時点が時間ステップΔｔずつ進む。機械学習モデルに入力される複数の粒子は、教師データで示される全粒子であってもよいし、一部の粒子のみであってもよい。一部の粒子は例えばランダムに選択されてもよい。一部の粒子のみを用いることで学習時間を短縮することができる。

粒子の総数をＮとし、粒子の位置をｒとし、粒子の速度をｕとし、さらに、時間ステップ毎に進む時刻を計算ステップｋで表すとする。この場合、入力ベクトルは｛（ｒ^ｋ _１，ｕ^ｋ _１），（ｒ^ｋ _２，ｕ^ｋ _２），…，（ｒ^ｋ _Ｎ，ｕ^ｋ _Ｎ）｝と表され、出力ベクトルは｛（ｒ^ｋ＋１ _１，ｕ^ｋ＋１ _１），（ｒ^ｋ＋１ _２，ｕ^ｋ＋１ _２），…，（ｒ^ｋ＋１ _Ｎ，ｕ^ｋ＋１ _Ｎ）｝と表される。第２時点における粒子の位置での圧力を求めてもよく、この場合には、出力ベクトルは｛（ｒ^ｋ＋１ _１，ｕ^ｋ＋１ _１，Ｐ^ｋ＋１ _１），（ｒ^ｋ＋１ _２，ｕ^ｋ＋１ _２，Ｐ^ｋ＋１ _２），…，（ｒ^ｋ＋１ _Ｎ，ｕ^ｋ＋１ _Ｎ，Ｐ^ｋ＋１ _Ｎ）｝と表される。

この機械学習では、教師データで示される第２時点での値が正解として用いられる。学習部１１は、入力ベクトルを機械学習モデルに与えることで出力ベクトルを生成する。そして、学習部１１は予測結果を示す出力ベクトルと正解との差（すなわち、誤差）に基づいて、バックプロパゲーション（誤差逆伝播法）などの手法を用いて機械学習モデル内のパラメータを更新する。なお、機械学習モデル内の更新されるパラメータの種類は限定されず、例えば、様々な設定値または制限値であってもよい。ニューラルネットワークの重みはパラメータの一例である。

ステップＳ２３において、学習部１１は機械学習を一つの計算ステップごとに実行してもよいし、すべての計算ステップについて一度に実行してもよい。すべての計算ステップについて一度に実行する場合には、入力ベクトルは｛（ｒ^１ _１，ｕ^１ _１），（ｒ^１ _２，ｕ^１ _２），…，（ｒ^１ _Ｎ，ｕ^１ _Ｎ），（ｒ^２ _１，ｕ^２ _１），（ｒ^２ _２，ｕ^２ _２），…，（ｒ^２ _Ｎ，ｕ^２ _Ｎ），…，（ｒ^ｋ _１，ｕ^ｋ _１），（ｒ^ｋ _２，ｕ^ｋ _２），…，（ｒ^ｋ _Ｎ，ｕ^ｋ _Ｎ）｝と表される。また、出力ベクトルは｛（ｒ^２ _１，ｕ^２ _１），（ｒ^２ _２，ｕ^２ _２），…，（ｒ^２ _Ｎ，ｕ^２ _Ｎ），（ｒ^３ _１，ｕ^３ _１），（ｒ^３ _２，ｕ^３ _２），…，（ｒ^３ _Ｎ，ｕ^３ _Ｎ），…，（ｒ^ｋ＋１ _１，ｕ^ｋ＋１ _１），（ｒ^ｋ＋１ _２，ｕ^ｋ＋１ _２），…，（ｒ^ｋ＋１ _Ｎ，ｕ^ｋ＋１ _Ｎ）｝と表される。圧力を含む出力ベクトルも同様に表現できる。

一つの教師データを用いた機械学習の終了条件は任意に設定されてよい。例えば、終了条件は誤差に基づいて設定されてもよいし、学習のサイクル数に基づいて設定されてもよい。

学習フェーズでは学習部１１は複数の液体モデルのそれぞれの教師データを用いて処理フローＳ２を実行し得る。任意の個数の液体モデルについて処理フローＳ２を実行した後に、学習部１１は処理された機械学習モデルを学習済みモデルとして取得する。この処理は、学習済みモデルが生成されたことを意味する。学習部１１はその学習済みモデルを出力する。学習済みモデルの出力方法は限定されない。例えば、学習部１１は学習済みモデルを、メモリ、データベースなどの記憶装置に格納してもよいし、他のコンピュータシステムに送信してもよい。本実施形態では、予測部１２がその学習済みモデルを用いる。

図５を参照しながら、予測部１２により実行される処理フローＳ３を説明する。予測部１２は、ステップＳ３１では予測しようとする液体モデルの計算条件を取得し、ステップＳ３２では該液体モデルの個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）を取得する。ステップＳ３１，Ｓ３２はそれぞれステップＳ１１，Ｓ１２に対応する。計算条件と粒子の初期条件とを示す入力データの取得方法は限定されない。例えば、予測部１２は流体解析システム１０のユーザ操作により入力されたデータを受け付けてもよいし、所与のデータベースにアクセスしてデータを読み出してもよいし、他のコンピュータシステムから送られてきたデータを受信してもよい。

ステップＳ３３では、予測部１２は、複数の粒子の位置および速度を示す入力ベクトルを学習済みモデルに入力することで出力ベクトルを生成する。具体的には、予測部１２は、第１時点における全粒子の位置および速度を示す入力ベクトルを学習済みモデルに入力し、第２時点における全粒子の位置および速度を示す出力ベクトルを得る。ステップＳ３３はステップＳ１３～Ｓ１７に対応する。第１時点は或る時刻ｔであり、第２時点は該時刻ｔから時間ステップΔｔだけ進んだ時刻（ｔ＋Δｔ）である。第１時点の初期値は粒子の初期条件に対応する時点であり、その後、第１時点は時間ステップΔｔずつ進む。入力ベクトルは、１ループ目では初期条件に対応する全粒子の位置および速度を含んで構成され、２ループ目以降では、前のループで得られた出力ベクトルで示される全粒子の位置および速度を含んで構成される。処理フローＳ２に対応して、出力ベクトルは圧力を含んでもよいし含まなくてもよい。

ステップＳ３４では、予測部１２は計算を終了する否かを判定する。具体的には、予測部１２は予測における所与の最終時刻に到達したか否かを判定する。第２時刻が最終時刻でない場合には（ステップＳ３４においてＮＯ）、処理はステップＳ３５に進む。ステップＳ３５では、予測部１２は時間ステップΔｔだけ時間を進める。そして、予測部１２はその進めた時間を第１時点として、ステップＳ３３の処理を再び実行する。

計算を終了すると判定した場合、すなわち所与の最終時刻に到達した場合には、処理はステップＳ３６に進む。ステップＳ３６では、予測部１２は予測データを出力する。予測データの構造は限定されず、例えば、予測部１２は、初期時刻（粒子の初期条件に対応する時刻）から所与の最終時刻までの各時点における、全粒子のそれぞれの位置および速度を示す予測データを出力してもよい。あるいは、予測部１２は初期時刻から最終時刻までの間の任意の一または複数の時点における、全粒子のそれぞれの位置および速度を示す予測データを出力してもよい。予測データの出力方法は限定されない。例えば、予測部１２は予測データを、モニタ上に表示してもよいし、所定のデータベースに格納してもよいし、他のコンピュータシステムに送信してもよい。予測データの表現方法も任意であり、例えば、予測データはコンピュータグラフィックス、表、数字列などのなどの任意の手法で表現されてよい。いずれにしても、流体解析システム１０のユーザはその予測データを参照することで液体の挙動の予測を得ることができる。

次に、図６および図７を参照しながら、ＭＰＳ法の一部を機械学習に置き換える手法の一例を説明する。図６は流体解析システム１０の学習処理の一例を処理フローＳ４として示すフローチャートである。図７は流体解析システム１０の予測処理の一例を処理フローＳ５として示すフローチャートである。図６および図７はそれぞれ、学習フェーズおよび予測フェーズに対応する。処理フローＳ４，Ｓ５のいずれも、本開示に係る流体解析方法に対応し得る。処理フローＳ４は学習済みモデルの生成方法の一例であるともいえる。

図６を参照しながら、学習部１１により実行される処理フローＳ４を説明する。学習部１１は、ステップＳ４１では一つの液体モデルを示す教師データに基づいて計算条件を取得し、ステップＳ４２ではその教師データに基づいて個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）を取得する。処理フローＳ２と同様に教師データの取得方法は限定されない。ステップＳ４１，Ｓ４２はそれぞれステップＳ１１，Ｓ１２に対応する。

学習部１１は、ステップＳ４３では個々の粒子の仮速度を求め、ステップＳ４４では個々の粒子の仮位置を求め、ステップＳ４５では個々の粒子について仮位置での圧力を求める。ステップＳ４３，Ｓ４４，Ｓ４５はそれぞれステップＳ１３，Ｓ１４，Ｓ１５に対応する。ステップＳ４３，Ｓ４４は、第１時点（時刻ｔ）における個々の粒子の位置および速度に基づいて、第２時点（時刻ｔ＋Δｔ）における個々の粒子の位置および速度を仮に求める処理である。処理フローＳ４では、ステップＳ４５で得られる圧力Ｐ^ｋ＋１が教師データとなり、正解として用いられる。

ステップＳ４６では、学習部１１は個々の粒子について特徴量を求める。本実施形態では、特徴量は以下に示す第１特徴量～第９特徴量を一つにまとめたベクトルとして表現される。この例では、学習部１１は第１特徴量～第９特徴量を以下の方針（Ａ）～（Ｄ）に基づいて計算する。

（Ａ）近傍粒子との距離に関連した値を用いる。
（Ｂ）圧力Ｐ^ｋ＋１を求める式（５）に関連した値を用いる。
（Ｃ）粒子密度に関連した値を用いる。
（Ｄ）相対密度に関連した値を用いる。

学習部１１は第１特徴量を計算する。第１特徴量は粒子ｉから距離ｒ_ｅ内に存在するｎ^０個の近傍粒子ｊとの距離である。第１特徴量はｎ^０次元のベクトルであり、各要素は昇順に配置される。第１特徴量のｎ^０個の要素Ｆ１_ｊは次式で定義される。

ここで、｜ｒ_ｊ ^＊－ｒ_ｉ ^＊｜は二つの粒子ｉ，ｊ間の距離（二つの仮位置の間の距離）を示す。ｒ_ｅは或る粒子を近傍粒子と見做す距離（近傍粒子の距離）の基準値である。距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第１特徴量の要素に設定される。

学習部１１は第２特徴量を計算する。第２特徴量は式（５）の行列Ａに関連する特徴量であり、ｎ^０次元のベクトルである。第２特徴量のｎ^０個の要素Ｆ２_ｊは次式で定義される。

ここで、ａ_ｉｉは式（６）で得られる。距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第２特徴量の要素に設定される。

学習部１１は第３特徴量を計算する。第３特徴量は式（５）の行列Ａに関連する特徴量であり、ｎ^０次元のベクトルである。第３特徴量のｎ^０個の要素Ｆ３_ｊは次式で定義される。

ここで、ｎ_ｊ ^＊はｊ番目の粒子の仮位置ｒ_ｊ ^＊での粒子密度である。距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第３特徴量の要素に設定される。

学習部１１は第４特徴量を計算する。第４特徴量は式（５）の行列Ａに関連する特徴量であり、ｎ^０次元のベクトルである。第４特徴量のｎ^０個の要素Ｆ４_ｊは次式で定義される。

ここで、ｎ_ｊ ^＊＊はｊ番目の粒子の仮位置ｒ_ｊ ^＊での相対密度であり、次式で定義される。

距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第４特徴量の要素に設定される。

学習部１１は第５特徴量を計算する。第５特徴量は式（５）のベクトルｂに関連する特徴量であり、ｎ^０次元のベクトルである。第５特徴量のｎ^０個の要素Ｆ５_ｊは次式で定義される。

距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第５特徴量の要素に設定される。

学習部１１は第６特徴量を計算する。第６特徴量は上記の粒子密度に関連する特徴量であり、ｎ^０次元のベクトルである。第６特徴量のｎ^０個の要素Ｆ６_ｊは次式で定義される。

距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第６特徴量の要素に設定される。

学習部１１は第７特徴量を計算する。第７特徴量は上記の相対密度に関連する特徴量であり、ｎ^０次元のベクトルである。第７特徴量のｎ^０個の要素Ｆ７_ｊは次式で定義される。

距離ｒ_ｅ内に存在する粒子がｎ^０個に満たない場合は、ゼロが代替として第７特徴量の要素に設定される。

学習部１１は第８特徴量を計算する。第８特徴量は式（５）のベクトルｂに関連する特徴量であり、スカラーである。第８特徴量は次式で定義される。

学習部１１は第９特徴量を計算する。第９特徴量は上記の粒子密度に関連する特徴量であり、スカラーである。第９特徴量は次式で定義される。

ステップＳ４７では、学習部１１は複数の粒子の特徴量を示す入力ベクトルを用いて機械学習を実行する。

この機械学習に対応する理論は次の通りである。計算ステップｋにおいて粒子ｉの圧力Ｐ_ｉ ^ｋ＋１は次式により得られる。Ｒｅｇｒｅｓｓｉｏｎ（）は回帰式を表し、Ｆ_ｉは粒子ｉの特徴量を示す。

第１特徴量～第７特徴量はそれぞれｎ^０次元のベクトルなので、それぞれの説明変数はｎ^０個である。第８特徴量および第９特徴量はスカラーなので、それぞれの説明変数は１個である。したがって、機械学習での説明変数は（７×ｎ^０＋２）個である。

学習部１１は、複数の粒子のそれぞれの特徴量を示す入力ベクトルを用いて機械学習を実行する。具体的には、学習部１１は、複数の粒子の特徴量を示す入力ベクトルを機械学習モデルに入力し、該複数の粒子の圧力Ｐ^ｋ＋１を示す出力ベクトルを得る。粒子の特徴量は、該粒子の第１時点における位置および速度に基づいて得られる。したがって、複数の粒子のそれぞれの特徴量を示す入力ベクトルも、複数の粒子のそれぞれの第１時点における位置および速度に基づく。処理フローＳ２と同様に、機械学習モデルに入力される複数の粒子は、教師データで示される全粒子であってもよいし、ランダム選択などの手法により選択された一部の粒子のみであってもよい。

粒子の総数をＮとし、粒子の特徴量をＦとし、さらに、時間ステップ毎に進む時刻を計算ステップｋで表すとする。この場合、入力ベクトルは｛Ｆ_１ ^ｋ，Ｆ_２ ^ｋ，…，Ｆ_Ｎ ^ｋ｝と表され、出力ベクトルは｛Ｐ_１ ^ｋ＋１，Ｐ_２ ^ｋ＋１，…，Ｐ_Ｎ ^ｋ＋１｝と表される。

この機械学習では、ステップＳ４５で得られた圧力Ｐ^ｋ＋１が正解として用いられる。学習部１１は出力ベクトルと正解との差（すなわち、誤差）に基づいて機械学習モデル内のパラメータを更新する。

ステップＳ４８では、学習部１１は、圧力勾配∇Ｐ^ｋ＋１による速度変化量ｕ´を式（８）によって計算する。ステップＳ４８はステップＳ１６に対応する。

ステップＳ４９では、学習部１１は粒子の速度および位置を更新する。これは、仮に求めた粒子の速度および位置を確定することを意味する。粒子の速度ｕ^ｋ＋１は式（１０）で示され、更新後の粒子の位置ｒ^ｋ＋１は式（１１）で示される。ステップＳ４９はステップＳ１７に対応する。

ステップＳ４３～Ｓ４９で示される一連の処理は、計算ステップｋ（第１時点）での値を示す入力から、計算ステップｋ＋１（第２時点）での値を示す出力を得る処理である。学習部１１はこの一連の処理を一つの計算ステップごとに実行してもよいし、すべての計算ステップについて一度に実行してもよい。すべての計算ステップについて一度に実行する場合には、入力ベクトルは｛Ｆ_１ ^１，Ｆ_２ ^１，…，Ｆ_Ｎ ^１,Ｆ_１ ^２，Ｆ_２ ^２，…，Ｆ_Ｎ ^２，…,Ｆ_１ ^ｋ，Ｆ_２ ^ｋ，…，Ｆ_Ｎ ^ｋ｝と表され、出力ベクトルは｛Ｐ_１ ^２，Ｐ_２ ^２，…，Ｐ_Ｎ ^２，Ｐ_１ ^３，Ｐ_２ ^３，…，Ｐ_Ｎ ^３，…，Ｐ_１ ^ｋ＋１，Ｐ_２ ^ｋ＋１，…，Ｐ_Ｎ ^ｋ＋１｝と表される。

学習フェーズでは学習部１１は複数の液体モデルのそれぞれの教師データを用いて処理フローＳ４を実行し得る。任意の個数の液体モデルについて処理フローＳ４を実行した後に、学習部１１は処理された機械学習モデルを学習済みモデルとして取得する。この処理は、学習済みモデルが生成されたことを意味する。処理フローＳ２と同様に、学習部１１はその学習済みモデルを任意の手法で出力する。

図７を参照しながら、予測部１２により実行される処理フローＳ５を説明する。予測部１２は、ステップＳ５１では予測しようとする液体モデルの計算条件を取得し、ステップＳ５２では該液体モデルの個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）を取得する。処理フローＳ３と同様に、液体モデルを示す入力データの取得方法は限定されない。ステップＳ５１，Ｓ５２はそれぞれステップＳ１１，Ｓ１２に対応する。

予測部１２は、ステップＳ５３では個々の粒子の仮速度を求め、ステップＳ５４では個々の粒子の仮位置を求める。ステップＳ５３，Ｓ５４はそれぞれステップＳ１３，Ｓ１４に対応する。ステップＳ５３，Ｓ５４は、第１時点（時刻ｔ）における個々の粒子の位置および速度に基づいて、第２時点（時刻ｔ＋Δｔ）における個々の粒子の位置および速度を仮に求める処理である。

ステップＳ５５では、予測部１２は個々の粒子について特徴量を求める。ステップＳ５５はステップＳ４６に対応し、したがって、予測部１２は、ステップＳ４６における学習部１１と同様の手法で特徴量（第１特徴量～第９特徴量を一つにまとめたベクトル）を求める。

ステップＳ５６では、予測部１２は、複数の粒子の特徴量を示す入力ベクトルを学習済みモデルに入力することで出力ベクトルを生成する。具体的には、予測部１２は、各粒子の特徴量を示す入力ベクトルを学習済みモデルに入力し、各粒子の圧力を示す出力ベクトルを得る。ステップＳ５６はステップＳ４７に対応する。粒子の特徴量は、該粒子の第１時点における位置および速度に基づいて得られる。したがって、複数の粒子のそれぞれの特徴量を示す入力ベクトルも、複数の粒子のそれぞれの第１時点における位置および速度に基づく。第１時点の初期値は粒子の初期条件に対応する時点であり、その後、第１時点は時間ステップΔｔずつ進む。１ループ目では、入力ベクトルは、初期条件に対応する全粒子の位置および速度に基づいて算出された全粒子の特徴量で構成される。２ループ目以降では、入力ベクトルは、前のループで得られた出力ベクトルで示される圧力に基づいて算出された全粒子の特徴量で構成される。

ステップＳ５７では、予測部１２は、出力ベクトルＰ^ｋ＋１を用いて、圧力勾配∇Ｐ^ｋ＋１による速度変化量ｕ´を式（８）により求める。ステップＳ５７はステップＳ１６に対応する。

ステップＳ５８では、予測部１２は個々の粒子の速度および位置を更新する。ステップＳ５８はステップＳ１７に対応する。

ステップＳ５９では、予測部１２は計算を終了する否かを判定する。具体的には、予測部１２は予測における所与の最終時刻に到達したか否かを判定する。第２時刻が最終時刻でない場合には（ステップＳ５９においてＮＯ）、処理はステップＳ５９１に進む。ステップＳ５９１では、予測部１２は時間ステップΔｔだけ時間を進める。そして、予測部１２はその進めた時間を第１時点として、ステップＳ５３～Ｓ５８の処理を再び実行する。

計算を終了すると判定した場合、すなわち所与の最終時刻に到達した場合には、処理はステップＳ５９２に進む。ステップＳ５９２では、予測部１２は予測データを出力する。処理フローＳ３と同様に、予測データのデータ構造、出力方法、および表現方法はいずれも限定されない。

次に、図８および図９を参照しながら、ＭＰＳ法の一部を機械学習に置き換える手法の別の例を説明する。図８は流体解析システム１０の学習処理の一例を処理フローＳ６として示すフローチャートである。図９は流体解析システム１０の予測処理の一例を処理フローＳ７として示すフローチャートである。図８および図９はそれぞれ、学習フェーズおよび予測フェーズに対応する。処理フローＳ６，Ｓ７のいずれも、本開示に係る流体解析方法に対応し得る。処理フローＳ６は学習済みモデルの生成方法の一例であるともいえる。

図８を参照しながら、学習部１１により実行される処理フローＳ６を説明する。学習部１１は、ステップＳ６１では一つの液体モデルを示す教師データに基づいて計算条件を取得し、ステップＳ６２ではその教師データに基づいて個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）を取得する。教師データの取得方法が限定されないのは処理フローＳ２と同様である。ステップＳ６１，Ｓ６２はそれぞれステップＳ１１，Ｓ１２に対応する。

ステップＳ６３では、学習部１１は個々の粒子について特徴量を求める。この例では、特徴量は以下の方針に基づいて計算される。

・物理量の平均場近似を考慮する。
・可能な変数をローカル空間にマッピングし、特徴量を計算する際に該変数を参照する。
・圧力場の計算方法はＭＰＳ陽解法に基づく。
・粒子間距離の計算では、平方根を使わない近似計算を用いる。
・重み関数として、距離の逆平方根を含まないモデルを用いる。

相互に左右し合う粒子の組合せを効率的に探すために、学習部１１はシミュレーション空間を、近傍粒子の距離の基準値ｒ_ｅを１辺の長さとする格子状の複数の部分空間に分割する。この部分空間をバケットともいう。

学習部１１は個々のバケットについて代表物理変数を定義してこの変数をバケットにマッピングする。このマッピングは平均場近似に基づく。或る一つのバケットをバケットＪと表すとすると、学習部１１は式（１２），（１３），（１４）で示される三つの代表物理変数ｒ_Ｊ ^ｋ，ｕ_Ｊ ^ｋ，Ｐ_Ｊ ^ｋを定義し、これらの変数をバケットＪにマッピングする。

ここで、ｎ_Ｊ ^ｋはバケットＪに含まれる粒子の数を示し、ｊは粒子番号を示す。

学習部１１は、相互に隣接するバケット内に属する二つの粒子ｉ，ｊ間の距離を、原点を共有するローカル座標系（ｘ_ｌ，ｙ_ｌ，ｚ_ｌ）で定義された点（ｘ_ｌｉ，ｙ_ｌｉ，ｚ_ｌｉ），（ｘ_ｌｊ，ｙ_ｌｊ，ｚ_ｌｊ）を用いて、式（１５）により近似する。

ここで、ｃ０～ｃ４は、相互に隣接するバケット同士の位置関係に応じた定数である。式（１５）ではｉ，ｊの成分が互いに独立なので、式（１６），（１７）によって各バケットへのマッピングが可能である。

学習部１１は、式（１８）で示される重み関数を用いる。式（１８）は、二つの粒子間の距離の逆平方根を含まないモデルを示し、このモデルを用いることで、計算を安定的に且つ高速に実行することが可能になる。

ここで、ｃ，Ｎはいずれも、特徴量に応じた定数である。

本実施形態では、学習部１１は、粘性力に対応する特徴量（粘性力特徴量）を求める。粘性力特徴量は、粒子に作用する力の特徴量の一例である。バケットＪからの寄与（∇^２ｕ_ｉ）_Ｊ ^ｋに対応する、粒子ｉの粘性力特徴量Φ^ｋ _{ｖｉｓｃ，ｉ，Ｊ}は、式（１２），（１３）の結果を用いて式（１９）で定義される。

学習部１１はｎ_Ｊ ^ｋ，ｕ_Ｊ ^ｋ，ｒ_Ｊ ^ｋをマッピング変数として参照することで、式（１９）を高速に計算することができる。

本実施形態では、学習部１１はさらに、圧力に対応する特徴量（圧力特徴量）を計算する。圧力特徴量も、粒子に作用する力の特徴量の一例である。バケットＪからの寄与（∇Ｐ_ｉ）_Ｊ ^ｋ＋１に対応する、粒子ｉの圧力特徴量Φ^ｋ _{ｐｒｅｓ，ｉ，Ｊ}は、式（２０）で定義される。

ここで、記号“＊”は、計算ステップｋ＋１の値が確定する前の仮変数であることを表す。

変数Ｐ_ｉ ^ｋ＋１，Ｐ_Ｊ ^ｋ＋１はＭＰＳ陽解法に基づいて得られる。

変数ｎ_ｉ ^＊，ｎ_Ｊ ^＊は、ＭＰＳ法の定義に従って重み関数を足し合わせたものに平均場近似を適用することで得られる。

ｎ_Ｊ ^ｋ，ｎ_Ｊ´ ^ｋはマッピング変数として参照される。ｒ_ｉ ^＊，ｒ_Ｊ ^＊は次式により表される。

Φ^ｋ _{ｖｉｓｃ，ｉ，Ｊ}，ｒ_Ｊ ^ｋ，ｕ_Ｊ ^ｋ，Φ^ｋ _{ｖｉｓｃ，Ｊ，Ｊ´}はマッピング変数として参照される。Φ^ｋ _{ｖｉｓｃ，Ｊ，Ｊ´}は、バケットＪに隣接するバケットＪ´を用いて式（２１）で表される。

ｎ_Ｊ´ ^ｋ，ｕ_Ｊ´ ^ｋ，ｕ_Ｊ ^ｋ，ｒ_Ｊ´ ^ｋ，ｒ_Ｊ ^ｋはマッピング変数として参照される。

ステップＳ６４では、学習部１１は複数の粒子のそれぞれの粘性力特徴量および圧力特徴量を示す入力ベクトルを用いて機械学習を実行する。

この機械学習に対応する理論は次の通りである。計算ステップｋにおいて粒子ｉに作用する加速度ａ_ｉ ^ｋは式（２２）により得られる。Ｒｅｇｒｅｓｓｉｏｎ（）は回帰式を表す。

速度変化量ｕ´_ｉは式（２３）で得られる。

粒子ｉの速度ｕ_ｉ ^ｋ＋１は式（２４）で示され、粒子ｉの位置ｒ_ｉ ^ｋ＋１は式（２５）で示される。

式（２２）で考慮されるバケットＪの範囲はシミュレーション空間の次元により異なり、これに応じて機械学習の説明変数の個数も異なる。シミュレーション空間が２次元であれば、一つの基準バケットは８個の隣接バケットで囲まれるので９個（＝３×３）のバケットが考慮される。個々のバケットについて粘性力特徴量および圧力特徴量が得られるので、機械学習の説明変数は１８個である。シミュレーション空間が３次元であれば、一つの基準バケットは２６個の隣接バケットで囲まれるので２７個（＝３×３×３）のバケットが考慮され、したがって、機械学習の説明変数は５４個である。

学習部１１は、複数の粒子のそれぞれの粘性力特徴量および圧力特徴量を示す入力ベクトルを用いて機械学習を実行する。具体的には、学習部１１は、複数の粒子の粘性力特徴量および圧力特徴量を示す入力ベクトルを機械学習モデルに入力し、第２時点における該複数の粒子の位置および速度を示す出力ベクトルを得る。粘性力特徴量および圧力特徴量は、粒子の第１時点における位置および速度に基づいて得られる。したがって、粘性力特徴量および圧力特徴量を示す入力ベクトルも、複数の粒子のそれぞれの第１時点における位置および速度に基づく。処理フローＳ２と同様に、機械学習モデルに入力される複数の粒子は、教師データで示される全粒子であってもよいし、ランダム選択などの手法により選択された一部の粒子のみであってもよい。

粒子の総数をＮ、時間ステップ毎に進む時刻を計算ステップｋ、粘性力特徴量をΦ^ｋ _ｖｉｓｃ、圧力特徴量をΦ^ｋ _ｐｒｅｓと表すとする。この場合、入力ベクトルは｛（Φ^ｋ _{ｖｉｓｃ，１}，Φ^ｋ _{ｐｒｅｓ，１}），（Φ^ｋ _{ｖｉｓｃ，２}，Φ^ｋ _{ｐｒｅｓ，２}），…，（Φ^ｋ _{ｖｉｓｃ，Ｎ}，Φ^ｋ _{ｐｒｅｓ，Ｎ}）｝と表され、出力ベクトルは｛（ｒ^ｋ＋１ _１，ｕ^ｋ＋１ _１），（ｒ^ｋ＋１ _２，ｕ^ｋ＋１ _２），…，（ｒ^ｋ＋１ _Ｎ，ｕ^ｋ＋１ _Ｎ）｝と表される。第２時点における粒子の位置での圧力を求めてもよく、この場合には、出力ベクトルは｛（ｒ^ｋ＋１ _１，ｕ^ｋ＋１ _１，Ｐ^ｋ＋１ _１），（ｒ^ｋ＋１ _２，ｕ^ｋ＋１ _２，Ｐ^ｋ＋１ _２），…，（ｒ^ｋ＋１ _Ｎ，ｕ^ｋ＋１ _Ｎ，Ｐ^ｋ＋１ _Ｎ）｝と表される。

この機械学習でも、教師データで示される第２時点での値が正解として用いられる。学習部１１は出力ベクトルと正解との差（すなわち、誤差）に基づいて機械学習モデル内のパラメータを更新する。

ステップＳ６３，Ｓ６４で示される一連の処理は、計算ステップｋ（第１時点）での値を示す入力から、計算ステップｋ＋１（第２時点）での値を示す出力を得る処理である。学習部１１はこの一連の処理を一つの計算ステップごとに実行してもよいし、すべての計算ステップについて一度に実行してもよい。すべての計算ステップについて一度に実行する場合には、入力ベクトルは｛（Φ^１ _{ｖｉｓｃ，１}，Φ^１ _{ｐｒｅｓ，１}，Φ^１ _{ｖｉｓｃ，２}，Φ^１ _{ｐｒｅｓ，２}，…，Φ^１ _{ｖｉｓｃ，Ｎ}，Φ^１ _{ｐｒｅｓ，Ｎ}）,（Φ^２ _{ｖｉｓｃ，１}，Φ^２ _{ｐｒｅｓ，１}，Φ^２ _{ｖｉｓｃ，２}，Φ^２ _{ｐｒｅｓ，２}，…，Φ^２ _{ｖｉｓｃ，Ｎ}，Φ^２ _{ｐｒｅｓ，Ｎ}）,…, （Φ^ｋ _{ｖｉｓｃ，１}，Φ^ｋ _{ｐｒｅｓ，１}，Φ^ｋ _{ｖｉｓｃ，２}，Φ^ｋ _{ｐｒｅｓ，２}，…，Φ^ｋ _{ｖｉｓｃ，Ｎ}，Φ^ｋ _{ｐｒｅｓ，Ｎ}）｝と表される。また、出力ベクトルは｛（ｒ^２ _１，ｕ^２ _１，ｒ^２ _２，ｕ^２ _２，…，ｒ^２ _Ｎ，ｕ^２ _Ｎ），（ｒ^３ _１，ｕ^３ _１，ｒ^３ _２，ｕ^３ _２，…，ｒ^３ _Ｎ，ｕ^３ _Ｎ），…，（ｒ^ｋ＋１ _１，ｕ^ｋ＋１ _１，ｒ^ｋ＋１ _２，ｕ^ｋ＋１ _２，…，ｒ^ｋ＋１ _Ｎ，ｕ^ｋ＋１ _Ｎ）｝と表される。圧力を含む出力ベクトルも同様に表現できる。

学習フェーズでは学習部１１は複数の液体モデルのそれぞれの教師データを用いて処理フローＳ６を実行し得る。任意の個数の液体モデルについて処理フローＳ６を実行した後に、学習部１１は処理された機械学習モデルを学習済みモデルとして取得する。この処理は、学習済みモデルが生成されたことを意味する。処理フローＳ２と同様に、学習部１１はその学習済みモデルを任意の手法で出力する。

図９を参照しながら、予測部１２により実行される処理フローＳ７を説明する。予測部１２は、ステップＳ７１では予測しようとする液体モデルの計算条件を取得し、ステップＳ７２では該液体モデルの個々の粒子の初期条件（ｒ^０，ｕ^０，Ｐ^０）を取得する。液体モデルの取得方法が限定されないのは処理フローＳ３と同様である。ステップＳ７１，Ｓ７２はそれぞれステップＳ１１，Ｓ１２に対応する。

ステップＳ７３では、予測部１２は個々の粒子について粘性力特徴量および圧力特徴量を求める。ステップＳ７３はステップＳ６３に対応する。

ステップＳ７４では、予測部１２は、複数の粒子の粘性力特徴量および圧力特徴量を示す入力ベクトルを学習済みモデルに入力することで出力ベクトルを生成する。具体的には、予測部１２は、全粒子の粘性力特徴量および圧力特徴量を示す入力ベクトルを学習済みモデルに入力し、第２時点における全粒子の位置および速度を示す出力ベクトルを得る。ステップＳ７４はステップＳ６４に対応する。粒子の特徴量は、該粒子の第１時点における位置および速度に基づいて得られる。したがって、複数の粒子のそれぞれの粘性力特徴量および圧力特徴量を示す入力ベクトルも、複数の粒子のそれぞれの第１時点における位置および速度に基づく。

第１時点の初期値は粒子の初期条件に対応する時点であり、その後、第１時点は時間ステップΔｔずつ進む。１ループ目では、入力ベクトルは、初期条件に対応する全粒子の位置および速度に基づいて算出された全粒子の粘性力特徴量および圧力特徴量を含んで構成される。２ループ目以降では、入力ベクトルは、前のループで得られた出力ベクトルで示される全粒子の位置および速度に基づいて算出された全粒子の粘性力特徴量および圧力特徴量を含んで構成される。処理フローＳ６に対応して、出力ベクトルは圧力を含んでもよいし含まなくてもよい。

ステップＳ７５では、予測部１２は計算を終了する否かを判定する。具体的には、予測部１２は予測における所与の最終時刻に到達したか否かを判定する。第２時刻が最終時刻でない場合には（ステップＳ７５においてＮＯ）、処理はステップＳ７６に進む。ステップＳ７６では、予測部１２は時間ステップΔｔだけ時間を進める。そして、予測部１２はその進めた時間を第１時点として、ステップＳ７３，Ｓ７４の処理を再び実行する。

計算を終了すると判定した場合、すなわち所与の最終時刻に到達した場合には、処理はステップＳ７７に進む。ステップＳ７７では、予測部１２は予測データを出力する。処理フローＳ３と同様に、予測データのデータ構造、出力方法、および表現方法はいずれも限定されない。

［効果］
以上説明したように、本開示の一側面に係る流体解析システムは少なくとも一つのプロセッサを備える。少なくとも一つのプロセッサは、液体を構成する複数の粒子のそれぞれの第１時点における位置および速度に基づく入力ベクトルを取得し、入力ベクトルを機械学習モデルに入力することで、第１時点から所与の時間ステップが経過した第２時点における、複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成する。

本開示の一側面に係る流体解析方法は、少なくとも一つのプロセッサを備える流体解析システムにより実行される。この流体解析方法は、液体を構成する複数の粒子のそれぞれの第１時点における位置および速度に基づく入力ベクトルを取得するステップと、入力ベクトルを機械学習モデルに入力することで、第１時点から所与の時間ステップが経過した第２時点における、複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成するステップとを含む。

本開示の一側面に係る流体解析プログラムは、液体を構成する複数の粒子のそれぞれの第１時点における位置および速度に基づく入力ベクトルを取得するステップと、入力ベクトルを機械学習モデルに入力することで、第１時点から所与の時間ステップが経過した第２時点における、複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成するステップとをコンピュータに実行させる。

このような側面においては、液体を構成する個々の粒子の位置および速度が機械学習を用いて計算される。粒子法を用いた予測の少なくとも一部に機械学習を導入することで、その予測を短時間で実行することが可能になる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、複数の粒子のそれぞれの第１時点における位置および速度から得られる、該複数の粒子のそれぞれに作用する力に基づいて、入力ベクトルを設定してもよい。粒子に作用する力に基づいて入力ベクトルを設定することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは力が粘性力および圧力を含んでもよい。粒子に作用する粘性力および圧力に基づいて入力ベクトルを設定することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、複数の粒子のそれぞれに作用する力の特徴量を示す入力ベクトルを設定し、特徴量が、粘性力に対応する粘性力特徴量と、圧力に対応する圧力特徴量とを含んでもよい。粒子に作用する力の特徴量（粘性力特徴量および圧力特徴量）を示す入力ベクトルを設定することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、流体が存在するシミュレーション空間を格子状に分割することで複数の部分空間を設定し、複数の粒子のそれぞれについて、該粒子に影響を及ぼす１以上の部分空間のそれぞれに対応する力を示す特徴量を算出し、複数の粒子のそれぞれについて算出された特徴量を示す入力ベクトルを設定してもよい。粒子に影響を与える可能性がある部分空間を設定して、その部分空間が粒子に及ぼす力を示す特徴量を考慮することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、複数の粒子のそれぞれに作用する力に基づいて、該複数の粒子のそれぞれの第２時点における仮位置を算出し、複数の粒子のそれぞれの仮位置に基づく特徴量を示す入力ベクトルを設定してもよい。第２時点における仮位置に基づく特徴量を示す入力ベクトルを設定することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは、複数の粒子のそれぞれの仮位置に基づく特徴量が、近傍粒子との距離に基づく特徴量と、粒子の仮位置での圧力を得るための連立一次方程式を示す行列またはベクトルに関連する特徴量と、仮位置での粒子密度に関連した特徴量と、仮位置での相対密度に関連した特徴量とを含んでもよい。これらの特徴量を示す入力ベクトルを設定することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、複数の粒子のそれぞれの第１時点における位置および速度を示す入力ベクトルを設定してもよい。粒子法を用いた予測の全体を機械学習で処理することで、その予測を短時間で実行することが可能になる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、複数の粒子のそれぞれの第１時点における圧力にさらに基づく入力ベクトルを取得してもよい。粒子の位置および速度だけでなく圧力も考慮することで、第２時点における粒子の位置および速度をより高い精度で予測することができる。

他の側面に係る流体解析システムでは、少なくとも一つのプロセッサが、複数の粒子のそれぞれの第２時点における圧力をさらに示す出力ベクトルを生成してもよい。粒子の位置および速度だけでなく圧力も出力することで、第２時点における粒子の状態をより詳細に予測することができる。

［変形例］
以上、本開示での実施形態に基づいて詳細に説明した。しかし、本開示は上記実施形態に限定されるものではない。本開示は、その要旨を逸脱しない範囲で様々な変形が可能である。

学習部および予測部の少なくとも一方は並列計算によって処理を実行してもよい。上述したように、少なくとも仮位置での圧力の計算を含む処理（すなわち、並列計算が困難な処理）が機械学習で実行されるので、処理フローＳ２～Ｓ７のそれぞれについて、その全体を並列に処理することが可能である。

上述したように、学習済みモデルはコンピュータシステム間で移植可能であり、流体解析システムは、学習フェーズおよび予測フェーズの双方を実行してもよいし、学習フェーズおよび予測フェーズのいずれか一方を実行しなくてもよい。したがって、流体解析システムは学習部１１および予測部１２のうちのいずれか一方に相当する機能要素を備えなくてもよい。

上記実施形態では粒子に作用する力の特徴量として粘性力特徴量および圧力特徴量を示したが、流体解析システムはこれら２種類の特徴量のうちのいずれか一つのみを用いてもよい。あるいは、流体解析システムはこれら２種類の特徴量の少なくとも一つと、別の種類の特徴とを用いてもよい。

少なくとも一つのプロセッサにより実行される流体解析方法の処理手順は上記実施形態での例に限定されない。例えば、上述したステップの一部が省略されてもよいし、別の順序で各ステップが実行されてもよい。また、上述したステップのうちの任意の２以上のステップが組み合わされてもよいし、ステップの一部が修正または削除されてもよい。あるいは、上記の各ステップに加えて他のステップが実行されてもよい。

流体解析システム内で二つの数値の大小関係を比較する際には、「以上」および「よりも大きい」という二つの基準のどちらを用いてもよく、「以下」および「未満」という二つの基準のうちのどちらを用いてもよい。このような基準の選択は、二つの数値の大小関係を比較する処理についての技術的意義を変更するものではない。

１０…流体解析システム、１１０…流体解析プログラム、１１…学習部、１２…予測部。

Claims (5)

1. 少なくとも一つのプロセッサを備え、
   前記少なくとも一つのプロセッサが、
   液体を構成する複数の粒子のそれぞれの第１時点における位置および速度から得られる、該複数の粒子のそれぞれに作用する力に基づいて、該複数の粒子のそれぞれの、前記第１時点から所与の時間ステップが経過した第２時点における仮位置を算出し、
   前記複数の粒子のそれぞれの前記仮位置に基づく特徴量を示す入力ベクトルを設定し、
   前記入力ベクトルを機械学習モデルに入力することで、前記第２時点における、前記複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成し、
   前記複数の粒子のそれぞれの前記仮位置に基づく特徴量が、近傍粒子との距離に基づく特徴量と、前記粒子の前記仮位置での圧力を得るための連立一次方程式を示す行列またはベクトルに関連する特徴量と、前記仮位置での粒子密度に関連した特徴量と、前記仮位置での相対密度に関連した特徴量とを含む、
   流体解析システム。
2. 前記少なくとも一つのプロセッサが、前記複数の粒子のそれぞれの前記第１時点における圧力にさらに基づく入力ベクトルを取得する、
   請求項１に記載の流体解析システム。
3. 前記少なくとも一つのプロセッサが、前記複数の粒子のそれぞれの前記第２時点における圧力をさらに示す出力ベクトルを生成する、
   請求項１または２に記載の流体解析システム。
4. 少なくとも一つのプロセッサを備える流体解析システムにより実行される流体解析方法であって、
   液体を構成する複数の粒子のそれぞれの第１時点における位置および速度から得られる、該複数の粒子のそれぞれに作用する力に基づいて、該複数の粒子のそれぞれの、前記第１時点から所与の時間ステップが経過した第２時点における仮位置を算出するステップと、
   前記複数の粒子のそれぞれの前記仮位置に基づく特徴量を示す入力ベクトルを設定するステップと、
   前記入力ベクトルを機械学習モデルに入力することで、前記第２時点における、前記複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成するステップと
   を含み、
   前記複数の粒子のそれぞれの前記仮位置に基づく特徴量が、近傍粒子との距離に基づく特徴量と、前記粒子の前記仮位置での圧力を得るための連立一次方程式を示す行列またはベクトルに関連する特徴量と、前記仮位置での粒子密度に関連した特徴量と、前記仮位置での相対密度に関連した特徴量とを含む、
   流体解析方法。
5. 液体を構成する複数の粒子のそれぞれの第１時点における位置および速度から得られる、該複数の粒子のそれぞれに作用する力に基づいて、該複数の粒子のそれぞれの、前記第１時点から所与の時間ステップが経過した第２時点における仮位置を算出するステップと、
   前記複数の粒子のそれぞれの前記仮位置に基づく特徴量を示す入力ベクトルを設定するステップと、
   前記入力ベクトルを機械学習モデルに入力することで、前記第２時点における、前記複数の粒子のそれぞれの位置および速度を示す出力ベクトルを生成するステップと
   をコンピュータに実行させ、
   前記複数の粒子のそれぞれの前記仮位置に基づく特徴量が、近傍粒子との距離に基づく特徴量と、前記粒子の前記仮位置での圧力を得るための連立一次方程式を示す行列またはベクトルに関連する特徴量と、前記仮位置での粒子密度に関連した特徴量と、前記仮位置での相対密度に関連した特徴量とを含む、
   流体解析プログラム。

Images