ホーム > 特許ランキング > 住友重機械工業株式会社
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-12-01
(45)【発行日】2023-12-11
(54)【発明の名称】シミュレーション装置、及びプログラム
(51)【国際特許分類】
G01N 11/00 20060101AFI20231204BHJP
G06F 30/25 20200101ALI20231204BHJP
【ＦＩ】
G01N11/00 A
G06F30/25
【請求項の数】
3
(21)【出願番号】P 2020191834
(22)【出願日】2020-11-18
(65)【公開番号】P2022080643
(43)【公開日】2022-05-30
【審査請求日】2023-02-15
(73)【特許権者】
【識別番号】000002107
【氏名又は名称】住友重機械工業株式会社
(74)【代理人】
【識別番号】100105887
【弁理士】
【氏名又は名称】来山 幹雄
(72)【発明者】
【氏名】廣田 真人
【審査官】野田 華代
(56)【参考文献】
【文献】特開２０１８－１４７３１５（ＪＰ，Ａ）
【文献】特開２０１１－２５２７４８（ＪＰ，Ａ）
【文献】米国特許出願公開第２０１７／０１６１４１３（ＵＳ，Ａ１）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｎ １１／００－１３／０４
Ｇ０６Ｆ ３０／２５
(57)【特許請求の範囲】
【請求項1】
粒子法を用いて流体の流れを解析するシミュレーション装置であって、解析すべき流体を定義する情報、解析の初期条件と境界条件、解析対象の空間に配置される壁面境界の形状を定義する壁情報が入力される入力部と、
前記入力部に入力された情報に基づいて、流体を複数の流体粒子で表し、複数の流体粒子の運動を解析する処理部と
を備え、
前記処理部は、
複数の流体粒子のそれぞれの運動への壁の寄与を、前記壁面境界の形状、及び前記壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求め、
複数の流体粒子のそれぞれについて、求められた壁の寄与及び他の流体粒子の寄与に基づいて複数の流体粒子の運動を解析し、
複数の流体粒子の運動を解析するときの支配方程式として、流体の連続の式及び運動方程式を用い、
複数の流体粒子のそれぞれの運動への壁の寄与を求める際に、
計算対象の流体粒子から前記壁面境界に下した垂線の延長線上であって、前記壁面境界からの距離が流体粒子の半径に等しい位置に、仮想的な壁を配置し、
前記支配方程式における仮想的な壁の速度を、計算対象の流体粒子の速度、計算対象の流体粒子から前記壁面境界までの距離に基づいて決定するシミュレーション装置。
【請求項2】
前記連続の式における仮想的な壁の速度の、前記壁面境界に平行な成分である接線成分の向き及び大きさを、計算対象の流体粒子の速度の接線成分の向き及び大きさと同一にし、前記連続の式における仮想的な壁の速度の、前記壁面境界に垂直な成分である法線成分の向きを、計算対象の流体粒子の速度の法線成分の向きと反対にし、
前記連続の式における仮想的な壁の速度の法線成分の大きさを、計算対象の流体粒子が前記壁面境界に接しているときに、計算対象の流体粒子の速度の法線成分の大きさに等しくし、計算対象の流体粒子が前記壁面境界から遠ざかるにしたがって、計算対象の流体粒子から前記壁面境界までの距離に反比例して減少させ、
前記運動方程式における仮想的な壁の速度の接線成分の向きを、計算対象の流体粒子の速度の接線成分の向きと反対にし、
前記運動方程式における仮想的な壁の速度の接線成分の大きさを、計算対象の流体粒子が前記壁面境界に接しているときに、計算対象の流体粒子の速度の接線成分の大きさに等しくし、計算対象の流体粒子が前記壁面境界から遠ざかるにしたがって、計算対象の流体粒子から前記壁面境界までの距離に反比例して減少させ、
前記運動方程式における仮想的な壁の速度の法線成分の向き及び大きさを、計算対象の流体粒子の速度の法線成分の向き及び大きさと同一にする請求項１に記載のシミュレーション装置。
【請求項3】
粒子法を用いて流体の流れを解析する手順をコンピュータに実行させるプログラムであって、解析すべき流体を定義する情報、解析の初期条件と境界条件、解析対象の空間に配置される壁面境界の形状を定義する壁情報を取得する手順と、
取得された情報に基づいて、流体を複数の流体粒子で表し、複数の流体粒子の運動を解析する手順と、
複数の流体粒子のそれぞれの運動への壁の寄与を、前記壁面境界の形状、及び前記壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求める手順と、
複数の流体粒子のそれぞれについて、求められた壁の寄与及び他の流体粒子の寄与に基づいて複数の流体粒子の運動を解析する手順と
をコンピュータに実行させ、
複数の流体粒子の運動を解析するときの支配方程式として、流体の連続の式及び運動方程式を用い、
複数の流体粒子のそれぞれの運動への壁の寄与を求める際に、
計算対象の流体粒子から前記壁面境界に下した垂線の延長線上であって、前記壁面境界からの距離が流体粒子の半径に等しい位置に、仮想的な壁を配置し、
前記支配方程式における仮想的な壁の速度を、計算対象の流体粒子の速度、計算対象の流体粒子から前記壁面境界までの距離に基づいて決定するプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、シミュレーション装置、及びプログラムに関する。
【背景技術】
【0002】
流体の流れ場を粒子系の運動として近似して、流体の挙動を解析するシミュレーション方法が知られている。このシミュレーション方法は、粒子法（ＳＰＨ法）と呼ばれる。ＳＰＨ法では、流体は複数の粒子として表される。ＳＰＨ法において壁面境界条件を課すための手法として、壁面境界の内部に複数の仮想粒子を配置する手法が公知である（特許文献１、非特許文献１参照）。また、壁面境界を仮想粒子ではなくポリゴンで再現する手法が公知である（非特許文献２参照）。
【先行技術文献】
【特許文献】
【0003】
【文献】特開２０１６―１２５９３４
【非特許文献】
【0004】
【文献】S. Marrone et. al., "An accurate SPH modeling of viscous flows around bodies at low and moderate Reynolds numbers", Journal of Computational Physics, 245 (2013), pp. 456-475
【文献】原田隆宏、越塚誠一、「ＳＰＨにおける壁境界計算手法の改良」、情報処理学会論文誌、Ｖｏｌ．４８、Ｎｏ．４、１８３８～１８４６頁、２００７年
【発明の概要】
【発明が解決しようとする課題】
【0005】
複数の仮想粒子を配置することによって壁面境界を再現する手法では、壁面境界の形状が複雑な場合や、壁面境界が変形したり運動したりする場合に、仮想粒子を適切に配置することが困難である。
【0006】
壁面境界をポリゴンで再現する手法では、流体を表す複数の粒子（流体粒子）のそれぞれから壁面境界に下した垂線に対して垂直な面内に均一に複数の壁粒子が並んでいると仮定して、重み関数や重み関数の勾配を、流体粒子から壁までの距離の関数として計算する。このように、壁面境界の形状がほぼ平坦であることが前提とされており、曲率を持つ壁が存在するときには、解析の精度が低下してしまう。
【0007】
本発明の目的は、ＳＰＨ法を用いたシミュレーションにおいて、流れ場内の壁の幾何学的形状に制約されることなく、高精度の解析を行うことが可能なシミュレーション装置、及びプログラムを提供することである。
【課題を解決するための手段】
【0008】
本発明の一観点によると
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
解析すべき流体を定義する情報、解析の初期条件と境界条件、解析対象の空間に配置される壁面境界の形状を定義する壁情報が入力される入力部と、
前記入力部に入力された情報に基づいて、流体を複数の流体粒子で表し、複数の流体粒子の運動を解析する処理部と
を備え、
前記処理部は、
複数の流体粒子のそれぞれの運動への壁の寄与を、前記壁面境界の形状、及び前記壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求め、
複数の流体粒子のそれぞれについて、求められた壁の寄与及び他の流体粒子の寄与に基づいて複数の流体粒子の運動を解析し、
複数の流体粒子の運動を解析するときの支配方程式として、流体の連続の式及び運動方程式を用い、
複数の流体粒子のそれぞれの運動への壁の寄与を求める際に、
計算対象の流体粒子から前記壁面境界に下した垂線の延長線上であって、前記壁面境界からの距離が流体粒子の半径に等しい位置に、仮想的な壁を配置し、
前記支配方程式における仮想的な壁の速度を、計算対象の流体粒子の速度、計算対象の流体粒子から前記壁面境界までの距離に基づいて決定するシミュレーション装置が提供される。
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
解析すべき流体を定義する情報、解析の初期条件と境界条件、解析対象の空間に配置される壁面境界の形状を定義する壁情報が入力される入力部と、
前記入力部に入力された情報に基づいて、流体を複数の流体粒子で表し、複数の流体粒子の運動を解析する処理部と
を備え、
前記処理部は、
複数の流体粒子のそれぞれの運動への壁の寄与を、前記壁面境界の形状、及び前記壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求め、
複数の流体粒子のそれぞれについて、求められた壁の寄与及び他の流体粒子の寄与に基づいて複数の流体粒子の運動を解析し、
複数の流体粒子の運動を解析するときの支配方程式として、流体の連続の式及び運動方程式を用い、
複数の流体粒子のそれぞれの運動への壁の寄与を求める際に、
計算対象の流体粒子から前記壁面境界に下した垂線の延長線上であって、前記壁面境界からの距離が流体粒子の半径に等しい位置に、仮想的な壁を配置し、
前記支配方程式における仮想的な壁の速度を、計算対象の流体粒子の速度、計算対象の流体粒子から前記壁面境界までの距離に基づいて決定するシミュレーション装置が提供される。
【0010】
本発明の他の観点によると、
粒子法を用いて流体の流れを解析する手順をコンピュータに実行させるプログラムであって、
解析すべき流体を定義する情報、解析の初期条件と境界条件、解析対象の空間に配置される壁面境界の形状を定義する壁情報を取得する手順と、
取得された情報に基づいて、流体を複数の流体粒子で表し、複数の流体粒子の運動を解析する手順と、
複数の流体粒子のそれぞれの運動への壁の寄与を、前記壁面境界の形状、及び前記壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求める手順と、
複数の流体粒子のそれぞれについて、求められた壁の寄与及び他の流体粒子の寄与に基づいて複数の流体粒子の運動を解析する手順と
をコンピュータに実行させ、
複数の流体粒子の運動を解析するときの支配方程式として、流体の連続の式及び運動方程式を用い、
複数の流体粒子のそれぞれの運動への壁の寄与を求める際に、
計算対象の流体粒子から前記壁面境界に下した垂線の延長線上であって、前記壁面境界からの距離が流体粒子の半径に等しい位置に、仮想的な壁を配置し、
前記支配方程式における仮想的な壁の速度を、計算対象の流体粒子の速度、計算対象の流体粒子から前記壁面境界までの距離に基づいて決定するプログラムが提供される。
粒子法を用いて流体の流れを解析する手順をコンピュータに実行させるプログラムであって、
解析すべき流体を定義する情報、解析の初期条件と境界条件、解析対象の空間に配置される壁面境界の形状を定義する壁情報を取得する手順と、
取得された情報に基づいて、流体を複数の流体粒子で表し、複数の流体粒子の運動を解析する手順と、
複数の流体粒子のそれぞれの運動への壁の寄与を、前記壁面境界の形状、及び前記壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求める手順と、
複数の流体粒子のそれぞれについて、求められた壁の寄与及び他の流体粒子の寄与に基づいて複数の流体粒子の運動を解析する手順と
をコンピュータに実行させ、
複数の流体粒子の運動を解析するときの支配方程式として、流体の連続の式及び運動方程式を用い、
複数の流体粒子のそれぞれの運動への壁の寄与を求める際に、
計算対象の流体粒子から前記壁面境界に下した垂線の延長線上であって、前記壁面境界からの距離が流体粒子の半径に等しい位置に、仮想的な壁を配置し、
前記支配方程式における仮想的な壁の速度を、計算対象の流体粒子の速度、計算対象の流体粒子から前記壁面境界までの距離に基づいて決定するプログラムが提供される。
【発明の効果】
【0011】
複数の流体粒子のそれぞれの運動への壁の寄与を、壁面境界の形状、及び壁面境界の近傍に存在する複数の流体粒子の空間分布を用いて求めることにより、流れ場内の壁の幾何学的形状に制約されることなく、高精度の解析を行うことが可能になる。また、壁面境界を再現する複数の仮想粒子を配置する必要がないため、容易に壁面境界を持つ流れ場のシミュレーションを行うことができる。
【図面の簡単な説明】
【0012】
【発明を実施するための形態】
【0013】
本発明の実施例について説明する前に、図１～図３Ｂを参照して、ＳＰＨ法を用いた従来のシミュレーション方法について簡単に説明する。ＳＰＨ法においては、流体を表す複数の流体粒子のそれぞれの位置にカーネル関数を配置し、カーネル関数の重ね合わせとして流体変数の空間分布を表す。
【0014】
図１は、解析モデルの一例を示す斜視図である。直方体状の解析空間１０内に流体を表す複数の流体粒子２０を配置し、流体粒子２０の各々に作用する力を求め、流体の支配方程式を数値的に解くことにより、流体粒子２０の挙動を解析する。解析空間１０の中に壁面境界１１をもつ壁が配置されている。
【0015】
ＳＰＨ法で用いられる流体の代表的な支配方程式は以下の式で表される。
【数1】
【数2】
【数3】
【数4】
【0016】
ここで、ｍｉ及びρｉは、それぞれｉ番目の流体粒子２０の質量及び密度を表す。ｐｉ、ｖｉ、ｒｉは、それぞれｉ番目の流体粒子２０の圧力、速度ベクトル（以下、単に「速度」という場合がある。）、及び位置ベクトル（以下、単に「位置」という場合がある。）を表す。ｃ０は音速、ρ０は基準密度、μは粘性係数である。また、ｄは空間の次元数を表す。Ｗｉｊは、ｉ番目の流体粒子２０とｊ番目の流体粒子２０との間のカーネル関数である。∇ｉＷｉｊは、ｉ番目の流体粒子２０の位置におけるカーネル関数Ｗｉｊの導関数を表すベクトルである。
【0017】
カーネル関数Ｗｉｊは、ｉ番目の粒子とｊ番目の粒子との距離ｒｉｊのみの関数であり、例えば以下の関数を用いることができる。
ここで、ｈはカーネル幅を示し、例えば初期状態の平均粒子間隔と同じ程度の値とすることができる。流体粒子２０は直径ｈの球体とみなすことができる。
【数5】
【0018】
式（１）は流体の連続の式を離散化したものであり、式（２）は流体粒子２０が従う運動方程式である。式（２）の右辺第１項は圧力勾配による力に相当し、右辺第２項は流体の粘性による力に相当する。
【0019】
流れ場に壁面境界１１が存在する場合、従来は、壁面境界１１の内部に複数の仮想粒子を配置していた。
【0020】
図２は、壁面境界１１の近傍の流体粒子２０、及び壁面境界１１を再現する仮想粒子２１を示す模式図である。壁面境界１１の外側に複数の流体粒子２０が配置されており、内側に複数の仮想粒子２１が配置されている。図２において、仮想粒子２１にハッチングを付している。流体粒子２０と仮想粒子２１とが相互作用すると考え、式（１）及び式（２）を下記のように書き直す。
式（６）及び式（７）のシグマ記号に付されたΩは、壁面境界１１の内側に配置されている仮想粒子２１の集合を意味する。すなわち、式（６）の右辺第１項、式（７）の右辺第１項及び第２項は、壁面境界１１の内側の複数の仮想粒子２１からｉ番目の流体粒子２０への寄与量の合計を表しており、式（６）の右辺第２項、式（７）の右辺第３項及び第４項は、複数の流体粒子２０からｉ番目の流体粒子２０への寄与量の合計を表している。
【数6】
【数7】
【0021】
仮想粒子２１からの寄与量を計算する際には、境界条件に応じて仮想粒子２１の質量ｍｊ、密度ρｊ、圧力ｐｊ、及び速度ベクトルｖｊの値を設定する必要がある。壁面境界１１において滑りなし条件（壁面で流速ゼロ）を課す場合、質量ｍｊ、密度ρｊ、及び圧力ｐｊを以下の式で与える。
【数8】
【0022】
連続の式（６）の速度ベクトルｖｊを以下の式で与える。
運動方程式（７）の速度ベクトルｖｊを以下の式で与える。
ここで、τ及びｎは、それぞれｉ番目の流体粒子２０に対して距離が最小となる壁面境界１１の位置における壁面境界１１の接線方向の単位ベクトル及び法線方向の単位ベクトルである。ｓｉは、ｉ番目の流体粒子２０と壁面境界１１との最短の距離を表し、ｓｊは、ｊ番目の仮想粒子２１と壁面境界１１との最短の距離を表す。
【数9】
【数10】
【0023】
次に、図３Ａ及び図３Ｂを参照して、式（９）及び式（１０）の物理的意味について説明する。図３Ａ及び図３Ｂは、流体粒子２０、仮想粒子２１、及び壁面境界１１の位置関係、及び速度ベクトルを示す模式図である。図３Ａは、式（９）の関係を持つ速度ベクトルを表しており、図３Ｂは、式（１０）の関係を持つ速度ベクトルを表している。
【0024】
図３Ａに示すように、式（９）は、ｊ番目の仮想粒子２１の速度ベクトルｖｊの接線方向の成分の大きさ及び向きが、ｉ番目の流体粒子ｉの速度ベクトルｖｉの接線方向の成分の大きさ及び向きと等しいことを意味している。さらに、式（９）は、ｊ番目の仮想粒子２１の速度ベクトルｖｊの法線方向の成分の向きが、ｉ番目の流体粒子ｉの速度ベクトルｖｉの法線方向の成分の向きと反対であり、法線方向の成分の大きさが、ｉ番目の流体粒子ｉの速度ベクトルｖｉの法線方向の成分の大きさに（ｓｊ／ｓｉ）を乗じた大きさと等しいことを意味している。
【0025】
図３Ｂに示すように、式（１０）は、ｊ番目の仮想粒子２１の速度ベクトルｖｊの法線方向の成分の大きさ及び向きが、ｉ番目の流体粒子ｉの速度ベクトルｖｉの法線方向の成分の大きさ及び向きと等しいことを意味している。さらに、式（１０）は、ｊ番目の仮想粒子２１の速度ベクトルｖｊの接線方向の成分の向きが、ｉ番目の流体粒子ｉの速度ベクトルｖｉの接線方向の成分の向きと反対であり、接線方向の成分の大きさが、ｉ番目の流体粒子ｉの速度ベクトルｖｉの接線方向の成分の大きさに（ｓｊ／ｓｉ）を乗じた大きさと等しいことを意味している。
【0026】
次に、図４～図７を参照して本願発明の実施例について説明する。
図４は、本実施例で用いる解析モデルの流体粒子２０及び壁面境界１１の位置関係を示す模式図である。壁面境界１１の外側に複数の流体粒子２０が配置されているが、壁面境界１１の内側に複数の仮想粒子は配置されない。壁面境界１１からｉ番目の流体粒子２０への寄与量を求めるとき、ｉ番目の流体粒子２０の近傍の壁面境界１１内に、壁面境界１１を再現する仮想的な壁２２を配置する。仮想的な壁２２は、解析モデルの空間内に常時配置される壁ではなく、流体粒子２０のそれぞれへの壁面境界１１からの寄与量を計算するときに、流体粒子２０ごとに暫定的に配置される。
図４は、本実施例で用いる解析モデルの流体粒子２０及び壁面境界１１の位置関係を示す模式図である。壁面境界１１の外側に複数の流体粒子２０が配置されているが、壁面境界１１の内側に複数の仮想粒子は配置されない。壁面境界１１からｉ番目の流体粒子２０への寄与量を求めるとき、ｉ番目の流体粒子２０の近傍の壁面境界１１内に、壁面境界１１を再現する仮想的な壁２２を配置する。仮想的な壁２２は、解析モデルの空間内に常時配置される壁ではなく、流体粒子２０のそれぞれへの壁面境界１１からの寄与量を計算するときに、流体粒子２０ごとに暫定的に配置される。
【0027】
【0028】
ここで、ｍｗ、ρｗ、ｐｗ、ｖｗ、ｒｗは、それぞれ仮想的な壁２２に与えられる質量、密度、圧力、速度、及び位置である。
【0029】
位置ｒｉに依存する任意の物理量Ａに対するカーネル関数の導関数の補間公式から、以下の式が成立する。
ここで、Ａ（ｒｉ）は、物理量Ａの空間上の任意の位置ｒｉにおける値であり、Ａｉは、ｉ番目の粒子が持っている物理量Ａの値である。式（１３）は、位置ｒｉの近傍に存在するｊ番目の粒子に対して、物理量Ａｊにカーネル関数の導関数を掛けた和を計算することにより、位置ｒｉにおける物理量Ａ（ｒｉ）の導関数が得られることを表している。
【数13】
【0030】
式（１３）において物理量Ａに「１」を代入すると、「１」の導関数はゼロであるから、以下の式が導出される。
【数14】
【0031】
式（１４）は、以下のように書き直すことができる。
式（１５）の左辺を、仮想的な壁２２からの寄与量としてまとめると、以下の式が導き出される。
【数15】
【数16】
【0032】
式（１６）は、カーネル関数の導関数を計算する際の仮想的な壁２２からの寄与量を、ｉ番目の流体粒子２０の近傍に存在する複数の流体粒子２０の空間分布から見積もることができることを意味する。連続の式（１１）及び運動方程式（１２）のｍｗ、ρｗ、∇ｉＷｉｗを含む項を、式（１６）の右辺に置き換えると、連続の式（１１）及び運動方程式（１２）を数値的に解くために決定すべきパラメータは、圧力ｐｗ、速度ｖｗ、及び位置ｒｗのみとなる。
【0033】
圧力ｐｗを、式（８）と同様に以下の式で与える。
【数17】
【0034】
連続の式（１１）の速度ｖｗを、式（９）と同様に以下の式で与える。
運動方程式（１２）の速度ｖｗを、式（１０）と同様に以下の式で与える。
【数18】
運動方程式（１２）の速度ｖｗを、式（１０）と同様に以下の式で与える。
【数19】
【0035】
ここで、ｓｗは、仮想的な壁２２（図４）から壁面境界１１までの距離を表しており、仮想的な壁２２の位置ｒｗが決まれば、距離ｓｗも決まる。
【0036】
次に、仮想的な壁２２の位置ｒｗについて説明する。本実施例において、位置ｒｗを以下の式で与える。
ここで、ｈは流体粒子２０のカーネル幅、すなわち粒子直径を表す。式（２０）のカーネル幅ｈは、式（５）のカーネル幅ｈと同一である。
【数20】
【0037】
次に、図５Ａ及び図５Ｂを参照して、式（１８）～式（２０）の物理的な意味について説明する。図５Ａ及び図５Ｂは、壁面境界１１、及び仮想的な壁２２の位置関係、及び速度を示す模式図である。図５Ａは、式（１８）の関係を持つ速度を表しており、図５Ｂは、式（１９）の関係を持つ速度を表している。
【0038】
図５Ａ及び図５Ｂに示すように、連続の式（１１）及び運動方程式（１２）を数値的に解く際の仮想的な壁２２の位置ｒｗとして、計算対象のｉ番目の流体粒子２０から壁面境界１１に下した垂線の延長線上であって、壁面境界１１からの距離が流体粒子２０のカーネル幅ｈの１／２（すなわち半径）に等しい位置を採用する。このとき、距離ｓｗは、以下の式で与えられる。
【数21】
【0039】
連続の式（１１）及び運動方程式（１２）における仮想的な壁２２の速度ｖｗは、それぞれ式（１８）及び式（１９）を用いて、計算対象のｉ番目の流体粒子２０の速度ｖｉ、及びｉ番目の流体粒子２０から壁面境界１１までの距離ｓｉに基づいて決定する。
【0040】
例えば、図５Ａに示すように、連続の式（１１）における仮想的な壁２２の速度ｖｗの接線成分の向き及び大きさを、計算対象のｉ番目の流体粒子２０の速度ベクトルｖｉの接線成分の向き及び大きさと同一にする。連続の式（１１）における速度ベクトルｖｗの法線成分の向きを、速度ベクトルｖｉの法線成分の向きと反対にする。距離ｓｉがｈ／２に等しいとき、すなわちｉ番目の流体粒子２０が壁面境界１１に接しているときに、速度ｖｗの法線成分の大きさを速度ｖｉの法線成分の大きさと等しくする。ｉ番目の流体粒子２０が壁面境界１１から遠ざかるにしたがって、速度ｖｗの法線成分の大きさを距離ｓｉに反比例して減少させる。
【0041】
例えば、図５Ｂに示すように、運動方程式（１２）における仮想的な壁２２の速度ｖｗの接線成分の向きを、ｉ番目の流体粒子２０の速度ｖｉの接線成分の向きと反対にする。さらに、距離ｓｉがｈ／２に等しいとき、すなわちｉ番目の流体粒子２０が壁面境界１１に接しているときに、速度ｖｗの接線成分の大きさを速度ｖｉの接線成分の大きさに等しくする。ｉ番目の流体粒子２０が壁面境界１１から遠ざかるにしたがって、速度ｖｗの接線成分の大きさを距離ｓｉに反比例して減少させる。運動方程式（１２）における仮想的な壁２２の速度ｖｗの法線成分の向き及び大きさは、ｉ番目の流体粒子２０の速度ｖｉの法線成分の向き及び大きさと同一にする。
【0042】
運動方程式式（１２）の右辺第２項の速度の差ｖｗ－ｖｉの法線方向の成分がゼロになる。また、位置の差ｒｗ－ｒｉの接線方向の成分はゼロである。このため、両者の内積が常にゼロになってしまう。言い換えると、仮想的な壁２２からの粘性力の寄与がゼロになってしまう。これを回避するために、運動方程式（１２）の右辺第２項を、以下の式に置き換える。
【数22】
【0043】
図６は、本実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部３０、処理部３１、出力部３２、及び記憶部３３を含む。入力部３０から処理部３１にシミュレーション条件等が入力される。さらに、オペレータから入力部３０に各種指令（コマンド）等が入力される。入力部３０は、例えば通信装置、リムーバブルメディア読取装置、キーボード、ポインティングデバイス等で構成される。
【0044】
処理部３１は、入力されたシミュレーション条件及び指令に基づいてＳＰＨ法によるシミュレーションを実行する。さらに、シミュレーション結果を出力部３２に出力する。シミュレーション結果には、シミュレーション対象物である粒子系の粒子の状態、粒子系の物理量の時間的変化等を表す情報が含まれる。処理部３１は、例えばコンピュータの中央処理ユニット（ＣＰＵ）を含む。ＳＰＨ法によるシミュレーションをコンピュータに実行させるためのプログラムが、記憶部３３に記憶されている。出力部３２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。
【0045】
図７は、実施例によるシミュレーション方法の手順を示すフローチャートである。
まず、ユーザが入力部３０からシミュレーション条件として、シミュレーション対象の流体を定義する情報、初期条件、境界条件、壁情報等を入力する。処理部３１は、入力部３０を介してこれらのシミュレーション条件を取得する（ステップＳ１）。シミュレーション対象の流体を定義する情報は、流体の密度、粘度等の物性値を含む。初期条件は、複数の流体粒子２０を解析空間１０（図１）に配置するための情報、流体粒子２０の初期速度を指定する情報等を含む。境界条件は、解析空間１０の形状及び大きさを指定する情報を含む。壁情報は、解析空間１０の内部に配置される壁面境界１１の幾何学的形状及び位置を定義する情報を含む。
まず、ユーザが入力部３０からシミュレーション条件として、シミュレーション対象の流体を定義する情報、初期条件、境界条件、壁情報等を入力する。処理部３１は、入力部３０を介してこれらのシミュレーション条件を取得する（ステップＳ１）。シミュレーション対象の流体を定義する情報は、流体の密度、粘度等の物性値を含む。初期条件は、複数の流体粒子２０を解析空間１０（図１）に配置するための情報、流体粒子２０の初期速度を指定する情報等を含む。境界条件は、解析空間１０の形状及び大きさを指定する情報を含む。壁情報は、解析空間１０の内部に配置される壁面境界１１の幾何学的形状及び位置を定義する情報を含む。
【0046】
処理部３１は、入力された初期条件に基づいて、解析空間１０に複数の流体粒子２０（図１）を配置する（ステップＳ２）。さらに、配置した流体粒子２０に初期速度を与える。
【0047】
次に、処理部３１は、流体粒子２０のそれぞれについて、流体の支配方程式（１１）、（１２）を数値的に解くことによって、時間刻み幅が経過した次状態の流体粒子２０の速度ｖｉを求める（ステップＳ３）。このとき、式（１６）を支配方程式（１１）、（１２）に代入するとともに、仮想的な壁２２に関する位置ｒｗ、速度ｖｗ、及び圧力ｐｗに、式（１７）～式（２０）を用いて得られる値を与える。
【0048】
その後、ステップＳ３で求めた複数の流体粒子２０のそれぞれの速度ｖｉに基づいて、複数の流体粒子２０のそれぞれの位置ｒｉを更新する（ステップＳ４）。ステップＳ３とステップＳ４とを、計算終了条件を満たすまで繰り返す（ステップＳ５）ことにより、流体粒子２０の位置ｒｉを時間発展させる。
【0049】
計算が終了すると、処理部３１は出力部３２に計算結果を出力する（ステップＳ６）。例えば、処理部３１は、解析により求められた流体粒子２０の位置及び速度、及び入力された壁情報により定義される壁面境界１１の形状を認識可能に画像として出力部３２に表示させる。
【0050】
上記実施例によるシミュレーション方法により、精度の高い解析結果が得られることを確認するために、実際に実施例による方法と比較例による方法とを用いてシミュレーションを行った。次に、図８及び図９を参照して、このシミュレーションの結果について説明する。
【0051】
図８は、実施例によるシミュレーション方法の精度を確認するために行ったシミュレーションの解析モデルを示す概略図である。シミュレーションの解析空間を２次元とし、円柱周りの流れ場のシミュレーションを行った。解析空間をｘ方向に長い長方形とし、解析空間の中心に円柱を配置した。解析空間の長辺の長さを２４Ｄとし、短辺の長さを１２Ｄとし、円柱の直径をＤとした。流体粒子のカーネル幅ｈをＤ／１２として、複数の流体粒子を間隔ｈで等間隔に配置した。解析空間の一方の短辺を流入境界とし、流入境界から解析空間内に、一様の速度Ｕで流体粒子を流入させた。解析空間の長辺には周期境界条件を課した。流体粒子の初期条件として、各流体粒子に、一様のｘ方向の速度Ｕを与えた。
【0052】
一様の速度Ｕを代表速度とするレイノルズ数Ｒｅが、３．２、６．４、及び１２．８の３つの条件でシミュレーションを行った。初期状態の時刻ｔをゼロとし、無次元時刻ｔＵ／Ｄ＝１５まで時間発展させ、無次元時刻ｔＵ／Ｄ＝１０から１５までの間に円柱に作用する抗力係数の時間平均を測定した。
【0053】
図９は、シミュレーション結果から求めた円柱に作用する抗力係数とレイノルズ数との関係を示すグラフである。横軸はレイノルズ数Ｒｅを表し、縦軸は抗力係数を表す。グラフ中の丸記号は、本実施例によるシミュレーション方法を用いて得られた抗力係数を示し、三角記号は、比較例によるシミュレーション方法を用いて得られた抗力係数を示す。比較例においては、円柱を再現する複数の仮想粒子２１（図２）を配置して、支配方程式（６）、（７）を用いて解析を行った。
【0054】
図９に示すように、実施例によるシミュレーション方法で得られた解析結果は、比較例によるシミュレーション方法で得られた解析結果にほぼ一致している。実施例によるシミュレーション方法により、壁面境界を複数の仮想的な粒子で再現する比較例によるシミュレーション方法の解析精度と同程度の精度で、物体周りの流れ場のシミュレーションを行うことが可能であることが確認された。
【0055】
次に、比較例によるシミュレーション方法と比べた上記実施例の優れた効果について説明する。
【0056】
比較例によるシミュレーション方法を用いる場合、シミュレーションの前工程として、図２に示すように、オペレータが壁面境界１１に沿って複数の仮想粒子２１を配置しなければならない。この作業は、壁面境界１１が図８に示した円柱のように、数学的に記述できる形状である場合には比較的容易である。ところが、壁面境界１１の形状が複雑で数学的に記述できない場合には、壁面境界１１に沿って複数の仮想粒子２１（図２）を適切に配置することは困難である。実用的な解析に用いられる壁面境界１１の形状は、ほとんどの場合、数学的に記述できず、ＣＡＤデータ等で与えられる。解析対象の壁面境界１１の形状が変わるたびに、オペレータは、壁面境界１１を再現する仮想粒子２１を配置する作業を行わなければならない。
【0057】
上記実施例によるシミュレーション方法では、壁面境界１１を再現する複数の仮想粒子２１（図２）を配置する必要がない。上記実施例では、式（１６）に示したように、流体の支配方程式を解く際に、流体粒子２０への壁面境界１１からの寄与量を、他の流体粒子２０の空間分布から見積もっている。さらに、図５Ａ及び図５Ｂを参照して説明したように、仮想的な壁２２の位置ｒｗ及び速度ｖｗは、計算対象のｉ番目の流体粒子２０の位置及び壁面境界１１の形状から求まる。このように、壁面境界１１の位置及び形状を定義する情報を用いて、支配方程式（１１）、（１２）を解くことができる。このため、壁面境界１１の形状が複雑になっても、壁面境界１１の位置及び形状を、シミュレーションの計算に簡便に反映させることができる。
【0058】
上記実施例では、壁面境界１１を再現する複数の仮想粒子２１を配置する必要がないため、オペレータの作業効率の向上を図ることができる。また、流れ場に配置される物体の形状（壁面境界の形状）を最適化するために、物体の形状を多様に変化させて流れ場の解析を行う場合に、上記実施例によるシミュレーション方法を用いることにより、特に顕著な効果が得られる。
【0059】
本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【0060】
１０ 解析空間
１１ 壁面境界
２０ 流体粒子
２１ 壁面境界を再現する仮想粒子
２２ 壁面境界を再現する仮想的な物体
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部
１１ 壁面境界
２０ 流体粒子
２１ 壁面境界を再現する仮想粒子
２２ 壁面境界を再現する仮想的な物体
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】