JP7008336B2

JP7008336B2 - 流体解析方法及び流体解析プログラム

Info

Publication number: JP7008336B2
Application number: JP2018559538A
Authority: JP
Inventors: 栄永田
Original assignee: E&IS CO., LTD.
Current assignee: E&IS CO., LTD.
Priority date: 2016-12-27
Filing date: 2017-12-26
Publication date: 2022-02-10
Anticipated expiration: 2037-12-26
Also published as: JPWO2018124125A1; WO2018124125A1

Description

本発明は、流体の流れや状態などの解析を行う流体解析方法及び流体解析プログラムに関する。

近年、新製品の開発現場における製品設計では、多様化されたニーズへの対応として多機種多仕様の開発が求められている。また、自動車などの製品設計においては、安全性能の向上に向けた新システム・構造の開発が必要である。しかしながら、現状の設計では、技術者の未経験領域で不具合事象が発生したり、関連部門とのコミュニケーション不足による最適仕様の確立が困難であるなどの問題がある。

これらの問題解決へのアプローチの一手段として、製品開発における仕様決定のためのシミュレーションおよび製品試験の手法を改善することが考えられる。すなわち、従来は見えなかったものを可視化することなどによって、製品の初期仕様の妥当性の判断をより容易かつ精度良く行うことができるようになり、革新的な構造の発想の着眼点を得ることができ、かつ開発現場での技術者のコミュニケーションの促進を図ることができる。このような革新的なシミュレーションの手法が必要とされている。

製品設計の際のシミュレーションの手法では、従来、流体等の連続体を解く数値計算手法として、格子を用いて微分方程式の近似解を求める有限差分法などの手法が知られている。近年では、これらの数値計算はＣＡＥ(Computer Aided Engineering)等の応用分野に適用され、流体と構造物とが相互作用する問題を解析するために用いられている。しかしながら、格子を用いる上記の手法では、自由表面等の界面が存在する問題や流体－構造連成問題等の移動境界が発生する場合には、取り扱いが複雑となる。また、格子を用いる上記の手法では、これらの場合についてプログラムの作成が困難となる場合がある。

すなわち、従来の既存ＣＡＥでは、ＣＡＤ設計→メッシュ化→分析、とメッシュ化のプロセスが必須であった。そのうえ、設計変更の反映があるごとに再度メッシュを作成する必要があり手間と時間がかかっていた。また、この方法では、製品における稼動部の取り扱いが複雑になるという問題がある。

これに対して、流体の挙動を解析する他の方法として粒子法が知られている（例えば、特許文献１を参照）。粒子法は、流体を複数の粒子の集まりとして表し、これら粒子間の相互作用の計算によって流体の挙動を解析する手法である。この粒子法は、有限体積法や有限要素法の計算において、手間のかかるメッシュ化のプロセスが不要であるという利点がある。

すなわち、従来のメッシュベース解析手法では、空間をメッシュで分割して解析していた。これは、各メッシュに流速・温度を計測するセンサーを配置するイメージである。これに対して、粒子法は、流体を粒子に変換して解析するため、空気・水など流体の粒子自体がセンサーを持つイメージである。このような粒子法によれば、各粒子の流速・温度変化がわかるので流体の流れを可視化することができる。

しかしながら、従来の粒子法を用いた流体の挙動解析(シミュレーション)では、複雑な演算処理が必要とされていたため、スーパーコンピュータなどの大型のコンピュータで多大な費用と時間をかけて解析を行わなければならなかった。そのため、実際の製品設計においては、粒子法を用いた流体の挙動解析を簡易かつ迅速に行うことができず、従来のメッシュ法による解析や実験室での実機の試験などを繰り返して製品の仕様を決めなければならなかった。そのため、製品設計に手間と時間がかかるうえ、必ずしも最適仕様を選択することができないおそれがあった。

特開２００８－１１１６７５号公報

本発明は上述の点に鑑みてなされたものであり、その目的は、より少ない演算負荷で簡便に流体の挙動を精度良く解析することが可能となることで、流体の解析に基づく製品の設計期間の大幅な短縮化を図ることができると共に、より最適な製品仕様を選択することができる流体解析方法及び流体解析プログラムを提供することにある。

上記課題を解決するため、本発明にかかる流体解析方法は、対象モデルに沿って流れる流体の流れ又は状態を解析するための流体解析方法であって、入力データとして、対象モデルの形状データと、前記対象モデルに沿って流れる流体を粒子に模した場合の当該粒子の特性に関する情報と、を入力するステップと、前記入力データを用いて粒子法により前記粒子の流れ又は状態の解析を行うステップと、前記流れ又は状態の解析の結果を用いて出力データを出力するステップと、前記出力データを用いた結果を画像処理することで、前記粒子の流れを可視化した画像情報を生成するステップと、を含むことを特徴とする。この場合、前記粒子の特性に関する情報は、前記粒子の座標、温度、速度、加速度の少なくともいずれかに関する情報であってよい。

そして、上記の粒子の流れ又は状態の解析を行うステップでは、コンピュータが備える制御手段（例えば、ＧＰＵやＣＰＵなど）の機能によってこの解析を行うに際して、当該解析における各演算に対して複数の演算処理を行うグリッド、及びこのグリッドの演算処理を管理するカーネルが個別に割り当てられられるようになっており、当該粒子の流れ又は状態の解析では、複数の粒子についての複数の演算を同時に行うことができるようになっているとよい。

また、この流体解析方法では、さらに、前記粒子が到達するターゲット領域の座標情報と、選択された前記粒子の識別情報及び座標情報と、前記出力データとを入力するステップと、前記ターゲット領域の座標情報を用いて前記粒子の識別情報のスキャニングを行うステップと、前記出力データ内の粒子の色分けを行うステップと、選択された粒子の識別情報と、前記出力データ内の色分けされた粒子のデータとを用いて、選択された粒子の前記ターゲット領域までの通過経路の追跡を行うステップと、を有していてもよい。

また、この流体解析方法では、さらに、前記出力ファイル内の粒子の色分けを行うステップは、任意に選択した平面状の焦点領域を前記粒子が通過したか否かに関する変数を作成するステップと、前記焦点領域を含む平面の方程式を求めるステップと、前記変数と前記方程式を用いて前記粒子が前記焦点領域を通過したか否かを判断するステップと、前記焦点領域を通過した後の前記粒子の色を他の色に変更するステップと、を含んでいてもよい。

また、本発明にかかる流体解析プログラムは、上記の流体解析方法をコンピュータに実行させるためのコンピュータプログラムである。

本発明にかかる流体解析方法及び流体解析プログラムによれば、粒子法アルゴリズムを用いた解析であるため、準備データ(入力データ)を用意する際にメッシュ化の処理が不要となる。また、従来の粒子法を用いた他の流体解析手法と比較して、演算(計算)処理の大幅な簡素化を図ることができる。

したがって、シミュレーションを行うためのハードウェアの構成の簡素化及び低価格化を実現することができる。すなわち、スーパーコンピュータなど大型のコンピュータが必要であった従来のシミュレーションと比較して、画像処理用演算プロセッサ(ＧＰＵ)を搭載可能なパーソナルコンピュータ(ＰＣ)など汎用のコンピュータで簡単にシミュレーションを実施することが可能となる。また、比較的に安価なパーソナルコンピュータ(ＰＣ)を複数台用いて並列計算を行うことが可能となるので、短時間で精度の高いシミュレーション結果を得ることができる。

また、本発明にかかる状態解析方法は、構造物の温度等の状態を解析するための状態解析方法であって、解析対象の構造物を有限個の静止粒子に置き換える処理を行う第１ステップと、前記解析対象の構造物に作用する流体を流体粒子に置き換える処理を行う第２ステップと、前記静止粒子と前記流体粒子それぞれの初期の状態量に関するデータを入力する第３ステップと、前記データを用いて粒子法により前記流体粒子の流れ又は状態の解析を行う第４ステップと、前記流体粒子から前記静止粒子に伝達される状態量を演算する第５ステップと、前記静止粒子に伝達された前記状態量に基づいて解析対象の構造物の状態量の変化を算出する第６ステップと、を含むことを特徴とする。

本発明にかかる状態解析方法では、解析対象の構造物を位置が変動しない静止粒子に置き換え、かつ解析対象の構造物に作用する流体を流体粒子に置き換えることによって、例えば、気相（流体粒子）から固相（静止粒子）への伝熱現象を模擬することが可能となる。この場合、データ入力ステップである第３ステップ及び演算プロセスである第５ステップなどを備えることにより、各流体粒子の時間幅毎の位置を精度良く推定することが可能となる。各流体粒子の位置を精度良く推定することが可能となることにより、各流体粒子と各静止粒子との間の粒子間距離を精度良く推定することが可能となる。粒子間距離を精度良く推定することが可能となることにより、所定時間幅内に有効半径内に存在する流体粒子から静止粒子へ伝達される伝熱量などの状態を粒子間距離に応じて精度良く重み付けすることが可能となる。これにより、気相（流体粒子）から固相（静止粒子）への熱量の総和（伝熱量）を精度良く推定することが可能となる。その結果、伝熱などの様々な状態量の変化の形態において解析対象物の状態及びその変化を精度良く推定することが可能となる。

また、上記状態解析方法では、解析対象物が固体である場合、その固体を有限個の静止粒子によって構成することが可能である。この場合、気相／固相間の伝熱形態と同様に固相／固相間の伝熱形態を精度良く模擬することが可能となる。これにより固相／固相間の伝熱形態において解析対象物の温度を精度良く推定することが可能となる。

また、上記状態解析方法では、解析対象の構造物の温度及びその変化を解析する場合には、当該構造物における流体が作用する表面及びその近傍の一部のみを静止粒子に置き換える処理を行うようにするとよい。これによれば、より少ない演算負荷で流体の作用による解析対象の構造物の状態変化を精度良く推定することが可能となるので、迅速かつ的確なシミュレーションを行うことが可能となる。

また、上記状態解析方法では、流体粒子は、静止粒子との間で粒子間距離に応じた相互作用を可能にする有効半径を有するようにすることが可能である。この場合、有効半径によって多数の流体粒子の中から静止粒子に作用することができる流体粒子を限定することが可能となる。静止粒子に作用する流体粒子が限定されることにより、流体粒子から静止粒子へ伝達される伝熱量を算出する際の演算負荷を軽減することが可能となる。

また、本発明にかかる状態解析プログラムは、上記状態解析方法をコンピュータに実行させるためのコンピュータプログラムであって、このコンピュータプログラムは複数の繰り返し演算プロセスを有し、各繰り返し演算プロセスはそのコンピュータ上で機能する専用のカーネル並びにその繰り返し演算プロセスを専門に行う専用の演算処理部を個別に有することを特徴とする。

上記状態解析プログラムでは、各繰り返し演算プロセスはそのコンピュータ上で機能する専用のカーネル並びにその繰り返し演算処理を専門に行う演算処理部を個別に有する。これにより、ホストコンピュータの演算処理装置（ＣＰＵ）の演算負荷を好適に軽減するのと同時に、汎用のコンピュータ上で複数の粒子についての複数の繰り返し演算プロセスを同時に実行することが可能となる。これにより、流体の解析に要する時間を大幅に短縮することが可能となり、最適な製品仕様を迅速に決定することができるようになる。

また、本発明にかかる粒子解析方法は、解析対象物に作用する流体を有限個の粒子で模擬し、該粒子の流れ又は状態を解析することで、流体の流れが解析対象物に対して及ぼす作用の影響を解析する粒子解析方法であって、粒子が作用する解析対象物に関する形状データ並びに粒子の初期値及び特性値に関するデータを入力するデータ入力プロセスと、データ入力プロセスで入力したデータを用いて粒子法により粒子の流れ又は状態の解析を行う粒子解析プロセスと、粒子のうちで解析対象物に作用を及ぼす影響粒子を選定する影響粒子選定プロセスと、解析対象物が影響粒子から受ける力である粒子力を算出する粒子力算出プロセスと、粒子力算出プロセスで算出された粒子力に基づいて、流体が解析対象物に対して及ぼす作用による影響を反映した指標を算出する影響指標算出プロセスと、を含むことを特徴とする。なお、ここでいう流体が解析対象物に対して及ぼす作用による影響を反映した指標は、後述する実施形態における「粒子抗力係数ＰＤ」及び「粒子揚力係数ＰＬ」が相当し、影響指標算出プロセスは、実施形態における「ＰＤ／ＰＬ計算プロセス（PD/PL Calculation Process）Ｐ３」が相当する。

本発明にかかる粒子解析方法によれば、上記のデータ入力プロセス及び粒子解析プロセスを備えることで、解析対象物に作用する流体を模擬した有限個の粒子それぞれの流れ状態（例えば、所定時間幅毎の位置、速度及び圧力など）を精度良く算出することができる。そのうえで、解析対象物が影響粒子から受ける力である粒子力を算出する粒子力算出プロセスと、粒子力算出プロセスで算出された粒子力に基づいて、流体が解析対象物に対して及ぼす作用による影響を反映した指標を算出する影響指標算出プロセスと、を備えることで、流体中における解析対象物の流体特性（例えば、空力特性）を粒子法による流体粒子の解析結果を用いて精度良く推定することが可能となる。これにより、設計対象物の流体特性（空力特性）に関する製品仕様などをより迅速かつ的確に決定することができるようになる。

また、格子を用いる従来の解析手法では、流体の流れに対する対象物の投影面の形状及び投影面積が変化する場合に対しては正確な解析を行うことができなかったところ、本発明の上記粒子解析方法によれば、静止している対象物や進行方向（ベクトル）が変化しない移動をする対象物など、流体に対する投影面の形状及び投影面積が変化しない対象物の周辺の流れに対する解析だけでなく、例えば、車両が旋回しながら走行する場合など流体に対する投影面の形状及び投影面積が刻々と変化する解析対象物の解析を精度良く行うことができるようになる。したがって、動的な解析対象物の周辺の流れ状態についても好適に模擬（解析）することが可能である。

また、本粒子解析方法の粒子選定プロセスでは、解析対象物の各表面から所定距離の範囲内に位置する粒子を影響粒子として選定することが可能である。この場合、所定距離を変えることにより、例えば、流体中における解析対象物の流体特性（例えば、空力特性、以下同じ。）に関する指標である上記流体が解析対象物に対して及ぼす作用による影響を反映した指標についての精度を向上させることが可能となる。

また、本粒子解析方法では、解析対象物の形状データは、解析対象物の表面を三角形で分割した形状モデルとすることが可能である。この場合、解析対象物の表面に接する接平面を容易に決定することができるため、解析対象物の表面と各粒子との距離を容易に算出することができる。これにより、解析対象物に作用を及ぼすことができる影響粒子を容易に選定することが可能となる。

さらに、影響粒子が解析対象物に作用する際の解析対象物の作用面積及び作用面の法線ベクトルを容易に算出することができるため、解析対象物が影響粒子から受ける粒子力、並びに上記の指標を容易に算出することが可能となる。

また、本粒子解析方法では、影響粒子の解析対象物に作用する作用面は、影響粒子の解析対象物への投影面とすることが可能である。影響粒子の解析対象物に作用する作用面が粒子の解析対象物への投影面である場合、解析対象物が影響粒子から受ける力が、解析対象物が流体から受ける流体力により近似するようになる。これにより、流体中おける解析対象物の流体特性を粒子法によって精度良く推定することが可能となる。

また、本発明にかかる粒子解析プログラムは、上記粒子解析方法をコンピュータに実行させるためのコンピュータプログラムであって、コンピュータプログラムは複数の繰り返し演算ステップを有し、各繰り返し演算ステップはそのコンピュータ上で機能する専用のカーネル並びにその繰り返し演算ステップを専門に行う専用の演算処理部を個別に有することを特徴とする。

上記粒子解析プログラムでは、各繰り返し演算ステップはそのコンピュータ上で機能する専用のカーネル並びにその繰り返し演算処理を専門に行う演算処理部を個別に有する。これにより、ホストコンピュータの演算処理装置（ＣＰＵ）の演算負荷を好適に軽減するのと同時に、汎用のコンピュータ上で複数の粒子についての複数の繰り返し演算ステップを同時に実行することが可能となる。これにより、流体の流れ又は状態解析に要する時間を大幅に短縮することが可能となり、最適な製品仕様を迅速に決定することができるようになる。

また、本発明にかかる粒子解析方法は、解析対象物に作用する流体を有限個の粒子で模擬し、該粒子の流れ又は状態を解析することで、前記流体の流れが前記解析対象物に対して及ぼす作用の影響を解析する粒子解析方法であって、前記粒子が作用する前記解析対象物に関する形状データ並びに前記粒子の初期値及び特性値に関するデータを入力するデータ入力処理（ＳＴ２－１）と、前記粒子のうちで前記解析対象物に作用を及ぼす影響粒子を選定する影響粒子選定処理（ＳＴ２－２）と、前記解析対象物が前記影響粒子から受ける力である粒子力を算出する粒子力算出処理（ＳＴ２－３）と、前記粒子力によって前記解析対象物の全体に作用する力を算出する全体力算出処理（ＳＴ２－４）と、前記解析対象物の各部に作用する力を算出する各部力算出処理（ＳＴ２－５）と、を含むことを特徴とする。

本発明にかかる粒子解析方法によれば、上記のデータ入力処理及び粒子解析処理を備えることで、解析対象物に作用する流体を模擬した有限個の粒子それぞれの流れ状態（例えば、所定時間幅毎の位置、速度及び圧力など）を精度良く算出することができる。そのうえで、解析対象物が影響粒子から受ける力である粒子力を算出する粒子力算出処理と、粒子力によって解析対象物の全体に作用する力を算出する全体力算出処理と、解析対象物の各部に作用する力を算出する各部力算出処理と、を有することで、流体中における解析対象物の流体特性（例えば、空力特性）を粒子法による流体粒子の解析結果を用いて精度良く推定することが可能となる。これにより、設計対象物の流体特性（空力特性）に関する製品仕様などをより迅速かつ的確に決定することができるようになる。

また、格子を用いる従来の解析手法では、流体の流れ中に有る解析対象物の全体に作用する力及び各部に作用する力の正確な解析を行うことができなかったところ、本発明の上記粒子解析方法によれば、影響粒子から解析対象物に作用する力を正確に算出できることで、解析対象物の全体及び各部に作用する力の解析を精度良く行うことができるようになる。したがって、流体の流れ中に有る解析対象物に作用する力をより正確に模擬（解析）することが可能である。

また、本粒子解析方法の粒子選定処理では、解析対象物の各表面から所定距離の範囲内に位置する粒子を影響粒子として選定することが可能である。この場合、所定距離を変えることにより、例えば、流体中における解析対象物の流体特性（例えば、空力特性、以下同じ。）に関する指標である上記流体が解析対象物に対して及ぼす作用による影響を反映した指標についての精度を向上させることが可能となる。

本発明にかかる流体解析方法及び流体解析プログラムによれば、より少ない演算負荷で簡便に流体の流れや状態変化を精度良く解析することが可能となることで、流体の流れや状態変化に基づく製品(対象モデル)の設計期間の大幅な短縮化を図ることができると共に、より最適な製品仕様を容易且つ迅速に選択することができるようになる。

本発明の一実施形態にかかる流体解析の処理を実行するためのコンピュータ装置の構成の一例を示すブロック図である。本発明の一実施形態にかかる流体の流れ状態解析（Simulation method）の全体の構成を示すブロック図である。図２のＳＴＬモデル作成プロセス（CAD model to STL）を示す説明図である。図２の粒子パラメータ設定プロセス（Initial fluid（particle） property and parameter setting）と、車両の特性パラメータ設定プロセス（Configuration setting for models）とを示す説明図である。図２の演算プロセスにおいてシミュレーションプログラムによって解かれるナビエ・スト－クス方程式（Navier stokes equation）を示す説明図である。ナビエ・ストークスの方程式（Navier stokes equation）を解くための演算プロセス（algorithm）を示すフロー図である。流体（粒子）の粒子可視化演算プロセス（Source tracking with color）を示すフロー図である。ＳＴＬ包絡平面の方程式を決定する算出方法を示す説明図である。粒子が，関連するＳＴＬ包絡平面に関し上流側（手前側）か、下流側（向こう側）のどちらの側に位置するかを決定する算出方法を示す説明図である。ＳＴＬ包絡平面上の粒子の通過点を求める算出方法を示す説明図である。粒子の通過点ＰがＳＴＬ部分平面の外部又は内部に位置する状態を示す図である。粒子の通過点（点Ｐ）がＳＴＬ三角平面上に位置しているか否かを判定する方法を示す説明図である。逆追跡演算プロセスの粒子ＩＤ特定プロセス（Particle ID scanning）を示す説明図である。逆追跡演算プロセスの逆追跡プロセス（Back tracking Process）を示す説明図である。上記の粒子の流れ解析（バックトラッキング）を用いた自動車のエンジンルーム内の最適設計の例を説明するための図である。本実施形態のシミュレーションプログラムとその演算処理を実行する画像演算処理装置を示す説明図である。図２の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）と、車両の特性パラメータ設定プロセス（Configuration setting for models）とを示す説明図である。本発明の静止状態解析方法（Static Particle analysis）の一実施形態にかかる静止粒子温度解析法の概要を示す説明図である。本実施形態にかかる静止粒子温度解析法の演算手順を示すフロー図である。本実施形態にかかる静止粒子温度解析法によるエンジンマウント部ブラケットの温度解析を示す説明図である。太陽からの輻射熱（Sun radiation）を考慮した気相／固相間伝熱における本静止粒子温度解析法による解析対象物の温度解析を示す説明図である。図２の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）と、車両の特性パラメータ設定プロセス（Configuration setting for models）とを示す説明図である。粒子解析方法（Particle analysis）における粒子影響解析の概要を示す説明図である。粒子影響解析において算出される粒子抗力係数を示す説明図である。粒子影響解析において算出される粒子揚力係数を示す説明図である。影響粒子が解析対象物に作用する際の解析対象物の作用面を示す説明図である。粒子抗力係数の検証結果を示すグラフである。粒子揚力係数の検証結果を示すグラフである。車両の旋回時における空力特性（空気抵抗力）の模擬（解析）を示す説明図である。本発明の粒子解析方法（Particle analysis）を用いた車両にかかる力の算出手順の全体の流れを示すフローチャートである。本発明の粒子解析方法（Particle analysis）を用いた車両にかかる力の算出手順で作成されるデータのファイルを示す図である。図３０のＳＴＬモデル作成プロセス（CAD model to STL）を示す説明図である。図３０の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）と、車両の特性パラメータ設定プロセス（Configuration setting for models）とを示す説明図である。車両にかかる力の中心の算出（CoFSolver）（ＳＴ２）の処理の流れを示すフローチャートである。車両にかかる力の算出（Force calculation）処理（ＳＴ２－３）の内容を示す図である。車両にかかる力の中心の算出（Center of Force calculation）処理（ＳＴ２－４）の内容を示す図である。車両にかかる力の釣り合い（Force Balance）を示す図である。車両の各車輪にかかる力（Force on each wheels）の算出処理（ＳＴ２－５）を説明するための図である。車両にかかる力の算出結果を示す図で、（ａ）は、車両の各車輪にかかる力を示す図、（ｂ）は、車体全体に作用する力を示す図である。車両の各車輪に作用する力の比較を示すグラフである。

以下、添付図面を参照して本発明の実施形態を詳細に説明する。
〔第１実施形態〕
図１は、本発明の一実施形態にかかる流体解析による処理を実行するためのコンピュータ装置の構成の一例を示すブロック図である。本実施形態の流体解析による処理は、予め用意されたシミュレーションプログラム１１２をコンピュータ装置１００で実行することによって実現できる。

コンピュータ装置１００は、各種演算処理を実行する中央演算処理装置ＣＰＵ(Central Processing Unit) (以下、「ＣＰＵ」と記す。)１０１と、主にグラフィック処理に関する演算を行うための画像演算処理装置 (Graphical Processing Unit) (以下、「ＧＰＵ」と記す。)１１５と、データ入力を受け付ける入力装置１０２と、表示装置１０３とを有する。また、コンピュータ装置１００は、記憶媒体からプログラム等を読み取る媒体読取装置１０４と、各種装置と接続するためのインタフェース装置１０５と、他の情報処理装置等と有線または無線により接続するための通信装置１０６とを有する。また、コンピュータ装置１００は、各種情報を一時記憶するＲＡＭ１０７と、ハードディスクドライブなどの記憶装置１０８とを有する。また、各装置１０１～１０８は、データが双方向に送信されるバス１０９に接続される。

記憶装置１０８には、後述する図２に示す各種解析処理を実行するシミュレーションプログラム１１２が記憶されている。また、記憶装置１０８には、シミュレーションプログラム１１２を実現するための各種パラメータデータ（流体の物性データ、流体が作用する物体の構造データ）が記憶される。入力装置１０２は、シミュレーションのための各種入力データの入力を受け付ける。表示装置１０３は、例えば液晶ディスプレイなどであって、シミュレーションの操作画面、結果画面等の各種画面を表示する。インタフェース装置１０５には、例えば印刷装置等が接続される。通信装置１０６は、例えば、他の情報処理装置と接続され、ユーザが遠隔地からリモート接続でコンピュータ装置１００を操作する情報等の各種情報をやりとりする。

なお、上記のシミュレーションプログラム１１２は、必ずしもコンピュータ装置１００内の記憶装置１０８に記憶されている必要はない。例えば、コンピュータ装置１００が読み取り可能な記憶媒体に記憶されたプログラムをコンピュータ装置１００が読み出して実行するようにしてもよい。ここでいうコンピュータ装置１００が読み取り可能な記憶媒体は、例えば、ＣＤ－ＲＯＭやＤＶＤディスク、ＵＳＢ(Universal Serial Bus)メモリ等の可搬型記録媒体、フラッシュメモリ等の半導体メモリ、外付けハードディスクドライブ等が対応する。

また、ＬＡＮ(Local Area Network)等の小規模(屋内)の通信ネットワーク１１３を介して接続されたデータベース１１０や、インターネットなどの大規模(屋外)の通信ネットワークＷを介して接続されたデータベース１１１にシミュレーションプログラムを記憶させておき、コンピュータ装置１００が通信装置１０６によるこれらデータベース１１０，１１１との通信によってシミュレーションプログラムを読み出して実行するようにしてもよい。

ＣＰＵ１０１及びＧＰＵ１１５は、ハードディスク装置１０８などに記憶されたシミュレーションプログラム１１２を読み出して、ＲＡＭ１０７に展開して実行することで各種の処理を行う。シミュレーションプログラム１１２は、コンピュータ装置１００を図２及び図４に示す本実施形態の流体の流れ解析に関する各種処理を実行する装置(ハードウェア)として機能させることができる。

すなわち、ＣＰＵ１０１又はＧＰＵ１１５は、シミュレーション対象(設計対象物)の形状データ、流体のデータ等を用いて、粒子法のアルゴリズムとして、例えばＭＰＳ (Moving particle semi-implicit) 法に基づいて、流体のシミュレーション計算を行う。ＭＰＳ法自体は公知であるためここではその詳細の説明は省略するが、ＭＰＳ法は、連続体の離散的な計算を、粒子間相互作用モデルを用いて行う。また、ＭＰＳ法では、粒子に対して所定の基準半径以内にある近接の粒子が、その粒子に対して相互作用を与えると仮定する。

図２は、本発明の一実施形態にかかる流体の流れ状態解析（Simulation method）の全体の構成を示すブロック図である。なお、ここでの解析対象は、車両（解析対象物）が一定速度で走行中にボディ正面のフロントグリル部（エアダクト）を通ってエンジンルーム内に流入する空気（粒子）の流れ状態（粒子の位置・速度・圧力）としている。

同図に示すように、本流れ解析は、後述する解析対象物の解析用形状モデル（以下、「ＳＴＬモデル（STL model）」という。）の準備ならびに粒子及び車両に関する各種パラメータの入力（設定）を行うデータ入力プロセスＳＴ１－１と、入力された各種パラメータに基づいてＧＰＵ１１５上でシミュレーションプログラム（solver）１１２を実行する演算プロセスＳＴ１－２と、シミュレーションプログラム１１２の実行結果（演算結果）を所定のコード及び所定のファイル形式で記憶装置１０８等の周辺装置に出力するデータ出力プロセスＳＴ１－３と、シミュレーションプログラム１１２の演算結果を基に流体（粒子）の流れ状態を可視化（着色化）する可視化プロセスＳＴ１－４と、を主に構成される。以下、各プロセスについて説明する。

データ入力プロセスＳＴ１－１は、ＳＴＬモデル作成プロセス（CAD model to STL）１２１と、流体（粒子）の初期状態量及び各種物性値を入力する粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）１２２と、解析対象物（車両）に関する各種特性パラメータを設定する解析対象物の特性パラメータ設定プロセス（Configuration setting for models）１２３と、から構成される。

図３は、図２のＳＴＬモデル作成プロセス（CAD model to STL）１２１を示す説明図である。図３（ａ）は、ＣＡＤ等の描画作成ソフトによって作成された車両に関する三次元ＣＡＤ形状モデル（３D CAD model）である。図３（ｂ）はＳＴＬ（STereoLithography）モデル（以下、「ＳＴＬモデル（STL model）」という。）である。このＳＴＬモデル（STL model）は、三次元ＣＡＤ形状モデル（３D CAD model）の表面（Surface）の一部又は全部を、図３（ｃ）に示される三角形平面（以下、「ＳＴＬ三角平面（STL triangle）という。」）で分割することによって作成される。

従って、解析対象物である車両のＳＴＬモデル（STL model）は、ＳＴＬ三角平面（STL triangle）を基に構成されている。そのため、詳細については後述するが、ＳＴＬモデルの表面に接する接平面の方程式を容易に決定することができ、その接平面によって粒子の追跡、可視化（着色化）が極めて容易となる。なお、１又は複数のＳＴＬ三角平面（STL triangle）を含むＳＴＬモデル（STL model）表面の一部分（Ａ部等）を、特にＳＴＬ部分平面（STL plate）と呼び、ＳＴＬ部分平面（STL plate）を含む接平面のことをＳＴＬ包絡平面（STL plane）と呼ぶことにする。

図４は、図２の粒子パラメータ設定プロセス（Initial fluid（particle） property and parameter setting）１２２と、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３とを示す説明図である。なお、ここでの粒子源（ソース）は、車両のボディ正面のフロントグリルを通ってエンジンルーム内に流入する粒子源（ソース１，２，３）としている。先ず粒子に関する初期値入力（Initial particle（fluid） setting）として、粒子間距離（Particle Distance）、半径影響係数（Radius Influence Coefficient）、各粒子の初期位置（Initial Particle Position）、初期速度（Initial Velocity）、初期圧力（Initial Pressure）、初期温度（Initial Temperature）及び初期粘度（Initial Viscosity）を入力する。なお、上記値は風洞試験結果（wind tunnel experiment）に基づいて決定される。これにより、各ソース１，２，３の初期条件が決定される。従って、後述する粒子の追跡演算（Source tracking）とは、これら各ソース１，２，３から流出する各粒子の流れ状態をシミュレーションプログラム１１２によって所定時間幅毎に推定・模擬することを意味している。

さらに、粒子に関するパラメータ（Parameter setting）として、流体の熱特性パラメータ（fluid Thermal Properties）、流体の移動パラメータ（Fluid Marching）、圧力計算パラメータ（Pressure Calculation）、熱特性パラメータに対する乗算係数（Fluid Thermal Properties LW）、空力特性パラメータ（Extra Aerodynamics）等を入力する。

さらに、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３として、粒子源データ（Bucket_dat）、エンジンルーム内部を構成する各部品の形状パラメータ、例えばラジエータ形状パラメータ（Radiator.stl）、ファン形状パラメータ（Mainfan.stl）、エンジン形状パラメータ（Engine.stl）、触媒形状パラメータ（CAT.stl）、ボディフレーム形状パラメータ（Framework_stl）、フロントグリル形状パラメータ（FR_Grille.stl）、等を設定する。

次に、粒子の初期値及びパラメータ、車両の特性パラメータを基に、粒子の流れ状態を追跡する図２の演算プロセスＳＴ１－２について説明する。

図５は、図２の演算プロセスＳＴ１－２においてシミュレーションプログラム１１２によって解かれるナビエ・スト－クス方程式（Navier stokes equation）を示す説明図である。図６はナビエ・ストークスの方程式（Navier stokes equation）を解くための演算プロセス（algorithm）を示すフロー図である。以下、各図について更に説明する。

図５（ａ）は、流体（粒子）の速度（ｕ）、密度（ρ）、圧力（Ｐ）及び粘度（ν）並びに重力（ｇ）の影響を考慮したナビエ・ストークスの方程式（以下、「一般式」という。）である。なお、この一般式の解である流体（粒子）の速度（ｕ）及び圧力（Ｐ）はベクトル量である。また、粒子の位置（ｒ）についてもベクトル量であり、粒子の速度（ｕ）を基に別途算出される。

図５（ｂ）は、一般式の右辺の第２項及び第３項を省略した、即ち流体（粒子）の圧力（Ｐ）のみを考慮したときのナビエ・ストークス方程式（以下、「第１特殊式」という。）である。また、図５（ｃ）は、一般式の右辺の第１項を省略した、即ち流体（粒子）の粘度（ν）及び重力（ｇ）の影響のみを考慮したときのナビエ・ストークス方程式（以下、「第２特殊式」という。）である。

本実施形態の演算プロセス（図６）では、一般式の解は、第１特殊式の解（ｕ’）と第２特殊式の解（ｕ_＊）とのベクトル和（＝ｕ_＊＋ｕ’）として求めている。以下、具体的な演算プロセスについて説明する。

図６に示されるように、シミュレーションプログラム１１２は、所定時間幅毎に、第ｋ回目の粒子の演算結果（ｒ_ｋ、ｕ_ｋ、Ｐ_ｋ）を基にして、最初に第２特殊式を解き、次に第１特殊式を解くことによって、第ｋ＋１回目の粒子の速度（ｕ_ｋ＋１）、圧力（Ｐ_ｋ＋１）及び位置（ｒ_ｋ＋１）を算出している。

従って、この演算プロセスは結果的に、所定時間幅毎に各粒子の流れ状態を追跡していることになる。従って、図６に示される演算プロセスは、粒子に関する「追跡演算プロセス（Source tracking Process）」である。以下、各ステップについて説明する。

ステップＳ１では、データ入力プロセスＳＴ１－１において入力された粒子のパラメータ及び解析モデルの形状パラメータを取り込む（Input Parameters）。

ステップＳ２では、粒子の初期位置（ｒ_０）、初期速度（ｕ_０）および初期圧力（Ｐ０）を取り込む（Input Initial r_0，u_0，P_0）。

ステップＳ３では、最新の粒子の位置（ｒ_ｋ）、速度（ｕ_ｋ）及び圧力（Ｐ_ｋ）を取り込む（r_k，u_k，P_k）。

ステップＳ４では、粒子の粘度（ν）と重力（ｇ）の影響を考慮した第２特殊式を解くことによって、粒子の速度（ｕ_＊）を算出する（Gravity and Viscosity factor calculation u_*）。

ステップＳ５では、ステップＳ４によって算出された速度（ｕ_＊）を基に、粒子の位置（ｒ_＊＝ｒ_ｋ＋Δｔ・ｕ_＊）を算出する（Particle Movement r_*=r_k+Δt・u_*）。

ステップＳ６では、粒子の圧力（Ｐ）のみの影響を考慮した第１特殊式を解くことによって、粒子の圧力（Ｐ_ｋ＋１）を算出する（Explicit Pressure Calculation）。

ステップＳ７では、ステップＳ６で第１特殊式の解として算出された粒子の圧力（Ｐ_ｋ＋１）を基に、粒子の速度（ｕ’）を算出する（Pressure gradient calculation）。

そして、ステップＳ４で第２特殊式の解として得られた粒子の速度（ｕ_＊）と、ステップＳ７で第１特殊式の解として得られた粒子の速度（ｕ’）との加算値（＝ｕ_＊＋ｕ’）を、第ｋ＋１番目の演算における粒子の速度（ｕ_ｋ＋１）とする（Velocity correction）。

同様に、ステップＳ５で第２特殊式の解（ｕ_＊）から算出された粒子の位置（ｒ_＊）と、ステップＳ６で第１特殊式の解（ｕ’）から算出された粒子の変位量（Δｔ・ｕ’）との加算値（＝ｒ_＊＋Δｔ・ｕ’）を、第ｋ＋１番目の演算における粒子の位置（ｒ_ｋ＋１）とする（Position correction）。

ステップＳ８では、データ出力の要求があるか否かを判定する（Is output ?）。データ出力の要求がある場合（ＹＥＳ）はステップＳ９に進み演算結果を出力する。データ出力の要求がない場合（ＮＯ）はステップＳ１０へ進む。

ステップＳ１０では、演算停止要求があるか否かを判定する（Is End ?）。演算停止要求がある場合（ＹＥＳ）は処理を終了する。演算停止要求がない場合（ＮＯ）はステップＳ１１へ進み、インデックスｋ＋１をｋに置換して次の時刻における粒子の位置、速度および圧力についての演算を実行する。

以上の通り、粒子追跡演算プロセス（Source tracking Process）では、所定時間幅毎に各粒子の位置、速度および圧力を算出し、その追跡演算結果を所定のコード及び所定のファイル形式で記憶装置１０８に保存している。

以降では、粒子追跡演算プロセス（Source tracking Process）で得られた各粒子の追跡演算結果（解析結果）を基に、流体（粒子）の流れを可視化（着色化）する粒子可視化演算プロセス（Source tracking Process with color）について説明する。

図７は、流体（粒子）の粒子可視化演算プロセス（Source tracking Process with color）を示すフロー図である。なお、説明を容易にするため、流体を構成する粒子の数は１個とする。また、各粒子の追跡演算結果は既に得られているものとする。

ステップＳ１では、粒子に関する追跡演算結果と、粒子に関連するＳＴＬ三角平面（STL triangles）を読み込む（Read the particles and STL（triangles））。

ステップＳ２では、ＳＴＬ三角平面（STL triangles）を包絡するＳＴＬ包絡平面の方程式（ａｘ＋ｂｙ＋ｃｚ＋ｄ＝０）を算出する（Find Equation of plane containing Plate:ax+bx+cz+d=0）。

図８は、ＳＴＬ包絡平面の方程式（ａｘ＋ｂｙ＋ｃｚ＋ｄ＝０）を決定する算出方法を示す説明図である。平面の方程式は、（１）平面に垂直な法線ベクトルＨ（＝（ａ，ｂ，ｃ））と、（２）平面上の任意の１点（ｘ１、ｙ１、ｚ１）を決定することにより、一意的に決定される。

図８（ａ）に示されるように、法線ベクトルＨは、ＳＴＬ三角平面（STL triangle）の頂点Ａ，Ｂ，Ｃの各座標を使用して求めることができる。すなわち、法線ベクトルＨは、ベクトルＡＢとベクトルＡＣの外積として求められる。

従って、図８（ｂ）に示されるように、ベクトルＨのｘ成分（＝ａ）は、ベクトルＡＢ＝頂点Ａ，Ｂ，Ｃの各座標を（Ａｘ、Ａｙ、Ａｚ），（Ｂｘ、Ｂｙ、Ｂｚ），（Ｃｘ、Ｃｙ、Ｃｚ）とするとき、
ａ＝（Ｂｙ－Ａｙ）*（Ｃｚ－Ａｚ）－（Ｃｙ－Ａｙ）*（Ｂｚ－Ａｚ）、として求められる。
同様に、ｙ成分（＝ｂ）およびｚ成分（＝ｃ）は、
ｂ＝（Ｂｚ－Ａｚ）*（Ｃｘ－Ａｘ）－（Ｃｚ－Ａｚ）*（Ｂｘ－Ａｘ）、
ｃ＝（Ｂｘ－Ａｘ）*（Ｃｙ－Ａｙ）－（Ｃｘ－Ａｘ）*（Ｂｙ－Ａｙ）、
としてそれぞれ求められる。

また、平面の方程式の定数項（＝ｄ）は、頂点Ａの座標（Ａｘ、Ａｙ、Ａｚ）を使用して、
ｄ＝－（ａ*Ａｘ＋ｂ*Ａｙ＋ｃ*Ａｚ）、として求められる。

再び図７に戻って、ステップＳ３では、粒子が、関連するＳＴＬ包絡平面に関し上流側（手前側）か、下流側（向こう側）のどちらの側に位置するかを決定する（Determine Position of particle wrt Plane）。

図９は、粒子が，関連するＳＴＬ包絡平面に関し上流側（手前側）か、下流側（向こう側）のどちらの側に位置するかを決定する算出方法を示す説明図である。
図９（ａ）は、点（ｘ０、ｙ０、ｚ０）と平面（ａｘ＋ｂｙ＋ｃｚ＋ｄ＝０）との距離を求める公式（｜ａ*ｘ０＋ｂ*ｙ０＋ｃ*ｚ０＋ｄ｜／√（ａ_２＋ｂ_２＋ｃ_２)から絶対値を外した式（以下、「相対位置判別式Ｄ」という。）に相当する。

相対位置判別式Ｄは、例えば、粒子（点）が平面の法線ベクトルＨに対向する側にある場合は、Ｄ＜０となり、その逆に粒子（点）が平面の法線ベクトルＨに対向しない側にある場合は、Ｄ＞０となる。

従って、図９（ｂ）に示されるように、相対位置判別式Ｄの符号を調べることにより、粒子が関連するＳＴＬ包絡平面に関し上流側（手前側）か、下流側（向こう側）のどちらの側に位置するかを判定することができる。同時に、粒子がＳＴＬ包絡平面を通過したか否かを判定することができる。

再び図７に戻って、ステップＳ４では、粒子がＳＴＬ包絡平面を通過したか否かを判定する（Particle Crossed the plane（position changes）?）。相対位置判別式Ｄの符号が負から正に変化するときに、粒子がＳＴＬ包絡平面を通過したと判定し、ステップＳ５を実行する。他方、相対位置判別式Ｄの符号が負のままで変化しない場合は、ステップＳ９へ進み、粒子は着色されない。そしてステップＳ１０（次の時刻のデータを読み込む）を経てステップＳ３を再度実行する。

ステップＳ５では、ＳＴＬ包絡平面上の粒子の通過点Ｐfootを求める（Find the projection（Pfoot）of Particle in Plane）。

図１０は、ＳＴＬ包絡平面上の粒子の通過点Ｐfootを求める算出方法を示す説明図である。

図１０（ａ）に示されるように、粒子の通過点Ｐfootは、直線（粒子の進行方向）と平面（ＳＴＬ包絡平面）との交点として求められる。今、粒子進行方向における粒子の座標をＰ＝（Ｐｘ、Ｐｙ、Ｐｚ）と、粒子の通過点Ｐfootの座標を（Ｐfoot.ｘ、Ｐfoot.ｙ、Ｐfoot.ｚ）と、粒子の方向ベクトルをＤ＝（Ｄｘ、Ｄｙ、Ｄｚ）とする。

図１０（ａ）に示されるように、粒子進行方向における粒子の座標Ｐは、媒介変数ｔを用いて、（１）Ｐ＝Ｐfoot－ｔ*Ｄ（ｔ：媒介変数）、と表すことができる。従って、Ｐfoot＝Ｐ＋ｔ*Ｄ＝（Ｐｘ、Ｐｙ、Ｐｚ）＋（ｔ*Ｄｘ、ｔ*Ｄｙ、ｔ*Ｄｚ）＝（Ｐｘ＋ｔ*Ｄｘ、Ｐｙ＋ｔ*Ｄｙ、Ｐｚ＋ｔ*Ｄｚ）＝（Ｐfoot.ｘ、Ｐfoot.ｙ、Ｐfoot.ｚ）、となる。
一方、ＰfootはＳＴＬ包絡平面上に位置するため、（２）ａ*Ｐfoot.ｘ＋ｂ*Ｐfoot.ｙ＋ｃ*Ｐfoot.ｚ＋ｄ＝０、を満足する。

従って、（１）を（２）に代入して整理すると、図１０（ｂ）１．に示されるように、Ｐfootに対する媒介変数ｔが求まる。そのｔを（１）に代入することにより、通過点Ｐfootのｘ座標Ｐfoot.ｘ、ｙ座標Ｐfoot.ｙ、ｚ座標Ｐfoot.ｚがそれぞれ求まる。

再び図７に戻って、ステップＳ６では、ステップＳ５で得られた粒子のＳＴＬ包絡平面上の通過点Ｐfootが、ＳＴＬ部分平面上の何れかのＳＴＬ三角平面上に位置するか否かを判定する。図１１（ａ）は、粒子の通過点ＰfootがＳＴＬ部分平面（STL plate）外に位置している場合を示している。図１１（ｂ）は、粒子の通過点ＰfootがＳＴＬ部分平面（STL plate）上の何れかのＳＴＬ三角平面上に位置している場合を示している。

図１２は、粒子の通過点Ｐfoot（点Ｐ）がＳＴＬ三角平面上に位置しているか否かを判定する方法を示す説明図である。図１２（ａ）は判定フローを示し、（ｂ）は（ａ）の各ステップの幾何学的意味を示す説明図であり、（ｃ）は（ａ）の各ステップで使用される演算要素の定義式である。

ステップＳ１では、ベクトルＶＷとベクトルＶＵの内積の符号（正、負）を判定する。（ｃ）に示されるように、ベクトルＶＷは、ベクトルＶ（ベクトルＡＣ）とベクトルＷ（ベクトルＡＰ）との外積を意味する。従って、ベクトルＶＷは、その大きさが三角形ＡＰＣの面積の２倍（２△ＡＰＣ）であり、その向きがベクトルＶからベクトルＷへ右ねじを回したときの進む方向である。ベクトルＶＵ、ベクトルＵＷ、ベクトルＵＶについても同様である。

ベクトルＶＵは定ベクトルで、その向きは紙面に垂直且つ紙面の反対側を指向する向きである。従って、ベクトルＶＷとベクトルＶＵの内積が負となるとは、ベクトルＶＷの方向が紙面に垂直且つ紙面の手前側を指向する向きであることと等価である。よって、ベクトルＶＷとベクトルＶＵの内積が負となるときの点Ｐの存在領域は、（ｂ）においてベクトルＶより上半分側の領域Ｒ１（縦線部分）となる。

ベクトルＶＷとベクトルＶＵの内積が負となる場合（ＹＥＳ）、点ＰはＳＴＬ三角平面（△ＡＢＣ）外に位置すると判定される。一方、ベクトルＶＷとベクトルＶＵの内積が負でない場合（ＮＯ）、ステップＳ２を実行する。

ステップＳ２では、ベクトルＵＷとベクトルＵＶの内積の符号（正、負）を判定する。ステップＳ１と同様に、ベクトルＵＶは定ベクトルで、その向きは紙面に垂直且つ紙面の手前側を指向する向きである。従って、ベクトルＵＷとベクトルＵＶの内積が負となるとは、ベクトルＵＷの方向が紙面に垂直且つ紙面の反対側を指向する向きであることと等価である。よって、ベクトルＵＷとベクトルＵＶの内積が負となるときの点Ｐの存在領域は、（ｂ）においてベクトルＵより左半分側の領域Ｒ２（横線部分）となる。

ベクトルＵＷとベクトルＵＶの内積が負となる場合（ＹＥＳ）、点ＰはＳＴＬ三角平面（△ＡＢＣ）外に位置すると判定される。一方、ベクトルＵＷとベクトルＵＶの内積が負でない場合（ＮＯ）、ステップＳ３を実行する。

ステップＳ３では、（ｃ）に示されるパラメータＲとパラメータＴとの加算値（＝Ｒ＋Ｔ）と、１との大小を判定する。

（ｃ）に示されるように、パラメータＲは、ベクトルＵＶの大きさ（＝２△ＡＢＣ）とベクトルＶＷの大きさ（＝２△ＡＰＣ）との比（＝△ＡＰＣ/△ＡＢＣ）である。同様に、パラメータＴは、ベクトルＵＶの大きさ（＝２△ＡＢＣ）とベクトルＵＷの大きさ（＝２△ＡＰＢ）との比（＝△ＡＰＢ/△ＡＢＣ）である。

従って、パラメータＲとパラメータＴとの加算値（＝Ｒ＋Ｔ）は、四角形ＡＢＰＣの面積（＝△ＡＰＣ＋△ＡＰＢ）と三角形ＡＢＣの面積（＝△ＡＢＣ）との比を意味している。従って、パラメータＲとパラメータＴとの加算値（＝Ｒ＋Ｔ）が１より大きいときの点Ｐの存在領域は、（ｂ）においてベクトルＶ及びベクトルＵによって挟まれた領域であって、線分ＢＣより右半分側の領域Ｒ３（斜線部分）となる。

従って、パラメータＲとパラメータＴとの加算値（＝Ｒ＋Ｔ）が１以下の場合（ＹＥＳ）、点ＰはＳＴＬ三角平面（△ＡＢＣ）内に位置すると判定される。他方、パラメータＲとパラメータＴとの加算値（＝Ｒ＋Ｔ）が１より大きい場合（ＮＯ）、点ＰはＳＴＬ三角平面（△ＡＢＣ）外に位置すると判定される。

再び、図７に戻って、ステップＳ６において、ＳＴＬ部分平面（STL plate）に含まれる全てのＳＴＬ三角平面（STL triangles）について、点ＰfootがＳＴＬ三角平面（STL triangles）内に位置するか否かを判定し、点Ｐfootが何れかのＳＴＬ三角形平面（STL triangles）内に位置する場合（ＹＥＳ）、ステップＳ７に進み、粒子は着色される。

一方、点Ｐfootが何れのＳＴＬ三角形平面（STL triangles）内にも位置しない場合（ＮＯ）、ステップＳ８に進み、粒子は着色されない。

なお、上記可視化演算結果は所定のコード及び所定のファイル形式で記憶装置１０８に保存される。

以上、図７～図１２を参照しながら説明した通り、粒子可視化演算プロセス（Source tracking Process with color）は、大きく４つの手順から構成される。一つめは、ＳＴＬ三角平面（STL triangles）を含むＳＴＬ包絡平面（STL plane）の方程式を決定する手順である。これは図７のステップＳ１～Ｓ２に相当する。

二つめは、粒子がＳＴＬ包絡平面を通過したか否かを判定する手順である。これは図７のステップＳ３～Ｓ４に相当する。

三つめは、粒子がＳＴＬ包絡平面を通過した場合、粒子のＳＴＬ包絡平面上の通過点Ｐfootを求める手順である。これは図７のステップＳ５に相当する。

四つめは、粒子の通過点ＰfootがＳＴＬ部分平面（STL plate）内に位置するか否かを判定し、粒子の通過点ＰfootがＳＴＬ部分平面（STL plate）内に位置する場合は粒子を着色し、それ以外の場合は着色しない手順である。これはステップＳ６～Ｓ８に相当する。

ところで、粒子追跡演算プロセス（Source tracking Process）では、流体（粒子）のソース（図４）を決定した後、すなわちソースを構成する各粒子の識別番号（以下、「粒子ＩＤ」という。）を決定した後に、粒子ＩＤに関する粒子の位置、速度および圧力を追跡して解析対象領域（Target region）における粒子の流れ状態を解析している。

以降では、粒子追跡演算プロセスとは逆に、既に取得された各粒子の流れ状態のデータ（位置、速度、圧力）から解析対象領域（Target region）を通過した粒子の粒子ＩＤを特定し、これら特定された粒子ＩＤに関する流れ状態のデータを出力データファイル（追跡演算結果）から抽出し、抽出した粒子ＩＤに関する流れ状態のデータを基に、解析対象領域（Target region）を通過した流体（粒子）のみによる流れ状態を新たに構築する演算プロセス（以下、「逆追跡演算プロセス（Reverse flow Process）」という。）について説明する。

この逆追跡演算プロセス（Reverse flow Process）は、解析対象領域（Target region）を通過した粒子のＩＤを特定する粒子ＩＤ特定プロセス（Particle ID scanning Process）と、各出力データファイルから特定された粒子ＩＤに関する流れ状態のデータを抽出して新たな出力データファイルを作成する逆追跡プロセス（Back tracking Process）と、から構成される。なお、逆追跡演算プロセス（Reverse flow Process）は、通常、粒子可視化演算プロセス（Source tracking Process with color）と一緒に実行される。粒子可視化演算プロセスが組み合わされた逆追跡演算プロセスのことを、特に「可視化逆追跡演算プロセス（Reverse flow Process with color coding）」と呼ぶことにする。

図１３は、逆追跡演算プロセスの粒子ＩＤ特定プロセス（Particle ID scanning Process）を示す説明図である。図１３（ａ）は粒子ＩＤ特定プロセスのフローを示し、（ｂ）は解析対象領域（Target region）の対角座標（diagonal coordinate）を示している。

ステップＳ１では、解析対象領域の対角座標を読み込む（Read diagonal coordinate of target region）と同時に、解析対象フレームを読み込む（Read Target Frame）。なお、ここで言う「対角座標（diagonal coordinate）」とは、解析対象領域（Target region）を包絡する直方体フレームの２つの対角点（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）を意味し、「解析対象フレーム（Target Frame）」とは、２つの対角点（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）によって形成された、直方体フレームを意味している。

ステップＳ２では、要求されたデータフレームから出力データファイル（追跡演算結果）を選定する（Load dat file of required frame）。ここでいう「データフレーム」とは、例えば粒子のソース（図４）に関する入力データファイル及び出力データファイルのリストを意味する。

ステップＳ３では、ステップＳ２で選定した出力データファイルから粒子の位置情報を抽出する（Get particle coordinate from file）。

ステップＳ４では、ステップＳ３で抽出した粒子が解析対象領域内部に位置しているか否かを判定する（If particle is inside target frame ?）。粒子が解析対象領域内部に位置している場合（ＹＥＳ）、ステップＳ５に進み、粒子のＩＤ（識別番号）をテキスト形式で保存する。

一方、粒子が解析対象領域（Target region）内部に位置していない場合（ＮＯ）、ステップＳ３を再度実行する。そして上記ステップＳ３からＳ５を全ての出力データファイルに対して実行する。

以上の通り、粒子ＩＤ特定プロセスでは、出力データファイル中の粒子の位置情報をスキャンすることにより、解析対象領域（Target region）を通過した粒子のＩＤ（識別番号）を特定する。以降では、特定された粒子ＩＤに関連する全ての出力データファイルをチェックして、粒子ＩＤに関する流れ状態のデータを抽出する逆追跡プロセス（Back tracking Process）について説明する。

図１４は、逆追跡演算プロセスの逆追跡プロセス（Back tracking Process）を示す説明図である。ステップＳ１では、ｃｓｖファイルをロードする（Load csv file）。なお、ここで言う「ｃｓｖファイル」とは、例えば粒子のソースに関する出力データファイルのリストを意味している。

ステップＳ２では、ｃｓｖファイルから１つの出力データファイルを読み込む（Read dat file from csv）。

ステップＳ３では、テキストファイルから粒子ＩＤを読み込む（Read particle id from txt file）。なお、ここで言う「テキストファイル」とは、上記粒子ＩＤ特定プロセスによって作成された、解析対象領域（Target region）を通過した粒子ＩＤのテキストリストである。

ステップＳ４では、ステップＳ３で読み込んだ粒子ＩＤが、ステップＳ２で読み込んだ出力データファイル内で有効であるか否かを判定する（If id is available in dat file ?）。なお、ここで言う「有効」とは、粒子ＩＤに関するデータがその出力データファイル内に存在することを意味している。粒子ＩＤがその出力データファイル内で有効である場合（ＹＥＳ）、ステップＳ５を実行する。一方、粒子ＩＤがデータファイル内で有効でない場合（ＮＯ）、ステップＳ３へ戻り、別の粒子ＩＤを読み込む。

ステップＳ５では、粒子ＩＤの情報を新たな出力データファイルに書き込む（Write particle info to the output dat file）。

ステップＳ６では、チェックすべき粒子ＩＤが残っているか否かを判定する（Any particle (id) remaining to check ?）。チェックすべき粒子ＩＤが残っている場合（ＹＥＳ）、ステップＳ３へ戻り、別の粒子ＩＤを読み込む。一方、チェックすべき粒子ＩＤが残っていない場合（ＮＯ）、ステップＳ７へ進み次の出力データファイルを読み込む（Read next dat file）。

ステップＳ８では、出力データファイルが目標フレーム（target frame）に達したか否かを判定する（If dat file reached to target frame ?）。出力データファイルが目標フレームに達した場合（ＹＥＳ）、処理を終了する。一方、データファイルが目標フレームに達していない場合（ＮＯ）、ステップＳ２に進み次の出力データファイルを読み込む。

以上の通り、逆追跡プロセスでは、各出力データファイルから、特定された粒子ＩＤに関する流れ状態のデータを抽出し、抽出した粒子ＩＤに関する流れ状態のデータを新たな出力データファイルに書き込んで、解析対象領域（Target region）を通過した粒子のみに関する流れ状態のデータを新たに構築している。

以上、本実施形態にかかる流体の流れ状態解析（Simulation method）によれば、解析対象物に関する解析用形状モデル（ＳＴＬモデル）が、ＳＴＬ三角平面（STL triangles）を基に構成されている。そのため、ＳＴＬ三角平面（STL triangles）を含む接平面であるＳＴＬ包絡平面（STL plane）の方程式が容易に決定される。ＳＴＬ包絡平面（STL plane）の方程式が容易に決定されることにより、粒子の追跡を各ＳＴＬ包絡平面（STL plane）を基準にして行うことができるようになる。その結果、粒子が各ＳＴＬ包絡平面（STL plane）に関し上流側または下流側のどちら側に位置しているのか、粒子がＳＴＬ包絡平面（STL plane）を通過したか否か、粒子がＳＴＬ包絡平面（STL plane）を通過したときの粒子の通過点Ｐfoot、並びに粒子がＳＴＬ三角平面（STL triangles）を通過したか否か、等が容易に判定することができるようになる。特に、粒子がＳＴＬ三角平面（STL triangles）を通過したか否かを容易に判定することができるため、粒子の可視化（着色化）が容易となる。

また、既に取得された解析対象領域（Target region）における粒子の流れ状態（位置、速度、圧力）から解析対象領域（Target region）を通過した粒子の粒子ＩＤを特定し、これら特定された粒子ＩＤに関する流れ状態のデータを出力データファイル（追跡演算結果）から抽出し、抽出した粒子ＩＤに関するデータを基に、解析対象領域（Target region）を通過した流体（粒子）のみによる流れ状態を新たに構築することができるようになる。従って、この新たに作成された出力データを用いることにより、例えば解析対象領域（Target region）を冷却するための熱設計等にフィードバックすることが可能となる。

すなわち、本実施形態にかかる流体の流れ状態解析（Simulation method）によれば、解析対象物の流れ状態を精度良く推定・模擬することができるようになる。

図１５は、上記の粒子の流れ解析(バックトラッキング)を用いた自動車のエンジンルーム内の最適設計の例を説明するための図である。同図（ａ）に示すように、エンジンルーム内の気体（空気）の解析計算を行う。この際、同図（ｂ）に示すように、各粒子の経路における温度・速度の変化を算出する。ここでは、ある粒子の経路における温度と速度の変化を解析したものを例示している。これにより、同図（ｃ）に示すように、流入源( Source1, Source2 )から流入してインテーク(ターゲット領域)に達した空気の経路の解析(経路の自動算出)が可能となる。

また、図１５（ｂ）の各粒子の速度・温度の経路の算出を用いて、図１５（ｄ）に示すように、部品形状の寄与度を算出することができる。ここでは、流入源 (Source1) から流入した（フロントグリルを出た）空気の流量の８０％がクーリングダクトに達し、クーリングダクトに達した空気はそこでＸ１［Ｗ］の熱交換が行われる。その後、クーリングダクトから出た空気の流量の６０％がインテーク（ターゲット領域）に達する。また、流入源 (Source2) から流入した（バンパーから出た）空気の流量の９０％がコンデンサに達する。コンデンサに達した空気はそこでＹ１［Ｗ］の熱交換が行われる。その後、コンデンサから出た空気の流量の９５％がラジエータに達する。ラジエータに達した空気はそこでＹ２［Ｗ］の熱交換が行われる。その後、ラジエータから出た空気の流量の６０％がトランスミッション(Ｔ／Ｍ)に達する。トランスミッションに達した空気はそこでＹ３［Ｗ］の熱交換が行われる。その後、トランスミッションから出た空気の流量の４０％がインテーク（ターゲット領域）に達する。このようなシミュレーション結果を得ることができる。したがって、このシミュレーション結果を用いてエンジンルーム１０内の部品の配置や形状の決定、あるいは温度耐性を考慮した材質選択などを行うことが可能となる。

このように、各粒子における部品との熱交換・速度変化を逐次記録することで、エンジンルーム１０内の対象箇所に対する部品形状の温度・速度変化の寄与度を算出ことに利用することができる。これにより、流れ・温度に対する部品形状及びレイアウトの寄与度を数値として算出することで、最適な設計案の選択を容易に行うことが可能となる。

図１６は、本実施形態のシミュレーションプログラムとその演算処理を実行する画像演算処理装置を示す説明図である。同図に示すように、シミュレーションプログラム(Solver)１１２の各演算ステップ（１）～（１２）は、ＣＰＵ１０１上で機能する専属のカーネル(Kernel)を個別に持っている。各カーネル(Kernel)は、ＧＰＵ１１５上の１つのグリッド(Grid)の演算処理を専門に管理するように構成されている。そして、各グリッドは、演算処理を行うスレッド(Thread)が１５個搭載されたブロック(Block)を６個搭載している。すなわち、シミュレーションプログラム１１２の各演算ステップ（１）～（１２）に対し、複数の演算処理を同時に行う専属のグリッド及びこのグリッドの演算処理を専門に管理するカーネルが個別に割り当てられている。そのため、複数の粒子についてシミュレーションプログラム１１２の複数の演算ステップを同時に実行することが可能となる。このように、シミュレーションプログラム１１２は、ＣＰＵ１０１と比較して圧倒的に多数の並列演算ユニットを有するＧＰＵ１１５を主に用いて、後述する粒子のトラッキングのための演算などを行うものである。このようなＧＰＵ１１５の並列演算ユニットでの並列計算によってパーソナルコンピュータなど汎用のコンピュータ装置で流体の粒子のトラッキングを短時間で精度良く行うことを可能としている。

すなわち、上記の粒子の流れ解析の処理（シミュレーション）は、ＧＰＵ１１５をＣＡＥ計算に応用することで行われている。ＧＰＵ１１５は、ＣＰＵ１０１と比較して圧倒的に多数の並列演算ユニットを備えている。これにより、ＧＰＵ１１５を搭載した汎用のパーソナルコンピュータ装置（ＰＣ）一台で上記のシミュレーションを実行可能となる。

したがって、上記のシミュレーションを行うためのハードウェアの構成の簡素化及び低価格化を実現することができる。すなわち、画像処理用演算プロセッサ（ＧＰＵ）を搭載可能なパーソナルコンピュータ（ＰＣ）で簡単にシミュレーションを実施することが可能となる。これにより、スーパーコンピュータなど大型のコンピュータが必要であった従来のシミュレーションと比較して、比較的に安価なパーソナルコンピュータ（ＰＣ）を複数台用いて並列計算を行うことが可能であるため、複数の形状に対するシミュレーションを同時並行で行うことができ、製品仕様の決定に要する時間を大幅に短縮することが可能となる。

以上、本発明の第１実施形態を説明したが、本発明は、上記実施形態に限定されるものではなく、特許請求の範囲、及び明細書と図面に記載された技術的思想の範囲内において種々の変形が可能である。例えば、上記実施形態では、本発明の流体解析方法を用いて自動車のエンジンルーム内の気体の流れを解析する例を示したが、本発明にかかる流体解析方法では、他にも様々な状況での流体の解析を行うことが可能である。例えば、パーソナルコンピュータの筐体内の気流の流れ、オフィス(室内)の空調設備による気流の流れの影響、自動車のトランスミッションケース内のオイルの流れ、自動車のボディの空気抵抗、などを解析することができる。さらには、空気などの流体だけでなく、部品など固体を粒子化し、その粒子の温度変化を解析することで、部品の温度変化などの解析を行うことも可能である。

〔第２実施形態〕
次に、本発明の第２実施形態について説明する。なお、第２実施形態の説明及び対応する図面においては、第１実施形態と同一又は相当する構成部分には同一の符号を付し、以下ではその部分の詳細な説明は省略する。また、以下で説明する事項以外の事項、及び図示する以外の事項については、第１実施形態と同じである。この点は、下記の他の実施形態についても同様である。
図１７は、図２の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）１２２と、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３とを示す説明図である。なお、ここでの粒子源（ソース）は、車両のデフロスタ吹出口を通ってフロントガラス内壁に吹き付けられる空気（温風）としている。先ず粒子に関する初期値入力（Initial particle（fluid） setting）として、粒子間距離（Particle Distance）、半径影響係数（Radius Influence Coefficient）、各粒子の初期位置（Initial Particle Position）、初期速度（Initial Velocity）、初期圧力（Initial Pressure）、初期温度（Initial Temperature）及び初期粘度（Initial Viscosity）を入力する。なお、上記値は風洞試験結果（wind tunnel experiment）に基づいて決定される。これにより、デフロスタ吹出口から流出する粒子（ソース）の初期条件が決定される。従って、後述する粒子の追跡演算（Source tracking）とは、デフロスタ吹出口から流出する各粒子の流れ状態をシミュレーションプログラム１１２によって所定時間幅毎に推定・模擬することを意味している。

さらに、粒子に関するパラメータ（Parameter setting）として、流体の熱特性パラメータ（fluid Thermal Properties）、流体の移動パラメータ（Fluid Marching）、圧力計算パラメータ（Pressure Calculation）、熱特性パラメータに対する乗算係数（Fluid Thermal Properties LW）、空力特性パラメータ（Extra Aerodynamics）等を入力する。なお、静止粒子についてのパラメータ入力は、図１９を参照しながら後述する。

さらに、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３として、例えば粒子源データ（Bucket_dat）、エンジンルーム内部を構成する各部品の形状パラメータ、例えばラジエータ形状パラメータ（Radiator.stl）、ファン形状パラメータ（Mainfan.stl）、エンジン形状パラメータ（Engine.stl）、触媒形状パラメータ（CAT.stl）、ボディフレーム形状パラメータ（Framework_stl）、フロントグリル形状パラメータ（FR_Grille.stl）、等を設定する。

次に、入力された粒子の初期値及びパラメータ並びに車両の特性パラメータを基に粒子の流れ状態を追跡する図２の演算プロセスＳＴ１－２について説明する。

以降では、粒子追跡演算プロセス（Source tracking process）によって得られた流体（粒子）の位置（ｒ^ｋ＋１）に基づき、流体（粒子）から車両の固体部位へ伝達される熱量の総和（伝熱量）を算出し、算出した伝熱量に基づいて、流体（粒子）が作用した車両の固体部位の時系列の温度分布を推定する静止粒子温度解析法について説明する。

図１８は、本発明の静止状態解析方法（Static Particle analysis）の一実施形態にかかる静止粒子温度解析法の概要を示す説明図である。図１８（ａ）に示されるように、ここでの温度解析対象は、車両のフロントガラス（Windshield）１３０の温度分布としている。そして、フロントガラス１３０に作用する粒子源は、デフロスタ吹出口からフロントガラス１３０へ吹き出される流体粒子１３２としている。

図１８（ｂ）に示されるように、温度解析対象物であるフロントガラス１３０のＳＴＬモデルを有限個の静止粒子（Static particles）１３１によって構成する。静止粒子（Static particles）１３１は、フロントガラス１３０のＳＴＬ三角平面（STL triangle）又はＳＴＬ部分平面（STL plate）を球体（Sphere）の粒子に変換したものである。なお、このＳＴＬモデルから静止粒子１３１への変換は、図２のＳＴＬモデル作成プロセス（CAD model to STL）１２１において行われる。

図１８（ｃ）に示されるように、各静止粒子１３１の温度Ｔ_ｉ、及び各流体粒子Ｔ_ｊは、有効半径（Effective radius）Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１へ伝達される伝熱量Ｈ_iを算出し、算出した伝熱量Ｈ_iを基に各静止粒子１３１及び各流体粒子１３２についての各熱方程式を解くことにより算出される。なお、伝熱量Ｈ_iを算出する際に必要となる流体粒子１３２の位置（ｒ_j）は、図６の粒子追跡演算プロセス（Source tracking process）から得られ、静止粒子１３１の位置（ｒ_i）は三次元ＣＡＤ形状モデル又はＳＴＬモデルから得られる。また、この熱方程式の演算は、図２の演算プロセスＳＴ１－２において行われる。

また、熱方程式を解く際に必要となる有効半径（Effective radius）Ｒｅ等の各種パラメータの入力は、図２の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）１２２において行われる。各種パラメータについては、図１９を参照しながら後述する。また、添字ｉは静止粒子１３１の帰属を意味し、添字ｊは流体粒子１３２の帰属を意味している。以降において同じ。

図１８（ｄ）に示されるように、静止粒子１３１についての熱方程式を解くことによって得られた各静止粒子１３１の温度｛Ｔ_i｝に基づいて、フロントガラス１３０の温度分布を着色化する。この着色化は図２の可視化プロセスＳＴ１－４において行われる。

図１９は、本実施形態にかかる静止粒子温度解析法の演算手順を示すフロー図である。なお、この静止粒子温度解析法では、有効半径Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１へ伝達される伝熱量Ｈ_iを用いて各静止粒子１３１の温度Ｔ_i，及び各流体粒子１３２の温度Ｔ_jを算出している。伝熱量Ｈ_iは粒子間距離ｒ_ijによって重み付けされて算出されるため、予め図６の粒子追跡演算プロセス（Source tracking process）を実行し、各流体粒子１３２の時系列位置（ｒ_j）を取得しておく必要がある。以下、各ステップについて説明する。

ステップＳ１では、静止粒子１３１及び流体粒子１３２についての熱パラメータを入力する（Input Thermal Parameters）。

熱パラメータとしては、流体粒子１３２（空気）の熱伝導率（fluid(air) thermal conductivity）λ_ｆ、静止粒子１３１（フロントガラス１３０）の熱伝導率（static wall thermal conductivity）λ_ｓ、静止粒子１３１（フロントガラス１３０）の熱抵抗（Static wall particle heat resistance）Ｒ、静止粒子１３１の粒子数密度（particle number density）ｎ_ｉ、空間次元数（number of space dimension）ｄ、時間幅（time step）Δｔ、有効半径（Effective radius）Ｒｅ、粒子間距離の重み関数（Weighting function）ω(ｒ_ij)、等を入力する。なお、この粒子間距離の重み関数ω(ｒ_ij)は、粒子間距離ｒ_ijに応じた重み係数を出力する関数であり、粒子間距離ｒ_ijが大きくなればなるほど重み係数はより小さくなり、粒子間距離ｒ_ijが有効半径Ｒｅ以上となるときに重み係数はゼロとなる。

ステップＳ２では、各静止粒子１３１の初期温度Ｔ^０ _i、ならびに各流体粒子１３２の初期温度Ｔ^０ _jをそれぞれ入力する（Assign Initial temperature）。

ステップＳ３では、ステップＳ１で入力された各熱パラメータ及びステップＳ２で入力された初期温度Ｔ^０ _i，Ｔ^０ _jを読み込んで、各静止粒子１３１及び各流体粒子１３２についての各熱方程式を解く（Calculate Heat equation）。熱方程式を解くことによって、第ｋ＋１時刻における各静止粒子温度Ｔ^ｋ＋１ _i、及び各流体粒子温度Ｔ^ｋ＋１ _jが算出される。熱方程式は公知であるため、ここではその説明は省略する。

熱方程式を解く際に必要となる、第ｋ＋１時刻までに有効半径Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１（添字ｉの静止粒子１３１）へ伝達された熱量の総和（伝熱量）Ｈ^ｋ＋１ _static-fluidは、下記式１で与えられる。

式１の右辺第１項は、第ｋ時刻までに有効半径Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１（添字ｉの静止粒子１３１）へ伝達された熱量の総和（伝熱量）Ｈ^ｋ _iを示し、同第２項は、時間幅Δｔ内（第ｋ時刻から第ｋ＋１時刻の間）に有効半径Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１へ伝達された熱量の総和（伝熱量）を示している。

式１の右辺第２項は、粒子間距離の重み関数ω(ｒ_ij)、第ｋ時刻における流体粒子温度Ｔ^ｋ _jと静止粒子温度Ｔ^ｋ _iとの温度差に比例する一方、粒子間距離ｒ_ijの２乗に反比例すると共に、静止粒子１３１の粒子数密度(particle number density)ｎ_ｉに反比例することが分かる。なお、１／λ_sf＋Ｒ／ｒ_ijは、いわゆる合成熱抵抗であり、値が大きいほど流体粒子１３２から静止粒子１３１へ伝達される熱量は小さくなる。

このように、第ｋ＋１時刻までに有効半径Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１へ伝達された伝熱量Ｈ^ｋ＋１ _static-fluidは、粒子間距離ｒ_ijに基づいて緻密に重み付けされて算出されることが分かる。

従って、ステップＳ３では、所定時間幅Δｔ毎に、第ｋ時刻における演算結果（Ｈ^ｋ _i、Ｔ^ｋ _i、Ｔ^ｋ _j）に基づいて、第ｋ＋１時刻までに有効半径Ｒｅ内に存在する流体粒子１３２から添字ｉの静止粒子１３１へ伝達された熱量の総和（伝熱量）Ｈ^ｋ＋１ _i、第ｋ＋１時刻における各静止粒子温度Ｔ^ｋ＋１ _i、及び各流体粒子温度Ｔ^ｋ＋１ _jをそれぞれ算出している。

ステップＳ４では、ステップＳ３で算出された各静止粒子温度Ｔ^ｋ＋１ _i、及び各流体粒子温度Ｔ^ｋ＋１ _jを更新する（Update static particle temperature）。

以上の通り、本実施形態にかかる静止粒子温度解析法によれば、各静止粒子温度Ｔ^ｋ＋１ _i、及び各流体粒子温度Ｔ^ｋ＋１ _jを算出する際に必要となる有効半径Ｒｅ内に存在する流体粒子１３２から各静止粒子１３１へ伝達された伝熱量が、粒子間距離ｒ_ijに基づいて緻密に重み付けされて算出されるため、各静止粒子温度Ｔ^ｋ＋１ _i、及び各流体粒子温度Ｔ^ｋ＋１ _jを精度良く推定することが可能となる。

このように、本実施形態にかかる静止粒子温度解析法によれば、流体の温度を精度良く推定・模擬することができると共に、流体が作用する物体の温度分布を精度良く推定・模擬することができる。

なお、上記実施例は、流体（流体粒子１３２）とフロントガラス１３０（静止粒子１３１）との間の気相／固相間伝熱におけるフロントガラス１３０の温度分布を本静止粒子温度解析法によって推定するものであった。以降では、固相／固相間伝熱における解析対象物の温度分布を本静止粒子温度解析法によって推定する実施例について説明する。

図２０は、本実施形態にかかる静止粒子温度解析法によるエンジンマウント部ブラケットの温度解析を示す説明図である。ここでの温度解析対象は、エンジン１５０と左右エンジンマウント部１５２Ｒ、１５２Ｌとをそれぞれ連結する左右ブラケット１５３Ｒ、１５３Ｌの各温度分布としている。

図２０（ａ）に示されるように、エンジン１５０はエンジン本体１５０ａと変速機１５０ｂとから構成される。エンジン１５０は、ボディ１５１に対し左右ブラケット１５３Ｒ、１５３Ｌを介して左右エンジンマウント部１５２Ｌ、１５２Ｒによって２点で支持されている。従って、エンジン１５０は、熱を供給する高温固体となり、左右エンジンマウント部１５２Ｌ、１５２Ｒ及びボディ１５１は熱を受ける低温固体となる。また以降において、左右エンジンマウント部１５２Ｌ、１５２Ｒ及び左右ブラケット１５３Ｌ、１５３Ｒは、特に区別する必要がある場合を除き左右の区別はしないこととする。

図２０（ｂ）に示されるように、ブラケット１５３を有限個の静止粒子１５４によって構成し、静止粒子ブラケット１５３’を形成する。エンジン１５０からの熱流は、有限個の静止粒子１５４を介してエンジンマウント部１５２及びボディ１５１へ伝導する。従って、エンジン１５０から静止粒子１５４に伝導される伝熱量Ｈを算出し、算出した伝熱量Ｈを公知の熱方程式に代入して解くことによって、静止粒子ブラケット１５３’（静止粒子１５４）の温度分布を推定することができる。

図２０（ｃ）は温度解析結果を示し、エンジン本体１５０ａ近傍に位置する左ブラケット１５３Ｌは高温となり、エンジン本体１５０ａから離れた右ブラケット１５３Ｒは低温となることを示している。この結果から、解析された温度分布と実際の温度分布とが近似していることが分かる。このように、本実施形態にかかる静止粒子温度解析法によれば、固相／固相間伝熱における解析対象物の温度を精度良く推定・模擬することができる。

続いて、太陽からの輻射熱（Sun radiation）を考慮した気相／固相間伝熱における本静止粒子温度解析法による解析対象物の温度解析について説明する。

図２１は、太陽からの輻射熱（Sun radiation）を考慮した気相／固相間伝熱における本静止粒子温度解析法による解析対象物の温度解析を示す説明図である。図２１（ａ）に示されるように、固相は車室（Cabin）のドアガラス（Door glass）１６０であり、気相（流体粒子１６２）は車室（Cabin）の空調機（A/C）１６３から吹き出される気流（Cabin air flow）である。従って、解析対象は、その気流（Cabin air flow）の温度としている。

図２１（ｂ）に示されるように、解析対象物であるドアガラス（Door glass）１６０を有限個の静止粒子（static particle）１６１によって構成する。

そして、静止粒子（static particle）１６１と流体粒子（fluid particle）１６２との間の伝熱量を算出し、算出した伝熱量を公知の熱方程式に代入して解くことによって各流体粒子（fluid particle）１６２の温度Ｔ_jを推定することができる。

なお、図１８に示されるフロントガラス１３０（静止粒子１３１）は、有効半径Ｒｅ内に存在する流体粒子１３２から静止粒子１３１へ伝達される伝熱量Ｈ^ｋ＋１ _static-fluidを主に受熱することとしていた。対する本実施例のドアガラス１６０（静止粒子１６１）は、有効半径内に存在する流体粒子１６２から静止粒子１６１に伝達される伝熱量以外に、太陽からの輻射熱（Sun radiation）を受熱することとし、気流（流体粒子１６２）も太陽光からの輻射熱（Sun radiation）を受熱することしている。

すなわち、太陽からの輻射熱（Sun radiation）の一部はドアガラス１６０（静止粒子１６１）によって吸収され、その残りはドアガラス１６０を透過して流体粒子１６２によって吸収される。従って、空調機（A/C）１６３から吹き出される気流（流体粒子１６２）は、空調機（A/C）１６３から受熱すると共に太陽からの輻射熱（Sun radiation）を受熱して熱エネルギーを増大させる。静止粒子１６１も太陽からの輻射熱（Sun radiation）を受熱して熱エネルギーを増大させる。従って、静止粒子１６１の温度Ｔ_iが、流体粒子１６２の温度Ｔ_jよりも高い場合はドアガラス１６０（静止粒子１６１）から気流（流体粒子１６２）へ熱エネルギーが伝達される。その逆に、ドアガラス１６０（静止粒子１６１）の温度Ｔ_iが、気流（流体粒子１６２）の温度Ｔ_jよりも低い場合は気流（流体粒子１６２）からドアガラス１６０（静止粒子１６１）へ熱エネルギーが伝達される。このように、本実施形態にかかる静止粒子温度解析法によれば、太陽からの輻射熱（Sun radiation）を考慮した気相／固相間伝熱において解析対象の温度を精度良く推定・模擬することができる。

以上の通り、本発明の状態解析方法によれば、解析対象物を構成する粒子は位置が変動する流体粒子１３２,１６２と、位置が変動しない静止粒子１３１,１５４,１６１とによって構成されることによって、対流を伴う伝熱現象、例えば気相（流体粒子）から固相（静止粒子）への伝熱現象を模擬することが可能となる。この場合、上記データ入力ステップＳＴ１－１及び上記演算プロセスＳＴ１－２を備えることにより、各流体粒子１３２,１６２の時間幅Δｔ毎の位置ｒ_ｊを精度良く推定することが可能となる。各流体粒子１３２,１６２の位置ｒ_ｊを精度良く推定することが可能となることにより、各流体粒子１３２,１６２と各静止粒子１３１,１６１との間の粒子間距離ｒ_ijを精度良く推定することが可能となる。粒子間距離ｒ_ijを精度良く推定することが可能となることにより、所定時間幅Δｔ内に有効半径Ｒｅ内に存在する流体粒子１３２,１６２から静止粒子１３１,１６１へ伝達される伝熱量を粒子間距離ｒ_ijに応じて精度良く重み付けすることが可能となる。これにより、気相（流体粒子）から固相（静止粒子）への熱量の総和（伝熱量）Ｈ^ｋ＋１ _static-fluidを精度良く推定することが可能となる。その結果、様々な伝熱形態において解析対象物の温度を精度良く推定することができる。

また、上記状態解析方法では、解析対象物が固体である場合、その固体を有限個の静止粒子１３１,１５４,１６１によって構成することとしている。この場合、気相／固相間の伝熱形態と同様に固相／固相間の伝熱形態を精度良く模擬することが可能となる。これにより固相／固相間の伝熱形態において解析対象物の温度を精度良く推定することが可能となる。

また、上記状態解析方法では、流体粒子１３２,１６２は、静止粒子１３１,１６１との間で粒子間距離ｒ_ijに応じた相互作用を可能にする有効半径Ｒｅを有していると良い。この場合有効半径Ｒｅによって多数の流体粒子の中から静止粒子に作用することが出来る流体粒子を限定することが可能となる。静止粒子に作用する流体粒子が限定されることにより、流体粒子から静止粒子へ伝達される伝熱量を算出する際の演算負荷を軽減することが可能となる。

更に、本発明にかかるシミュレーションプログラム（状態解析プログラム）１１２では、各繰り返し演算プロセス（例えば、図６のステップＳ４、ステップＳ５、ステップＳ６及びステップＳ７）は、ＣＰＵ１０１上で機能する専用のカーネル（図１６）を個別に有すると共に、そのカーネルにリンクされた一の繰り返し演算プロセスを専門に実行する専用の演算処理部（グリッド）（図１６）をＧＰＵ１１５上に個別に有する。これにより、コンピュータ装置１００（ホストコンピュータ）のＣＰＵ１０１の演算負荷を好適に軽減するのと同時に、汎用のコンピュータ上で複数の粒子についての複数の繰り返し演算プロセスを同時に実行することが可能となる。これにより、流体の解析に要する時間を大幅に短縮することが可能となり、最適な製品仕様を迅速に決定することができるようになる。

〔第３実施形態〕
次に、本発明の第３実施形態について説明する。
図２２は、図２の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）１２２と、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３とを示す説明図である。ここでは、先ず粒子に関する初期値入力（Initial particle（fluid）setting）として、粒子間距離（Particle Distance）、半径影響係数（Radius Influence Coefficient）、各粒子の初期位置（Initial Particle Position）、初期速度（Initial Velocity）、初期圧力（Initial Pressure）、初期温度（Initial Temperature）及び初期粘度（Initial Viscosity）等を入力する。なお、上記値は風洞試験結果（Wind tunnel experiment）に基づいて決定される。これにより、流体（粒子）の源となるソースが決定される。従って、後述する粒子の追跡演算（Source tracking）とは、このソースから流出する各粒子の流れ状態（位置、速度、圧力又は温度）をシミュレーションプログラム１１２によって所定時間幅毎に推定・模擬することを意味している。

さらに、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３として、使用される粒子源データ（bucket_dat）、バンパーパラメータ（Bumper_stl）、可動平面パラメータ（moving_plate_stl）、空気抵抗パラメータ（drag1_stl）、蒸発に関するパラメータ（evaporator_stl）、凝縮に関するパラメータ（condenser_stl）等を設定する。

次に、入力された粒子の初期値及びパラメータ並びに車両の特性パラメータを基に粒子の流れ状態を追跡する図２の演算プロセスＳＴ１－２の内容について説明する。

以降では、粒子追跡演算プロセス（Source tracking process）によって得られた流体（粒子）の速度（ｕ^ｋ＋１）及び圧力（Ｐ^ｋ＋１）に基づき、流体粒子（Fluid particles）中に置かれた物体（解析対象物）についての粒子抗力係数ＰＤ（Coefficient of drag force in Particle method）および粒子揚力係数ＰＬ（Coefficient of lift force in Particle method）を算出する粒子抵抗解析法について説明する。そして、本粒子抵抗解析法によって算出された流体粒子中における解析対象物の粒子抗力係数ＰＤ／粒子揚力係数ＰＬと、風洞試験によって得られた流体中における解析対象物の抗力係数ＣＤ／揚力係数ＣＬとの比較検証結果について説明する。

図２３は、本発明の粒子解析方法（Particle analysis）における粒子影響解析の概要を示す説明図である。この粒子影響解析は、流体中に置かれた解析対象物の空力特性（従来の抗力係数ＣＤ及び揚力係数ＣＬに相当する値）をシミュレーションプログラム１１２上で推定・模擬するためのプロセスである。本解析によって算出される「粒子抗力係数ＰＤ」は、流体中に置かれた解析対象物の抗力係数ＣＤに相当するものであり、詳細については図２４を参照しながら後述する。同様に、本解析によって算出される「粒子揚力係数ＰＬ」は、流体中に置かれた解析対象物の揚力係数ＣＬに相当するものであり、詳細については図２５を参照しながら後述する。

図２３（ａ）に示されるように、本解析は、全流体粒子（All particles）のうちで解析対象物の近傍に位置し、力学的作用を及ぼすことができる影響粒子（Influence particles）を選定する粒子選定プロセス（Particle Filtering Process）Ｐ１と、解析対象物が影響粒子（Influence particles）から受ける全抗力（Total drag force）を算出する抗力測定プロセス（Force Measurement Process）Ｐ２と、算出した全抗力（Total drag force）に基づいて粒子抗力係数ＰＤ及び粒子揚力係数ＰＬを算出するＰＤ／ＰＬ計算プロセス(PD/PL Calculation Process)Ｐ３とから主として構成される。なお、ＰＤ／ＰＬ計算プロセス（PD/PL Calculation Process）Ｐ３を実行する際に必要となる定数の入力は、図２のデータ入力プロセスＳＴ１－１において行われる。

図２３（ｂ）に示されるように、データ入力プロセスＳＴ１－１では定数（Constant factors）として、流体粒子の初期速度（Initial Velocity）、流体の密度（Density）、投影面積（Projection area）をユーザが入力する。「投影面積（Projection area）」とは、解析対象物の前面投影面積のことである。

図２３（ｃ）に示されるように、シミュレーションプログラム１１２を実行することによって算出された解析対象物の粒子抗力係数ＰＤ又は粒子揚力係数ＰＬについては、所定のコード及び所定のファイル形式で記憶装置１０８（図１）に保存されると共に、必要に応じて表示装置１０３（図１）において時系列毎に表示される。

図２４は、本解析において算出される粒子抗力係数ＰＤを示す説明図である。図２４（ａ）は、粒子抗力係数ＰＤにかかる各抗力及び上記影響粒子（Influence particles）１３２並びに影響粒子（Influence particles）１３２が作用する解析対象物１３０の作用面（ＳＴＬ三角平面）ｄＳをそれぞれ示している。また、図２４（ｂ）は、粒子抗力係数ＰＤを算出する手順の概要を示している。

図２４（ａ）に示されるように、ここでの解析対象物１３０は立方体（Cube）形状のＳＴＬモデル（STL model）としている。解析対象物１３０の各表面は、有限個のＳＴＬ三角平面（STL Triangles）ｄＳによって構成されている。このＳＴＬ三角平面（STL Triangle）ｄＳは、影響粒子（Influence particles）１３２が解析対象物１３０に作用する際の作用面となる。影響粒子（Influence particles）１３２は、解析対象物１３０の各面から所定距離（例えば、係数Ｋ*初期粒子間距離（initial particle distance）ｄ０）内に存在する流体粒子であり、解析対象物１３０に対する抗力源となる流体粒子である。

解析対象物１３０が影響粒子（Influence particles）１３２から受ける抗力（Drag）としては、影響粒子（Influence particles）１３２の圧力Ｐに起因する圧力抗力（Presure drag)Ｆ_Ｐと、影響粒子（Influence particles）１３２と解析対象物１３０との間のせん断力τ_ｗ（shear stress）に起因する粘性抗力（Viscous drag）Ｆ_Ｖとが挙げられる。

従って、粒子抗力係数Ｐ_Ｄは、解析対象物１３０が影響粒子（Influence particles）１３２の流れ場から受ける圧力抗力（Presure drag)Ｆ_Ｐと粘性抗力（Viscous drag）Ｆ_Ｖとの総和である全抗力（Total Drag force）Ｆ_Ｄの大きさを表す指標となる。なお、圧力抗力（Presure drag)Ｆ_Ｐ及び粘性抗力（Viscous drag）Ｆ_Ｖを算出する際の単位法線ベクトルｎの符号は、全抗力（Total Drag force）Ｆ_Ｄと同方向がプラスとなり、逆方向がマイナスとなる。

図２４（ｂ）に示されるように、圧力抗力Ｆ_Ｐは、圧力ＰがＳＴＬ三角平面（STL Triangle）ｄＳに作用する際の向きと大きさを示す微小面圧ベクトル（＝Ｐ・単位法線ベクトルｎ・ｄＳ）を全表面（Total surface area）、すなわち解析対象物１３０の全ＳＴＬ三角平面について積分することにより得られる。

粘性抗力Ｆ_Ｖについても同様に、せん断力τ_ｗがＳＴＬ三角平面ｄＳに作用する際の向きと大きさを示す微小せん断ベクトル（＝τ_ｗ・単位法線ベクトルｎ・ｄＳ）を全表面（Total surface area）、すなわち解析対象物１３０の全ＳＴＬ三角平面について積分することにより得られる。

算出した圧力抗力Ｆ_Ｐと粘性抗力Ｆ_Ｖとを合算して、全抗力（Total Drag force）Ｆ_Ｄを算出し、算出した全抗力（Total Drag force）Ｆ_Ｄを下記式２に代入して粒子抗力係数ＰＤを求める。なお、下記Ａは解析対象物１３０の前面投影面積（Projected frontal area）であり、下記ρは流体の密度（density of fluid）であり、下記Ｖは影響粒子（Influence particles）１３２の平均速度である。
式２：ＰＤ＝２*｜Ｆ_Ｄ｜/（Ａ*ρ*Ｖ^２）

図２５は、本解析において算出される粒子揚力係数ＰＬを示す説明図である。図２５（ａ）は、粒子揚力係数ＰＬにかかる揚力（Lift force）及び上記影響粒子（Influence particles）１３２並びに影響粒子（Influence particles）１３２が作用する解析対象物１３０の作用面（ＳＴＬ三角平面）ｄＳをそれぞれ示している。また、図２５（ｂ）は、粒子揚力係数ＰＬを算出する手順の概要を示している。また、説明の都合上、影響粒子（Influence particles）１３２の流れ方向（Flow direction）をＸ方向とし、それに直交する方向をＹ方向としている。

図２５（ａ）に示されるように、影響粒子（Influence particles）１３２の流れ場に置かれた解析対象物１３０は、流れ方向（Flow direction）に対し交差方向（transverse direction）に作用する揚力（Lift force）Ｆ_ｙを受けることになる。この揚力（Lift force）Ｆ_ｙは、圧力抗力Ｆ_ＰのＹ方向成分Ｆ_Ｐｙと、粘性抗力Ｆ_ＶのＹ方向成分Ｆ_Ｖｙとの総和に等しくなる。従って、粒子揚力係数Ｐ_Ｌは、影響粒子（Influence particles）１３２の流れ方向（Flow direction）と交差方向（transverse direction）に作用する力Ｆ_ｙの大きさを表す指標となる。なお、この揚力（Lift force）Ｆ_ｙを算出する際の単位法線ベクトルｎの符号は、Ｙ方向の正負方向と同じある。

図２５（ｂ）に示されるように、上記圧力抗力Ｆ_ＰのＹ方向成分Ｆ_Ｐｙを算出する。同様に、上記粘性抗力Ｆ_ＶのＹ方向成分Ｆ_Ｖｙを算出する。そして、圧力抗力のＹ方向成分Ｆ_Ｐｙと、粘性抗力のＹ方向成分Ｆ_Ｖｙとを合算した揚力Ｆ_ｙを算出する。そして、その揚力Ｆ_ｙを下記式３に代入して粒子揚力係数ＰＬを求める。なお、下記Ａは解析対象物１３０のＹ軸方向の投影面積（Projected frontal area）であり、下記ρは流体の密度（density of fluid）であり、下記Ｖは影響粒子（Influence particles）１３２の平均速度である。
式３：ＰＬ＝２*｜Ｆ_ｙ｜/（Ａ*ρ*Ｖ^２）

図２６は、影響粒子（Influence particles）１３２が解析対象物１３０に作用する際の作用面を示す説明図である。上記図２４及び図２５に示される粒子抗力係数ＰＤ及び粒子揚力係数ＰＬの算出（積分の計算）において、微小面圧ベクトル（＝Ｐ・単位法線ベクトルｎ・ｄＳ）又は微小せん断ベクトル（＝τ_ｗ・単位法線ベクトルｎ・ｄＳ）を全ＳＴＬ三角平面（＝ｄＳ_１＋ｄＳ_２＋・・・＋ｄＳ_ｎ）について積分している。つまり、図２６（ａ）に示されるように、積分計算における作用面（Surface area）として、全ＳＴＬ三角平面（＝ｄＳ_１＋ｄＳ_２＋・・・＋ｄＳ_ｎ）を用いている。

上記の手法以外にも、図２６（ｂ）に示されるように、積分計算における作用面（Surface area）として、影響粒子（Influence particles）１３２の解析対象物１３０への投影面積（Projected particle area）ｄｐ×ｄｐ（＝ｄｐ^２）を用いることも可能である。この場合、影響粒子（Influence particles）１３２が解析対象物１３０に及ぼす全抗力Ｆ_Ｄ及び揚力Ｆ_ｙをより精度良く算出することができるようになる。これにより、本粒子抵抗解析法によって算出される粒子抗力係数ＰＤ及び粒子揚力係数ＰＬを、従来の抗力係数ＣＤ又は揚力係数ＣＬと比較して、流体中における解析対象物１３０に作用する流体の影響をより正確に模擬（シミュレーション）した値に近付けることが可能となる。なお、「ｄｐ」は影響粒子（Influence particles）１３２の粒子径（Particle distance）に相当する。

図２７は、上記の粒子影響解析によって算出された粒子抗力係数ＰＤの検証結果を示すグラフである。粒子抗力係数ＰＤの検証は、風洞試験によって得られた解析対象物の実抗力係数ＣＤと、本粒子抵抗解析法によって得られた粒子抗力係数ＰＤとの間の誤差率（%ERROR）を評価することによって行われた。グラフの縦軸は各抗力係数（Drag Coefficient）を表し、横軸はレイノルズ数（Reynolds number）ＮＲを表している。また、曲線は風洞試験結果（実抗力係数ＣＤ）を表し、四角印（■）は本粒子抵抗解析法による粒子抗力係数ＰＤを表している。

また、解析対象物は直径１０ｍｍの球体（Sphere）とした。流体は、レイノルズ数（Reynolds number）ＮＲが１０^５未満では空気（Air）とし、レイノルズ数（Reynolds number）ＮＲが１０^５以上では水（Water）とした。

同図のグラフから明らかな通り、実抗力係数ＣＤと粒子抗力係数ＰＤとの間の誤差率は、小さいことが分かる。これは、本粒子抵抗解析法により得られた粒子抗力係数ＰＤは、流体にかかる実抗力係数ＣＤに近似することを示している。すなわち、流体の流れ場に置かれた解析対象物の抗力係数ＣＤは、本粒子影響解析によって算出される粒子抗力係数ＰＤによって評価可能であることを示している。

図２８は、上記の粒子影響解析によって算出された粒子揚力係数ＰＬの検証結果を示すグラフである。この粒子揚力係数ＰＬの検証は、風洞試験によって得られた解析対象物の実揚力係数ＣＬと、本粒子影響解析によって得られた粒子揚力係数ＰＬとの間の誤差率（%ERROR）を評価することによって行われた。グラフの縦軸は各揚力係数(Lift Coefficient)を表し、横軸は迎え角（Angle of attack）を表している。

図２８（ａ）に示されるように、解析対象物は、全米航空諮問委員会（＝National Advisory Committee for Aeronautics）によって開発された翼型（Airfoil）のうちで、「ＮＡＣＡ００１２」と通称されている翼型（以下、「本試験翼」という。）とした。本試験翼の寸法（dimension）は翼弦長が１２２ｍｍであり、最大翼厚が１４.５ｍｍである。

図２８（ｂ）に示されるように、本試験翼は翼幅を１０ｍｍとし、横寸法３７０ｍｍ×縦寸法３８０ｍｍの平板間に、翼前縁から平板間入口までの距離が９０ｍｍ、且つ翼後縁から平板底辺までの距離が１９０ｍｍとなるように設置されている。本試験翼の迎え角（Angle of attack）は、０°、２°、４°、６°、８°、及び１０°に変化させた。また、レイノルズ数（Reynolds number）ＮＲは１０００００と、流体粒子の初期速度（Initial velocity）は１２.８０４ｍｍ／ｓｅｃと、流体粒子の動粘度（Kinematic Viscosity）は１５.６２ｍｍ^２／ｓｅｃと、なるようにそれぞれ設定されている。

図２８（ｃ）は、揚力係数ＣＬと粒子揚力係数ＰＬとの比較結果を示す。菱形印（◆）が揚力係数ＣＬを表し、四角印（■）が粒子揚力係数ＰＬを表している。同図に示されるように、揚力係数ＣＬと粒子揚力係数ＰＬとの間の誤差率は、小さいことが分かる。これは、本粒子影響解析により得られた粒子揚力係数ＰＬは、揚力係数ＣＬに近似することを示している。すなわち、流体の流れ場に置かれた解析対象物の流れからの影響は、本粒子影響解析によって算出される粒子揚力係数ＰＬで評価が可能であることを示している。

図２９は、上記の粒子影響解析による車両の旋回時における空力特性（空気抵抗力）の解析を示す説明図である。なお、車両の速度は一定の１４４ｋｍ／ｈ（＝４０ｍ／ｓ）とし、旋回方向は左方向とし、旋回角度θは０°、５°、１０°、１５°及び２０°(degree)とした。また、車両の先端部の斜線部分は影響粒子１３２の集まりを表している。

図２９（ａ）に示されるように、車両の直進時（旋回角度θ＝０°）では、空気抵抗力（影響粒子１３２）は車両先端部において車幅方向に均等に作用している。つまり、車両先端部における空気抵抗力の作用領域は車幅方向に均等に分布している。

図２９（ｂ）から（ｅ）に示されるように、旋回角度θが大きくなるにつれて、車両先端部おける空気抵抗力の作用領域は左旋回の場合、車幅方向右側部分に集約されるようになる。つまり、車両が左旋回する場合、空気抵抗力は車両先端部において車幅方向右側部分に局所的に作用するようになる。これは、車両が左旋回する場合、旋回角度θが大きくなるにつれて空気抵抗力に起因した、車両を左方向に回転させようとする回転モーメントを車両が受けることを示している。

従来の抗力係数ＣＤ、揚力係数ＣＬを用いた解析では、流体の流れに対する対象物の投影面の形状及び投影面積が変化する場合に対しては正確な解析を行うことができなかった。これに対して、上記の粒子抗力係数ＰＤ及び粒子揚力係数ＰＬを用いた本実施形態の流体影響解析は、静止している対象物や進行方向（ベクトル）が変化しない移動をする対象物など、流体に対する投影面の形状及び投影面積が変化しない対象物の周辺の流れに対する解析だけでなく、車両が旋回しながら走行する場合など流体に対する投影面の形状及び投影面積が刻々と変化する対象物の周辺の流れに対する解析を行うことができる。したがって、動的な対象物の周辺の流れ状態（Aerodynamics）についても好適に模擬（解析）することが可能である。

以上の通り、本発明の粒子解析方法（Particle analysis）における粒子影響解析は、データ入力プロセスＳＴ１－１及び演算プロセスＳＴ１－２を備えるため、各流体粒子の時間幅Δｔ毎の位置、速度及び圧力を精度良く算出することができる。さらに、この粒子影響解析は、全流体粒子の内で解析対象物１３０に作用することができる影響粒子（Influence particles）１３２を選定する粒子選定プロセス（Particle Filtering Process）Ｐ１、影響粒子（Influence particles）１３２が解析対象物１３０に及ぼす全抗力Ｆ_Ｄ及び揚力Ｆ_ｙを算出する抗力測定プロセス（Force Measurement Process）Ｐ２、並びに全抗力Ｆ_Ｄ及び揚力Ｆ_ｙの各大きさに関する指標を算出するＰＤ／ＰＬ計算プロセス（PD/PL Calculation Process）Ｐ３を有する。これにより、流体中において解析対象物１３０が流体から受ける力の大きさの指標となる抗力係数ＣＤ及び揚力係数ＣＬを、本粒子影響解析により算出される粒子抗力係数ＰＤ及び粒子揚力係数ＰＬによって評価することが可能となる。これにより、解析対象物１３０の空力特性に関する製品仕様を迅速に決定することができるようになる。

また、本粒子解析方法の粒子選定プロセス（Particle Filtering Process）Ｐ１では、解析対象物１３０の各表面から所定距離の範囲内に位置する流体粒子を影響粒子（Influence particles）１３２として選定することが可能である。この場合、所定距離を変えることにより、算出される全抗力Ｆ_Ｄ及び揚力Ｆ_ｙ並びに上記粒子抗力係数ＰＤ及び粒子揚力係数ＰＬについての精度を向上させることが可能となる。

また、本粒子影響解析では、解析対象物１３０の形状データは、解析対象物１３０の表面を三角形で分割したＳＴＬモデル（STL model）とすることが可能である。この場合、解析対象物１３０の表面に接する接平面を容易に決定することができるため、解析対象物１３０の表面と各流体粒子との距離を容易に算出することができる。これにより、解析対象物に作用する影響粒子（Influence particles）１３２を容易に選定することが可能となる。

さらに、影響粒子（Influence particles）１３２が解析対象物１３０に作用する際の解析対象物の作用面ｄＳの法線ベクトルｎを容易に算出することができるため、解析対象物１３０が影響粒子（Influence particles）１３２から受ける全抗力Ｆ_Ｄ及び揚力Ｆ_ｙ、並びに上記粒子抗力係数ＰＤ及び粒子揚力係数ＰＬを容易に算出することが可能となる。

また、本粒子影響解析では、影響粒子（Influence particles）１３２が解析対象物１３０に作用する作用面ｄＳは、影響粒子（Influence particles）１３２の解析対象物１３０（ＳＴＬ三角平面）への投影面（ｄ_ｐ×ｄ_ｐ）とすることが可能である。影響粒子（Influence particles）１３２が解析対象物１３０に作用する作用面ｄＳが影響粒子（Influence particles）１３２の解析対象物１３０への投影面（ｄ_ｐ×ｄ_ｐ）である場合、影響粒子（Influence particles）１３２が解析対象物１３０に及ぼす全抗力Ｆ_Ｄ及び揚力Ｆ_ｙをより精度良く算出することができるようになる。これにより、本粒子影響解析によって算出される上記粒子抗力係数ＰＤ及び粒子揚力係数ＰＬを、流体における解析対象物１３０の抗力係数ＣＤ又は揚力係数ＣＬにより近似させることが可能となる。

更に、本発明にかかるシミュレーションプログラム（粒子解析プログラム）１１２では、各繰り返し演算ステップ（例えば、図６のステップＳ４、ステップＳ５、ステップＳ６及びステップＳ７）は、ＣＰＵ１０１上で機能する専用のカーネル（図１６）を個別に有すると共に、そのカーネルにリンクされた一の繰り返し演算ステップを専門に実行する専用の演算処理部（グリッド）（図１６）をＧＰＵ１１５上に個別に有する。これにより、コンピュータ装置１００（ホストコンピュータ）のＣＰＵ１０１の演算負荷を好適に軽減するのと同時に、汎用のコンピュータ上で複数の粒子についての複数の繰り返し演算ステップを同時に実行することが可能となる。これにより、流体の解析に要する時間を大幅に短縮することが可能となり、最適な製品仕様を迅速に決定することができるようになる。

〔第４実施形態〕
次に、本発明の第４実施形態について説明する。
図３０は、本発明の粒子解析方法（Particle Analysis）の全体の処理の流れを示すフローチャートである。また、図３１は、図３０の各処理で作成されるデータのファイルを模式的に示す図である。なお、本実施形態では、解析対象物である車両の走行中に該車両の外周面に沿って流れる空気（粒子）の流れから受ける力（揚力）の解析を行う場合について説明する。

図３０に示すように、本解析は、解析対象物である車両の解析用形状モデル（以下、「ＳＴＬモデル（STL model）」という。）の準備ならびに粒子及び車両に関する各種パラメータの入力（設定）を行うデータ入力プロセス（ＳＴ１）と、入力された各種パラメータに基づいてＧＰＵ１１５上でシミュレーションプログラム（solver）１１２を実行する演算プロセス（ＳＴ２）と、シミュレーションプログラム１１２の実行結果（演算結果）を所定のコード及び所定のファイル形式で記憶装置１０８等の周辺装置に出力するデータ出力プロセス（ＳＴ３）と、シミュレーションプログラム１１２の演算結果を基に流体（粒子）の流れが解析対象物に及ぼす力の状態を可視化する可視化プロセス（ＳＴ４）と、を有して構成される。以下、各プロセスについて説明する。

データ入力プロセス（ＳＴ１）は、ＳＴＬモデル作成プロセス（CAD model to STL）１２１と、流体（粒子）の初期状態量及び各種物性値を入力する粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）１２２と、解析対象物（車両）に関する各種特性パラメータを設定する解析対象物の特性パラメータ設定プロセス（Configuration setting for models）１２３と、から構成される。

図３２は、図３０のＳＴＬモデル作成プロセス（CAD model to STL）１２１を示す説明図である。図３２（ａ）は、ＣＡＤ等の描画作成ソフトによって作成された車両に関する三次元ＣＡＤ形状モデル（3D CAD model）である。図３２（ｂ）はＳＴＬ（STereoLithography）モデル（以下、「ＳＴＬモデル（STL model）」という。）である。このＳＴＬモデル（STL model）は、三次元ＣＡＤ形状モデル（3D CAD model）の表面（Surface）の一部又は全部を、図３２（ｃ）に示される三角形平面（以下、「ＳＴＬ三角平面（STL triangle）」という。）で分割することによって作成される。

従って、解析対象物である車両のＳＴＬモデル（STL model）は、ＳＴＬ三角平面（STL triangle）を基に構成されている。そのため、ＳＴＬモデルの表面に接する接平面の方程式を容易に決定することができ、その接平面によって粒子の追跡、可視化（着色化）が極めて容易となる。なお、１又は複数のＳＴＬ三角平面（STL triangle）を含むＳＴＬモデル（STL model）表面の一部分（Ａ部等）を、特にＳＴＬ部分平面（STL plate）と呼び、ＳＴＬ部分平面（STL plate）を含む接平面のことをＳＴＬ包絡平面（STL plane）と呼ぶことにする。上記で作成されたＳＴＬモデル（STL model）のデータは、図３１に示すＳＴＬファイル１５（Stl File）に格納されている。

図３３は、図３０の粒子パラメータ設定プロセス（Initial fluid（particle）property and parameter setting）１２２と、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３とを示す説明図である。ここでは、先ず粒子に関する初期値入力（Initial particle（fluid）setting）として、粒子間距離（Particle Distance）、半径影響係数（Radius Influence Coefficient）、各粒子の初期位置（Initial Particle Position）、初期速度（Initial Velocity）、初期圧力（Initial Pressure）、初期温度（Initial Temperature）及び初期粘度（Initial Viscosity）等を入力する。なお、上記値は風洞試験結果（wind tunnel experiment）に基づいて決定される。これにより、流体（粒子）の源となるソースが決定される。また、粒子に関するパラメータ（Parameter setting）として、初期粒子間距離（Initial Particle Distance）、自由表面の判定（Free Surface Threshold）、勾配モデル影響半径係数（Gradient model influence radius coefficient）、ラプラシアンモデル影響半径係数（Laplacian model influence radius coefficient）、温度計算用影響半径係数（Temperature calculation influence radius coefficient）、粒子壁面間距離オフセット（Wall Particle Distance Offset）などを含めてもよい。上記のデータは、図３１に示す粒子ファイル１１（Lift.csv file）及びパラメータファイル１３（Parameter file）に格納されている。

さらに、車両の特性パラメータ設定プロセス（Configuration setting for models）１２３で設定される車両の特性パラメータとして、車両の各車輪の接地点（接触点）の位置（座標）データ（.FRxyz, .FLxyz, .RRxyz, .RLxyz）や、車両の中心（重心）の位置（座標）データなどを設定する。上記の車両の特性パラメータは、図３１に示す構成ファイル１２（Config file or Configuration file）に格納されている。

次に、上記のデータ入力プロセスＳＴ１で入力された粒子の初期値及びパラメータ並びに車両の特性パラメータを用いて車両にかかる力を算出する図３０の演算プロセス（ＳＴ２）の内容について説明する。

図３４は、車両にかかる力の算出（CoFSolver）を行う演算プロセス（ＳＴ２）での処理の流れを示すフローチャートである。この演算プロセス（ＳＴ２）での処理には、粒子ファイル１１に格納されたデータ及びＳＴＬファイル１５に格納されたデータの読込（Read particle and stl）処理（ＳＴ２－１）、粒子選別（Filter particles）処理（ＳＴ２－２）、力の算出（Force calculation）処理（ＳＴ２－３）、力の中心算出（Center of Force calculation）処理（ＳＴ２－４）、各車輪にかかる力の算出（Force on each wheels）処理（ＳＴ２－５）の各処理が含まれている。以下、各処理について詳細に説明する。

粒子ファイル１１に格納されたデータ及びＳＴＬファイル１５に格納されたデータの読込（Read particle and stl）処理（ＳＴ２－１）では、構成ファイル１２（config.txt）に格納されているすべてのＳＴＬデータを読み込む（Read all .stl from config.txt）（ＳＴ２－１ａ）。また、粒子ファイル１１にリストされている各データから粒子の情報を読み込む（Read particles from each output.dat listed in .csv）（ＳＴ２－１ｂ）。また、この際、パラメータファイル１３、計算ゾーンファイル１４（Compute Zone File）及びＳＴＬファイル１５に格納されているデータも読み込む。

次に、粒子選別（Filter particles）処理（ＳＴ２－２）では、まず、ＳＴＬファイル１５に格納されている解析対象物である車両の各三角形の重心の面からの距離に基づいて粒子を選択する（Select particle based on distance from the face of centroid of each triangle）（ＳＴ２－２ａ）。そして、選択した粒子以外の他の粒子を除外する（Neglect other particles）（ＳＴ２－２ｂ）。また、車両の各三角形の平均圧力を格納する（Store an average pressure for each triangle）（ＳＴ２－２ｃ）。こうして、全流体粒子（All particles）のうち解析対象物である車両の近傍に位置し、力学的作用を及ぼすことができる影響粒子（Influence particles）を選定する。

力算出処理（Force calculation）（ＳＴ２－３）では、まず、解析対象物である車両が有する各三角形の領域（面積）を計算する（Calculate area of triangles）（ＳＴ２－３ａ）。また、各三角形にかかる（作用する）力を計算する（ＳＴ２－３ｂ）。そして、算出したこれらデータを別の粒子タイプまたは色として同じ出力データファイル（第一出力データファイル１６）に格納する（ＳＴ２－３ｃ）。

図３５は、力算出処理（ＳＴ２－３）を説明するための図である。力算出処理（ＳＴ２－３）では、図３５（ｂ）に示すように、各三角形について、粒子の選択（Particle selection）を行う。この粒子の選択では、（各三角形の）中心を算出（Centroid calculation）し、重心と各頂点との距離Ｌ１、Ｌ２、Ｌ３（Distances between centroid and each vertex）を算出し、これら距離Ｌ１、Ｌ２、Ｌ３の中で最大値を確認（Check for Maximum value among L1,L2,L3）し、影響半径は最大長+ 2.1（Influence radius is maximum length+2.1）とし、通常ベクトル（0<cosθ<1）方向の粒子のみを選択（Selecting particle only at direction of normal vector）する。また、図３５（ｃ）に示すように、力の計算（Force Calculation）として、平均圧力Average Pressure（P）、Ｘ方向分力Force（Fx, i）、Ｙ方向分力Force（Fy, i）、Ｚ方向分力Force（Fz, i）を算出する。こうして、最終的に図３５（ｄ）に示す合力（Resultant Force）Fx,Fy,Fz,Fが算出される。

図３４に戻り、力の中心算出処理（Center of Force calculation）（ＳＴ２－４）では、すべての（三角形の）重心の位置ベクトルを読み込む（Read the position vector of all centroids）（ＳＴ２－４ａ）。また、力ベクトルと位置ベクトルのクロス積を足し合わせる（Add the cross product of force and position vectors）（ＳＴ２－４ｂ）。そして、各力の大きさを加える（Add the magnitude of each force）（ＳＴ２－４ｃ）。さらに、クロス積を和で除して重心位置を求める（Divide the cross product by sum to find the centroids position）（ＳＴ２－４ｄ）。上記のデータを出力データファイルに別々の粒子タイプとして保存する（Store this data as separate particle type in output.dat）（ＳＴ２－４ｅ）。

図３６は、力の中心算出処理（Center of Force calculation）（ＳＴ２－４）で算出した力の中心の値を示す図である。同図に示すように、力の中心算出処理（Center of Force calculation）（ＳＴ２－４）では、i番目の三角形のＸ軸、Ｙ軸、Ｚ軸それぞれのモーメントであるMx, i, My, i, Mz, iと、Ｘ軸、Ｙ軸、Ｚ軸それぞれにおける力の中心であるCOFx, COFy, COFzとが算出される。ここで、Rx, i = i番目の三角形に対する重心ベクトルのX成分であり、Fx = x軸の合力の単位ベクトルであり、F =総合力の大きさである。

図３７は、車両にかかる力の釣り合い（Force Balance）を示す図である。各車輪にかかる力の算出処理（Force on each wheels）（ＳＴ２－５）では、構成ファイル１２から車輪の位置情報を読み込む（Read the wheel points from config）（ＳＴ２－５ａ）。また、力平衡法を用いて各車輪にかかる力を計算する（Calculate forces on each wheels using force balance method）（ＳＴ２－５ｂ）。図３７では、FRは右前輪接地点、FLは左前輪接地点、RRは右後輪接地点、RLは左後輪接地点である。また、COFは力の中心である。そして、dx=RR.x-COF.x、dy=FR.y-COF.y、dz=FR.z-COF.z、X=RR.x-FR.x、Y=FR.y-FL.yである。

図３８は、各車輪にかかる力の算出（Force on each wheels）処理（ＳＴ２－５）を説明するための図で、同図（ａ）は、処理の手順を示すフローチャートであり、同図（ｂ）は、力の算出手順を説明するための図である。図３８（ａ）に示すように、各車輪にかかる力の算出（Force on each wheels）処理（ＳＴ２－５）では、まず、リフト成分のみによって中間平面（x-z平面に平行な平面）に分配された力を計算する（Step1）。ここで、図３８（ｂ）に示すように、（右）前輪にかかる力FR lift （L1）及び（右）後輪にかかる力RR lift （L2）と、車両の中心にかかる力（Center of Force）CoF Lift=Total lift force（L_tot）との関係は、前輪と後輪との間の距離Length=lとし、CoFと後輪との間の距離をXとすると、
L1=L_tot*（x/l）
L2=L_tot*（l-x）/l
すなわち、
L_tot=L1+L2
となる。
次に、各車輪にかかる力を算出するために前平面（y-z平面に平行な平面）に分配された力を計算する（Step2）。こうして、車両の各車輪にかかる力が算出される。

図３１のデータ出力プロセス（ＳＴ３）では、上記の演算プロセス（ＳＴ２）で算出した車両の全体にかかる力及び各車輪にかかる力のデータを第一出力データファイル１６（Output data File from COF solver）及び第二出力データファイル１７（Output files）として出力する。

ここで、第一出力データファイル１６には、例えば、複数（ここでは４つ）の異なるタイプ（タイプ０～３）の粒子データが格納されている。タイプ０の粒子は、各三角形の重心位置（力位置）で生成される通常の粒子である。タイプ１の粒子は、車体（全体）の空力の中心に位置する粒子である。タイプ２の粒子は、三角形の中心に位置する粒子であり、この粒子は力の大きさがゼロの粒子である。タイプ３の粒子は、地面に接触する車輪上の点に位置する粒子である。これら４つの異なるタイプの粒子データは、例えば色分けなどされて区別されている。これら各タイプの粒子のデータのうち、タイプ１の粒子のデータによって、車両の全体にかかる力が算出される。また、タイプ３のデータによって、各車輪（の接地点）にかかる力が算出される。

このように、第一出力データファイル１６に格納されたデータによって、粒子から車両の各部に及ぼされる力の分布が定まる。このデータを元に車両の全体にかかる力及び各車輪にかかる力が算出される。算出されたデータは、第二出力データファイル１７に格納されている。

図３９は、車両に作用する力の算出結果を示す図で、（ａ）は、車両の各車輪にかかる力（Lift forces on each wheel）を示す図、（ｂ）は、車体全体に作用する力（Resultant Force acting on car body）を示す図である。同図は、演算プロセス（ＳＴ２）における算出結果を図３０の可視化プロセス（ＳＴ４）で可視化した結果を示している。このように、本実施形態のシミュレーションプログラム１１２によれば、車両の各車輪にかかる力（揚力）と車体の全体に作用する力（揚力）を算出して視覚化することができる。

図４０は、車両の各車輪に作用する力（揚力）の比較（Comparison of lift forces acting on wheel）を示すグラフで、横軸はタイムステップ（Time step）、縦軸は車両にかかる力（Force）の大きさである。同図のグラフには、車両の左前輪FL、右前輪FR、左後輪RL、右後輪RRそれぞれにかかる揚力の経時変化、すなわち車両の四輪それぞれにかかる揚力の分布状態の経過時間に応じた変化が示されている。

以上説明したように、本発明にかかる粒子解析方法によれば、解析対象物である車両に作用する流体（車両の周囲を流れる気体）を有限個の粒子で模擬し、該粒子の流れ又は状態を解析することで、流体の流れが解析対象物に対して及ぼす作用の影響を解析することで、解析対象物に作用する流体を模擬した有限個の粒子それぞれの流れ状態（例えば、所定時間幅毎の位置、速度及び圧力など）を精度良く算出することができる。そのうえで、解析対象物が影響粒子から受ける力である粒子力を算出する粒子力算出処理である力の算出（Force calculation）処理（ＳＴ２－３）と、粒子力によって解析対象物の全体に作用する力を算出する全体力算出処理である力の中心算出（Center of Force calculation）処理（ＳＴ２－４）と、解析対象物の各部に作用する力を算出する各部力算出処理である各車輪にかかる力の算出（Force on each wheels）処理（ＳＴ２－５）と、を有することで、流体中における解析対象物の流体特性（例えば、空力特性）を粒子法による流体粒子の解析結果を用いて精度良く推定することが可能となる。これにより、設計対象物の流体特性（空力特性）に関する製品仕様などをより迅速かつ的確に決定することができるようになる。

また、格子を用いる従来の解析手法では、流体の流れ中に有る解析対象物の全体に作用する力及び各部に作用する力の正確な解析を行うことができなかったところ、本発明の上記粒子解析方法によれば、影響粒子から解析対象物に作用する力を正確に算出できるので、解析対象物の全体及び各部に作用する力の解析を精度良く行うことができるようになる。したがって、流体の流れ中に有る解析対象物に作用する力をより正確に模擬（解析）することが可能である。

また、本粒子解析方法の影響粒子選定処理である粒子選別（Filter particles）処理（ＳＴ２－２）では、解析対象物である車両の各表面から所定距離の範囲内に位置する空気の粒子を影響粒子として選定することが可能である。この場合、所定距離を変えることにより、例えば、流体中における解析対象物の流体特性（例えば、空力特性、以下同じ。）に関する指標である上記流体が解析対象物に対して及ぼす作用による影響を反映した指標についての精度を向上させることが可能となる。

また、本粒子解析方法では、解析対象物である車両の形状データは、当該車両の表面を三角形で分割した形状モデルとしている。この場合、解析対象物の表面に接する接平面を容易に決定することができるため、解析対象物の表面と各粒子との距離を容易に算出することができる。これにより、解析対象物に作用を及ぼすことができる影響粒子を容易に選定することが可能となる。

更に、本実施形態のシミュレーションプログラム（粒子解析プログラム）１１２では、各繰り返し演算ステップ（例えば、図３０のステップＳＴ２）は、ＣＰＵ１０１上で機能する専用のカーネル（図１６）を個別に有すると共に、そのカーネルにリンクされた一の繰り返し演算ステップを専門に実行する専用の演算処理部（グリッド）（図１６）をＧＰＵ１１５上に個別に有する。これにより、コンピュータ装置１００（ホストコンピュータ）のＣＰＵ１０１の演算負荷を好適に軽減するのと同時に、汎用のコンピュータ上で複数の粒子についての複数の繰り返し演算ステップを同時に実行することが可能となる。これにより、流体の解析に要する時間を大幅に短縮することが可能となり、最適な製品仕様を迅速に決定することができるようになる。

以上、本発明の第４実施形態を説明したが、本発明は、上記実施形態に限定されるものではなく、特許請求の範囲、及び明細書と図面に記載された技術的思想の範囲内において種々の変形が可能である。例えば、上記実施形態では、本発明の解析対象物が車両（自動車）であり、当該車両の走行中に該車両の外周面に沿って流れる空気（粒子）の流れから受ける力（揚力）の解析を行う場合について説明したが、本発明の解析対象物は車両（自動車）には限られず、それ以外にも、例えば、飛行機や船舶など他の乗物類や、ソーラパネルや家屋などの建築物などを解析対象物とすることができ、それらの周囲を流れる空気や水などの流体の流れからの影響による力の計算（解析）が必要な設計全般に応用することが可能である。

Claims

対象モデルに沿って流れる流体の流れ又は状態を解析するための流体解析方法であって、
入力データとして、対象モデルの形状データと、前記対象モデルに沿って流れる流体を粒子に模した場合の当該粒子の特性に関する情報と、を入力するステップと、
前記入力データを用いて粒子法により前記粒子の流れ又は状態の解析を行うステップと、
前記流れ又は状態の解析の結果を用いて出力データを出力するステップと、
前記出力データを用いた結果を画像処理することで、前記粒子の流れを可視化した画像情報を生成するステップと、
前記粒子が到達するターゲット領域の座標情報と、選択された前記粒子の識別情報及び座標情報と、前記出力データとを入力するステップと、
前記ターゲット領域の座標情報を用いて前記粒子の識別情報のスキャニングを行うステップと、
前記出力データ内の粒子の色分けを行うステップと、
選択された粒子の識別情報と、前記出力データ内の色分けされた粒子のデータとを用いて、選択された粒子の前記ターゲット領域までの通過経路の追跡を行うステップと、を含み、
前記出力ファイル内の粒子の色分けを行うステップは、
任意に選択した平面状の焦点領域を前記粒子が通過したか否かに関する変数を作成するステップと、
前記焦点領域を含む平面の方程式を求めるステップと、
前記変数と前記方程式を用いて前記粒子が前記焦点領域を通過したか否かを判断するステップと、
前記焦点領域を通過した後の前記粒子の色を他の色に変更するステップと、を含むことを特徴とする流体解析方法。
前記粒子の特性に関する情報は、前記粒子の座標、温度、速度、加速度の少なくともいずれかに関する情報であることを特徴とする請求項１に記載の流体解析方法。
請求項１又は２に記載の流体解析方法をコンピュータに実行させるためのコンピュータプログラム。