特開 2023-117319 シミュレーション方法、シミュレーション装置、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2023117319

(43)【公開日】2023-08-23

(54)【発明の名称】シミュレーション方法、シミュレーション装置、及びプログラム

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20230816BHJP

【ＦＩ】

G16Z99/00

【審査請求】未請求

【請求項の数】5

【出願形態】ＯＬ

(21)【出願番号】P 2022019957

(22)【出願日】2022-02-10

(71)【出願人】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105887

【弁理士】

【氏名又は名称】来山 幹雄

(72)【発明者】

【氏名】廣田 真人

【テーマコード（参考）】

5L049

【Ｆターム（参考）】

5L049DD02

(57)【要約】

【課題】液滴の飛散等の現象を、ナビエストークス方程式を用いて解析可能で、かつ計算資源を削減することが可能なシミュレーション方法を提供する。
【解決手段】解析空間内の流れ場を、粒子法を用いて解析する際に、解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分する。２つの領域のそれぞれに複数の計算粒子を配置する。２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行う。２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解く。
【選択図】図３

【特許請求の範囲】

【請求項1】

解析空間内の流れ場を、粒子法を用いて解析するシミュレーション方法であって、
前記解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分し、
前記２つの領域のそれぞれに複数の計算粒子を配置し、
前記２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行い、
前記２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解くシミュレーション方法。

【請求項2】

前記境界面を横切って移動する計算粒子の位置及び個数に基づいて、計算粒子が流入した側の領域に計算粒子を追加し、追加された計算粒子に、流出側の領域の計算粒子の持つ物理量を、流入側の領域に適用される支配方程式で扱う物理量に交換した物理量を与える請求項１に記載のシミュレーション方法。

【請求項3】

前記２つの領域のそれぞれの領域から前記境界面を超えて相手側の領域に重なるバッファ領域を定義し、
前記バッファ領域のそれぞれに計算粒子を配置し、
前記２つの領域の境界面を含むある範囲内で物理量の交換を行う際に、前記バッファ領域のそれぞれに存在する計算粒子に付与される物理量を、前記バッファ領域に重なる相手側の領域の計算粒子に付与された物理量に基づいて計算する請求項１または２に記載のシミュレーション方法。

【請求項4】

粒子法を用いて液膜の流れを解析するシミュレーション装置であって、
解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分する情報を含むシミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて、前記２つの領域のそれぞれに計算粒子を配置し、前記２つの領域のそれぞれについて支配方程式を解くことにより、複数の計算粒子の位置を時間発展させる処理部と
を備え、
前記処理部は、
前記２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行い、
前記２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解くシミュレーション装置。

【請求項5】

粒子法を用いて液膜の流れを解析する手順をコンピュータに実行させるプログラムであって、
解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分する情報を含むシミュレーション条件を取得する手順と、
取得したシミュレーション条件に基づいて、前記２つの領域のそれぞれに計算粒子を配置する手順と、
前記２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行う手順と、
前記２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解くことにより、複数の計算粒子の位置を時間発展させる手順と
を実行させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。

【背景技術】

【0002】

流れ場を、スタッガードグリッドにより差分化して浅水方程式を解くことで、河川の氾濫のシミュレーションを行う技術が、下記の特許文献１に開示されている。特許文献１に開示されたシミュレーション方法では、粗い格子を用いて得られる浸水深分布を、細かい地形データで補間して、最終的に細かい浸水深分布を得る。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特許第４９７９３２２号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

河川の流れや液膜の流れのように、流路幅に対して水深の浅い流れは、壁面の形状等によってしぶきを発生し、液滴となって飛散することがある。この現象は、治水事業、排水装置及び噴霧化装置等の幅広い技術分野に関わる工学上重要な現象である。液膜の流れを二次元的に取り扱う浅水方程式では、しぶきの発生や液滴となって飛散する現象を解析することができない。

【0005】

ナビエストークス方程式に重力や表面張力等の適切な物理モデルを組み込むことで、しぶきの発生や液滴となって飛散する現象を解析することが可能である。ナビエストークス方程式を、粒子法を用いて離散化し、水深の浅い流れや薄い液膜の流れを解析するためには、計算粒子の直径を水深または液膜の厚さに対して十分小さい値にしなければならない。例えば、１００ｍｍ×１００ｍｍの領域に存在する厚さ０．１ｍｍの液膜の流れを解析するためには、計算粒子の直径を液膜の厚さの１／２０程度、すなわち０．００５ｍｍ程度にしなければならない。

【0006】

計算粒子の直径を０．００５ｍｍにしたとき、解析に必要な計算粒子の個数は、８×１０^９個程度になり、膨大な計算資源を必要とすることがわかる。本発明の目的は、液滴の飛散等の現象を、ナビエストークス方程式を用いて解析可能で、かつ計算資源を削減することが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。

【課題を解決するための手段】

【0007】

本発明の一観点によると、
解析空間内の流れ場を、粒子法を用いて解析するシミュレーション方法であって、
前記解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分し、
前記２つの領域のそれぞれに複数の計算粒子を配置し、
前記２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行い、
前記２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解くシミュレーション方法が提供される。

【0008】

本発明の他の観点によると、
粒子法を用いて液膜の流れを解析するシミュレーション装置であって、
解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分する情報を含むシミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて、前記２つの領域のそれぞれに計算粒子を配置し、前記２つの領域のそれぞれについて支配方程式を解くことにより、複数の計算粒子の位置を時間発展させる処理部と
を備え、
前記処理部は、
前記２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行い、
前記２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解くシミュレーション装置が提供される。

【0009】

本発明のさらに他の観点によると、
粒子法を用いて液膜の流れを解析する手順をコンピュータに実行させるプログラムであって、
解析空間を、支配方程式として浅水方程式を用いる領域とナビエストークス方程式を用いる領域との２つの領域に区分する情報を含むシミュレーション条件を取得する手順と、
取得したシミュレーション条件に基づいて、前記２つの領域のそれぞれに計算粒子を配置する手順と、
前記２つの領域の境界面を含むある範囲内の計算粒子に付与する物理量の値を、相手側の計算粒子に付与されている物理量の値から計算する物理量の交換を行う手順と、
前記２つの領域のそれぞれにおいて、物理量の交換によって物理量の値が計算された計算粒子の物理量の値を用いて支配方程式を解くことにより、複数の計算粒子の位置を時間発展させる手順と
を実行させるプログラムが提供される。

【発明の効果】

【0010】

ナビエストークス方程式を用いる領域で、液滴の飛散等の現象を解析することができる。一部の領域で浅水方程式を用いることにより、計算資源を削減することができる。

【図面の簡単な説明】

【0011】

【図1】図１は、液滴が発生する解析モデルの一例を示す模式図である。

【図2】図２は、解析空間に複数の計算粒子を配置した解析モデルを示す模式図である。

【図3】図３は、２つの領域間で物理用を交換する手法を説明するための模式図である。

【図4】図４は、液面高さＨを交換する手法を説明するための模式図である。

【図5】図５は、ナビエストークス方程式を適用する領域のバッファ領域内の壁面、計算粒子、及び液膜の表面を示す模式図である。

【図6】図６は、浅水方程式の領域からナビエストークス方程式の領域へ移動する計算粒子の制御を説明するための模式図である。

【図7】図７は、一実施例によるシミュレーション装置のブロック図である。

【図8】図８は、一実施例によるシミュレーション方法の手順を示すフローチャートである。

【図9】図９は、シミュレーション対象の解析空間を示す模式図である。

【図10】図１０Ａ及び図１０Ｂは、解析結果を示すグラフである。

【発明を実施するための形態】

【0012】

［流体の運動を記述するナビエストークス方程式］
本願発明の実施例について説明する前に、流体の運動を記述するナビエストークス方程式について説明する。

【0013】

流れ場の解析手段として、流れの支配方程式を離散化して数値的に解く手法が広く用いられている。ナビエストークス方程式に重力や表面張力等の適切な物理モデルを組み込むことで、水深の浅い流れ、その崩壊、及び液滴の発生現象を捉えることが可能である。例えば、重力と流体の表面張力を考慮する場合の、流体の圧縮性を無視した非圧縮性流れの連続の式及びナビエストークス方程式は次のように書ける。

【数1】

【数2】

ここで、ｕは速度ベクトル、ｔは時間、ρは密度、ｐは圧力、及びμは粘性係数を表す。ｇは重力ベクトル、σは表面張力係数、κは自由表面の曲率、及びｎは自由表面の単位法線ベクトルである。

【0014】

自由表面を有する流れ場の数値解析には、計算点が流れとともに移動するラグランジュ型解法（粒子法と表現される場合もある。）が適している。粒子法では、計算点の移動によって液面の形状及び液滴の形状等の自由表面の形状変化が自動的に表現される。このため、自由表面が大きく変化したり崩壊したりするような流れ場でも計算が破綻せず、安定した解析を行うことができる。流れ場のラグランジュ型解法として、例えばＳＰＨ（Smoothed Particle Hydrodynamics）法を用いることができる。

【0015】

［ＳＰＨ法の概要］
ＳＰＨ法は、連続体を粒子の集合を用いて離散化し、流れ場の時間発展を計算するラグランジュ型の数値計算手法である。粒子１個の質量分布として、連続かつ微分可能なカーネル関数を用いるとともに、カーネル関数の微分係数の重ね合わせを空間微分とする離散化モデルである。

【0016】

ＳＰＨ法では、粒子の質量ｍはカーネル関数Ｗによってその周囲に分配され、連続的に分布しており、密度場は、その重ね合わせとして与えられる。すなわち、ｉ番目の粒子の位置ｘ_ｉにおける密度ρ_ｉは、以下の式で表される。

【数3】

カーネル関数Ｗとして、例えば、以下の関数を用いることができる。以下の関数は、Ｃ２Ｗｅｎｄｌａｎｄカーネルと呼ばれる。

【数4】

【0017】

ここで、ｘ_ｉｊは位置ベクトルｘ_ｉと位置ベクトルｘ_ｊとの差の絶対値であり、ｑ＝ｘ_ｉｊ／ｈである。ｈは、平滑化長さであり、カーネル幅ともいわれる。カーネル幅ｈは、複数の粒子を同一の質量及び物性密度の球体に置き換えた等価球の直径に等しい。本明細書において、カーネル幅ｈを粒子の直径という場合がある。

【0018】

任意の位置ｘ_ｉにおける任意のスカラー量Ａ及びベクトル量Ａの値は、カーネル関数の重ね合わせによって以下の式で表すことができる。

【数5】

ここで、Ｖ_ｊ＝ｍ_ｊ／ρ_ｊである。下付き文字ｊは、ｊ番目の粒子の持つ物理量の値であることを意味する。

【0019】

任意のスカラー量Ａの勾配は、カーネル関数の微分係数の重ね合わせによって以下の式で表される。

【数6】

【0020】

ここで、∇Ｗ（ｘ_ｉｊ，ｈ）はカーネル関数の導関数であり、カーネル関数としてＣ２Ｗｅｎｄｌａｎｄカーネルを用いる場合、以下の式で表される。

【数7】

【0021】

任意のベクトル量Ａに対する発散は、式（６）と同様に、以下の式で表すことができる。

【数8】

【0022】

任意のスカラー量Ａに対する２次導関数は、例えば以下の公式を用いて算出することができる。

【数9】

【0023】

［ＳＰＨ法を用いたナビエストークス方程式の離散化手法］
次に、ＳＰＨ法を用いたナビエストークス方程式の離散化手法について説明する。
非圧縮性流れの支配方程式である式（１）及び式（２）を、ＳＰＨ法を用いて離散化すると、例えば次のように表すことができる。

【数10】

【数11】

【数12】

【数13】

ここで、下付き文字ｉ及びｊは、粒子に付与される粒子番号を表す。式（１２）は、流体が弱い圧縮性を持つという仮定により導出される方程式である。ｃは仮想的な音速、ρ_０は基準密度を表す。音速ｃの値は、例えば流れ場の代表速度の１０倍程度に設定することができる。基準密度ρ_０は、例えば作動流体の密度に設定することができる。式（１０）～式（１３）に示す支配方程式は、非圧縮性流れのＳＰＨ法による解法として広く用いられている。

【0024】

［ＳＰＨ法を用いた浅水方程式の離散化手法］
次に、ＳＰＨ法を用いた浅水方程式の離散化手法について説明する。浅水方程式を用いることにより、液膜流れを二次元的に解析することができる。

【0025】

浅水方程式は、以下の式で与えられる。

【数14】

【数15】

【数16】

【数17】

ここで、Ｈは液膜の厚さ（液面の高さ）、ｔは時間、Ｄ／Ｄｔは時間での全微分、ρは密度、ｕは速度ベクトル、ｐは圧力、μは粘性係数、ｇは重力ベクトル、ｎは液面が流れる壁面（曲面）の法線ベクトル、σは表面張力係数、κは液膜の曲率、ｇ_ｔは重力ベクトルの曲面に対する接線方向成分、θは液膜の接触角、ｎ_ｃｌは液膜の縁の法線ベクトルを表す。式（１６）の右辺第１項は重力による静水圧力を表し、第２項は液膜の表面張力によるラプラス圧を表す。式（１７）の右辺第１項は、液膜に作用する重力加速度、第２項は粘性摩擦力、第３項は液膜の縁に作用する表面張力を表す。

【0026】

液面の曲率κは、以下の式で表すことができる。

【数18】

ここで、下付き文字ｘ、ｙは、対応する空間方向の導関数を表す。例えば、Ｈ_ｘ、Ｈ_ｘｘは、以下の式で表される。

【数19】

【0027】

液膜の厚さＨが十分薄い場合、その１次導関数Ｈ_ｘ、Ｈ_ｙは１より十分小さくなるため、Ｈ_ｘの２乗の項及びＨ_ｙの２乗の項は無視することができる。したがって、式（１８）は、以下のように簡略化することができる。

【数20】

【0028】

ＳＰＨ法を用いて浅水方程式を離散化すると、次のように書ける。

【数21】

【数22】

【数23】

【数24】

【数25】

【0029】

浅水方程式は、液膜の厚さＨが、計算粒子に付与される物理量として表現される二次元的な方程式であるため、支配方程式の離散化の際に、液膜の厚さを基準として計算粒子の直径を設定する必要がない。このため、１００ｍｍ×１００ｍｍの領域に存在する厚さ０．１ｍｍの液膜を、浅水方程式を用いて離散化する場合、計算粒子の直径は、例えば液膜の厚さより大きな１ｍｍ程度に設定することができる。このとき、計算粒子の個数は１×１０^４個程度になる。これに対してナビエストークス方程式を離散化する場合、計算粒子の直径を液膜の厚さに比べて十分小さくする必要がある。このため、計算粒子の個数が多くなってしまう。浅水方程式を用いることにより、計算粒子の個数が大幅に減り、計算資源を削減することができる。

【0030】

［浅水方程式を用いる場合の課題］
ナビエストークス方程式に代えて浅水方程式を用いることにより、計算資源を削減することができる。ところが、浅水方程式では、液面からしぶきが発生し、液滴となって飛散する現象は、浅水方程式の前提条件から逸脱しているため、浅水方程式を用いて解くことはできず、ナビエストークス方程式を用いなければならない。

【0031】

次に、図１を参照して、液滴が発生する解析モデルの一例について説明する。
図１は、液滴が発生する解析モデルの一例を示す模式図である。水平面に対して傾斜した壁面１０が液膜１１で覆われている。液膜１１は、重力によって斜め下方に向かって流れ、壁面１０の下端に達した液膜１１は、液滴１２となって落下する。

【0032】

落下する液滴１２を含む流れ場を解析するためには、浅水方程式を用いることができず、ナビエストークス方程式を用いなければならない。流れ場全域を、ナビエストークス方程式を離散化して解くことは、有限の計算資源の観点から困難である。以下に説明する実施例では、流れ場の解析において、解析空間の液滴１２が発生しないことが予想できる領域では浅水方程式を用い、液滴１２が発生すると予想される領域ではナビエストークス方程式を用いる。

【0033】

図１に示した例では、解析空間を、傾斜した壁面１０の下端の近傍及び液滴が落下する領域Ω_ｎｓと、下端からある程度離れた壁面１０を含む領域Ω_ｓｗとに区分する。領域Ω_ｓｗでは、浅水方程式を用いて流れ場の解析を行い、領域Ω_ｎｓでは、ナビエストークス方程式を用いて液滴１２を含む流れ場の解析を行う。

【0034】

［領域間での物理量の交換］
次に、図２を参照して、解析空間を２つの領域に区分した場合の計算粒子の分布について説明する。

【0035】

図２は、解析空間Ωに複数の計算粒子２０を配置した解析モデルを示す模式図である。壁面１０が液膜で覆われている。壁面１０をｘｚ面とし、壁面１０の法線方向をｙ方向とするｘｙｚ直交座標系を定義する。解析対象の流体はｘ軸の正の向きに流れている。解析空間Ωを、ｘ方向に２つの領域Ω_ｓｗとΩ_ｎｓとに区分する。一方の領域Ω_ｓｗには浅水方程式を適用し、他方の領域Ω_ｎｓにはナビエストークス方程式を適用する。

【0036】

浅水方程式が適用される領域Ω_ｓｗでは、複数の計算粒子２０が壁面１０に沿って二次元的に配置される。すなわち、高さ方向（ｙ方向）には、計算粒子が分布しない。ナビエストークス方程式が適用される領域Ω_ｎｓでは、複数の計算粒子２０が三次元的に配置される。さらに、浅水方程式が適用される領域Ω_ｓｗ内の計算粒子２０に付与される物理量（浅水方程式が取り扱う物理量）と、ナビエストークス方程式が適用される領域Ω_ｎｓ内の計算粒子２０に付与される物理量（ナビエストークス方程式が取り扱う物理量）とは異なる。このため、２つの領域Ω_ｓｗとΩ_ｎｓとの境界面において、物理量の交換を行わなければならない。

【0037】

次に、図３を参照して領域間で物理量を交換する手法について説明する。
図３は、２つの領域間で物理量を交換する手法を説明するための模式図である。解析空間Ωが、境界面ＢＰで２つの領域Ω_ｓｗとΩ_ｎｓとに区分されている。図３では、２つの領域Ω_ｓｗとΩ_ｎｓとを、液膜の高さ方向（ｙ方向）に関して異なる位置に記載している。

【0038】

一方の領域Ω_ｓｗから境界面ＢＰを越えて他方の領域Ω_ｎｓ内まで広がり、他方の領域Ω_ｎｓに重なるバッファ領域Ω_ｓｗｂを定義する。バッファ領域Ω_ｓｗｂ内に、領域Ω_ｓｗ内の計算粒子２０と同じ大きさの複数の計算粒子２０を配置する。もう一方の領域Ω_ｎｓについても同様に、バッファ領域Ω_ｎｓｂを定義し、バッファ領域Ω_ｎｓｂ内に、領域Ω_ｎｓ内の計算粒子２０と同じ大きさの計算粒子２０を配置する。図３において、領域Ωｓｗ、Ωｎｓ内の計算粒子２０に相対的に濃いハッチングを付し、バッファ領域Ω_ｓｗｂ、Ω_ｎｓｂ内の計算粒子２０に相対的に淡いハッチングを付している。バッファ領域Ω_ｓｗｂ、Ω_ｎｓｂのそれぞれのｘ方向の寸法は、そのバッファ領域Ω_ｓｗｂ、Ω_ｎｓｂ内の計算粒子２０の直径の４倍以上、例えば約４倍とする。

【0039】

一方のバッファ領域Ω_ｓｗｂ内では、領域Ω_ｓｗに適用される支配方程式である浅水方程式を解かず、相手側の領域Ω_ｎｓのうちバッファ領域Ω_ｓｗｂと重なる領域内の計算粒子２０に付与されている物理量に基づいて求められた値を、バッファ領域Ω_ｓｗｂ内の計算粒子２０の物理量として与える。同様に、他方のバッファ領域Ω_ｎｓｂ内では、領域Ω_ｎｓに適用される支配方程式であるナビエストークス方程式を解かず、相手側の領域Ω_ｓｗのうちバッファ領域Ω_ｎｓｂと重なる領域内の計算粒子２０に付与されている物理量に基づいて求められた値を、バッファ領域Ω_ｎｓｂ内の計算粒子２０の物理量として与える。

【0040】

一方の領域のバッファ領域内の計算粒子２０の物理量として、相手側の領域の計算粒子に付与されている物理量に基づいて計算した値を与えることを、物理量の交換という。図３において、物理量の交換を行う向きを白抜き矢印で表している。このように、２つの領域Ω_ｓｗとΩ_ｎｓとの境界面ＢＰを含むある範囲内で物理量の交換を行う。流体の圧縮性を無視できる場合、支配方程式を解くために必要な物理量は、液面の高さＨ、速度ベクトルｕ、及び圧力ｐである。したがって、領域間でこれらの物理量を交換しなければならない。以下、これらの物理量を効果する手法について説明する。

【0041】

［液面高さＨの交換］
次に、図４を参照して、液面高さＨを交換する手法について説明する。
図４は、液面高さＨを交換する手法を説明するための模式図である。図４に示すように、浅水方程式を解く領域Ω_ｓｗでは、液面の高さＨが式（２１）を用いて計算されるため、高さ方向（ｙ方向）の粒子分布は存在しない。領域Ω_ｓｗ内の計算粒子の直径をｈ_ｓｗと標記する。

【0042】

ナビエストークス方程式を解く領域Ω_ｎｓでは、高さ方向にも計算粒子２０が分布しており、液膜の表面を構成している計算粒子２０Ｓのｙ方向の位置によって、液面の高さが与えられる。図４において、液膜の表面を構成している計算粒子２０Ｓにハッチングを付している。領域Ω_ｎｓ内の計算粒子２０の直径をｈ_ｎｓと標記する。

【0043】

浅水方程式を解く領域Ω_ｓｗのバッファ領域Ω_ｓｗｂ内のｊ番目の計算粒子２０の位置において、ｊ番目の計算粒子２０に付与される液面の高さＨ_ｊは以下の式で計算することができる。

【数26】

式（２６）においてｉは、領域Ω_ｎｓ及びバッファ領域Ω_ｎｓｂ内の計算粒子２０に付された通し番号を表し、ｊは、領域Ω_ｓｗ及びバッファ領域Ω_ｓｗｂ内の計算粒子２０に付された通し番号を表す。

【0044】

式（２６）の右辺のΣ記号は、ナビエストークス方程式を解く領域Ω_ｎｓ及びそのバッファ領域Ω_ｎｓｂの流体の表面を構成する計算粒子２０Ｓについて合計することを意味する。式（２６）のｙ_ｉは、ｉ番目の計算粒子２０Ｓの高さ方向の位置を表す。（ｘ_ｉ－ｘ_ｊ）_ｘｚは、バッファ領域Ω_ｓｗｂ内のｊ番目の計算粒子２０と、相手側の領域Ω_ｎｓ及びバッファ領域Ω_ｎｓｂ内のｉ番目の計算粒子２０との間の距離ベクトルをｘｚ面に射影したベクトルを表す。Ｗ_２ｄは、二次元的なカーネル関数であり、Ｗ_２ｄとして、例えば浅水方程式を解く領域Ω_ｓｗで用いられるカーネル関数と同じ関数を選ぶことができる。

【0045】

バッファ領域Ω_ｎｓｂ内のｉ番目の計算粒子２０の位置において、液面の高さＨ_ｉは、以下の式を用いて算出することができる。

【数27】

式（２７）の右辺のΣ記号は、領域Ω_ｓｗ及びバッファ領域Ω_ｓｗｂ内の計算粒子２０について合計することを意味する。

【0046】

式（２７）で算出される値を用いてバッファ領域Ω_ｎｓｂ内の表面の計算粒子２０Ｓの高さを変化させると、バッファ領域Ω_ｎｓｂ内の質量保存則に不整合を生じ、計算が不安定になることがある。バッファ領域Ω_ｎｓｂ内の表面の計算粒子２０Ｓの高さは、バッファ領域Ω_ｎｓｂ内に流入する流れの高さ及び速度によって自動的に決まるため、バッファ領域Ω_ｎｓｂ内の計算粒子２０及び２０Ｓの高さを、式（２７）を用いて変化させる必要はない。なお、式（２７）は、後述の計算粒子の流入量の制御を行う際に用いる。

【0047】

［速度ベクトルｕの交換］
次に、速度ベクトルｕを交換する手法について説明する。
浅水方程式を解く場合には、高さ方向の速度分布は前提条件によって決まる。例えば、高さ方向の平均流速ベクトルをｕと標記すると、速度ベクトルの高さ方向の分布ｕ（ｙ）として、次の放物線分布を用いることができる。

【数28】

ここで、ｙは高さ方向の位置を表す。式（２８）の速度ベクトルの分布を用いたとき、式（１７）の右辺の第２項の粘性係数μに関する項が得られる。

【0048】

ナビエストークス方程式を解く場合には、高さ方向の速度分布は、各計算粒子２０に与えられている速度ベクトルの分布によって決まる。

【0049】

浅水方程式の計算粒子２０が持つ速度ベクトルは、流れの高さ方向の速度ベクトルの平均値であるから、バッファ領域Ω_ｓｗｂ（図３）内に存在する計算粒子２０のそれぞれの速度ベクトルは、その計算粒子２０の位置において、相手側の領域Ω_ｎｓ内の計算粒子２０が持つ速度ベクトルの高さ方向の平均値を計算することにより求めることができる。バッファ領域Ω_ｓｗｂ（図３）内に存在する計算粒子２０のそれぞれの速度ベクトルｕ_ｊは、例えば、以下の式で計算することができる。

【数29】

式（２９）の右辺のΣ記号は、領域Ω_ｎｓ内の計算粒子２０について合計することを意味する。

【0050】

次に、ナビエストークス方程式を解く領域Ω_ｎｓのバッファ領域Ω_ｎｓｂ内の計算粒子２０の速度ベクトルｕの値を求める方法について説明する。計算粒子２０のそれぞれのｘｚ面内の位置における速度ベクトルを、相手側の領域Ω_ｓｗ内の計算粒子２０の速度ベクトルから算出し、式（２８）を用いて、バッファ領域Ω_ｎｓｂ内の計算粒子２０のｙ方向の位置に応じた速度ベクトルに補正すればよい。例えば、バッファ領域Ω_ｎｓｂ内の位置（ｘ_ｉ，ｙ_ｉ，ｚ_ｉ）の計算粒子２０の速度ベクトルｕ_ｉは、以下の式で計算することができる。

【数30】

式（３０）の右辺のΣ記号は、領域Ω_ｓｗ内の計算粒子２０について合計することを意味する。

【0051】

［圧力ｐの交換］
次に、図５を参照して圧力ｐを交換する手法について説明する。
浅水方程式を解く場合には、表面張力に起因する圧力（式（２３）の右辺第２項）が液膜の表面に作用しており、それに加えて、重力に起因する圧力（式（２３）の右辺第１項）が静水圧力として重力方向に分布していると考えることができる。これらの値は、式（２６）によって得られるバッファ領域Ω_ｓｗｂ内の計算粒子２０が持つ高さＨ_ｊの値を用いて算出することができる。

【0052】

図５は、バッファ領域Ω_ｎｓｂ内の壁面１０、計算粒子２０、及び液膜の表面１３を示す模式図である。傾斜した壁面１０が、複数の計算粒子２０で表される液膜１１で覆われている。壁面１０に平行な面をｘｚ面とするｘｙｚ直交座標系を定義する。計算粒子２０に重力加速度ｇが作用する。ｉ番目の計算粒子２０の位置（ｘ_ｉ，ｙ_ｉ，ｚ_ｉ）における壁面１０の単位法線ベクトルをｎ_ｉと標記する。液膜の表面１３の高さＨ_ｉは、その位置のｙ座標で定義される。

【0053】

ナビエストークス方程式を解く場合には、高さ方向の圧力分布は、計算粒子２０のそれぞれに付与されている圧力の分布によって決まる。例えば液膜の表面１３における圧力をゼロとしてバッファ領域Ω_ｎｓｂ内の計算粒子２０に付与される圧力ｐ_ｉは以下の式で計算することができる。

【数31】

ここで、ρ_０は式（１２）の基準密度を表す。Ｈ_ｉは、液面１３の高さであり、式（２６）のｘ_ｊに計算粒子２０の位置を代入することによって計算することができる。

【0054】

［浅水方程式の領域からナビエストークス方程式の領域へ移動する計算粒子の制御］
次に、図６を参照して、浅水方程式を適用する領域Ω_ｓｗからナビエストークス方程式を適用する領域Ω_ｎｓへ移動する計算粒子２０の制御について説明する。

【0055】

図６は、浅水方程式の領域からナビエストークス方程式の領域へ移動する計算粒子の制御を説明するための模式図である。図３に示した模式図と同様に、領域Ω_ｓｗに接続されたバッファ領域Ω_ｓｗｂが定義され、領域Ω_ｎｓに接続されたバッファ領域Ω_ｎｓｂが定義されている。領域Ω_ｓｗから領域Ω_ｎｓへの流体の流れがある場合、その速度に応じてバッファ領域Ω_ｎｓｂの端面に計算粒子２０を追加する必要がある。

【0056】

まず、バッファ領域Ω_ｎｓｂの端面に複数の格子３０を配置する。ここで、「端面」とは、領域Ω_ｎｓとの境界面以外の表面のうち、ｘ軸と交差する面を意味する。複数の格子３０のそれぞれを、計算粒子２０を追加する候補位置とする。

【0057】

次に、格子３０の配置方法の一例について説明する。まず、解析空間Ω全体に、格子の一辺の長さが、領域Ω_ｎｓ内の計算粒子２０の直径ｈ_ｎｓとほぼ等しい直交等間隔格子を張る。この複数の格子のうちバッファ領域Ω_ｎｓｂの内部にあって、バッファ領域Ω_ｎｓｂの端面からの距離が２^１／２×ｈ_ｎｓよりも近い位置に中心点を持つ格子を、バッファ領域Ω_ｎｓｂの端面に配置する複数の格子３０として選択する。

【0058】

複数の格子３０に通し番号ｋを付す。それぞれの格子３０において、以下の２つの式を満たす場合、その格子からバッファ領域Ω_ｎｓｂ内に流体が流入していると判定することができる。

【数32】

【数33】

ここで、添え字ｋは、ｋ番目の格子３０を意味する。Ｈ_ｋは、ｋ番目の格子３０の位置において、式（２７）によって得られる液面の高さを表す。ｕ_ｓｗｋは、ｋ番目の格子３０の位置において、式（３０）によって得られる速度ベクトルを表す。ｎ_ｋは、ｋ番目の格子３０の位置におけるバッファ領域Ω_ｎｓｂの端面の外向きの単位法線ベクトルを表す。

【0059】

式（３２）は、ｋ番目の格子３０の高さｙ_ｋが、液面の高さＨ_ｋより低いことを意味する。式（３３）は、ｋ番目の格子３０の位置における速度ベクトルが、バッファ領域Ω_ｎｓｂの内側を向いていること、すなわちバッファ領域Ω_ｎｓｂの端面を通してバッファ領域Ω_ｎｓｂ内に流体が流入していることを意味する。

【0060】

流体が流入していると判定された格子３０において、単位時間あたりに領域Ω_ｎｓへ流入する計算粒子２０の個数は、流速の高さ方向の分布を考慮して、以下の式で表すことができる。

【数34】

【数35】

ここで、Ｖ_ｋは、ｋ番目の格子３０の体積を表す。Ｎ_ｋは、格子３０ごとに定義される変数である。

【0061】

領域Ω_ｎｓへ流入させる計算粒子２０の位置、個数、及び時刻は、式（３５）で与えられる微分方程式を解くことにより算出することができる。具体的には、格子３０のそれぞれに定義された変数Ｎ_ｋの初期値を０とし、式（３５）にしたがって変数Ｎ_ｋを時間発展させる。変数Ｎ_ｋの値が１になった時刻において、ｋ番目の格子３０の位置に計算粒子２０を追加し、変数Ｎ_ｋの値を０にする。新たに追加した計算粒子２０に与える物理量は、相手側の領域Ω_ｓｗ内の計算粒子２０が持つ物理量に基づいて、物理量を交換する手法を用いて求めることができる。

【0062】

領域Ω_ｓｗから外に流出した計算粒子２０は、領域Ω_ｓｗにおける解析対象から除外する。この操作を繰り返すことで、流体の流れに応じて領域Ω_ｎｓに計算粒子２０が追加される。

【0063】

［ナビエストークス方程式の領域から浅水方程式の領域へ移動する計算粒子の制御］
次に、ナビエストークス方程式を適用する領域Ω_ｎｓから浅水方程式を適用する領域Ω_ｓｗへ流体が流れている場合の計算粒子２０の制御について説明する。

【0064】

ナビエストークス方程式の領域Ω_ｎｓから浅水方程式の領域Ω_ｓｗに流体が流入している場合も、上述の場合と同様に、計算粒子２０を追加する位置、個数、及び時刻を算出することができる。ただし、浅水方程式は二次元的であるため、バッファ領域Ω_ｓｗｂの端面に配置する格子を高さ方向に分布させる必要がない。また、液面の高さＨ_ｋは、式（２６）を用いて計算する。速度ベクトルｕ_ｋは、式（２９）を用いて計算する。

【0065】

［実施例によるシミュレーション装置及びシミュレーション方法］
次に、図７及び図８を参照して、一実施例によるシミュレーション装置及びシミュレーション方法について説明する。

【0066】

図７は、一実施例によるシミュレーション装置のブロック図である。本実施例によるシミュレーション装置は、入力部４０、処理部４１、出力部４２、及び記憶部４３を含む。入力部４０から処理部４１にシミュレーション条件等が入力される。さらに、ユーザから入力部４０に各種指令（コマンド）等が入力される。入力部４０は、例えば通信装置、リムーバブルメディア読取装置、キーボード、ポインティングデバイス等で構成される。

【0067】

処理部４１は、入力されたシミュレーション条件及び指令に基づいてシミュレーションを実行する。さらに、シミュレーション結果を出力部４２に出力する。シミュレーション結果には、例えば、シミュレーション対象物を構成する粒子の位置の時間変化等を表す情報が含まれる。処理部４１は、例えばコンピュータの中央処理ユニット（ＣＰＵ）を含む。本実施例によるシミュレーション方法の各手順をコンピュータに実行させるためのプログラムが、記憶部４３に記憶されている。出力部４２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

【0068】

図８は、本実施例によるシミュレーション方法の手順を示すフローチャートである。処理部４１（図７）が記憶部４３（図７）に格納されているプログラムを実行することにより、図８に示した各ステップの処理が実行される。

【0069】

処理部４１が、入力部４０に入力されたシミュレーション条件を取得する（ステップＳ１）。シミュレーション条件には、シミュレーション対象の流体を定義する情報、境界条件、初期条件等が含まれる。シミュレーション対象の流体を定義する情報には、流体の密度、粘度等の物性値が含まれる。境界条件には、解析空間Ω（図３）の形状及び大きさを指定する情報、解析空間Ωを複数の領域（例えば領域Ω_ｓｗ、Ω_ｎｓ（図３））に区分するための情報、各領域に適用される支配方程式を指定するための情報、区分された領域のそれぞれにおける計算粒子２０（図３）の直径を指定する情報、バッファ領域（例えばバッファ領域Ω_ｓｗｂ、Ω_ｎｓｂ（図３））を定義するための情報等が含まれる。初期条件には、流体の初期物性値を指定する方法、複数の計算粒子２０を解析空間Ωに配置するための情報等が含まれる。

【0070】

次に、入力されたシミュレーション条件に基づいて、解析空間Ωを複数の領域（例えば領域Ω_ｓｗ、Ω_ｎｓ（図３））に区分し、各領域内で用いる支配方程式を設定する（ステップＳ２）。次に、各領域に、バッファ領域（例えばバッファ領域Ω_ｓｗｂ、Ω_ｎｓｂ（図３））を定義する（ステップＳ３）。次に、区分された領域及びバッファ領域に、それぞれの領域に適用される支配方程式に応じた複数の計算粒子２０を配置する（ステップＳ４）。領域Ω_ｓｗ、Ω_ｎｓ内の計算粒子２０のそれぞれの物理量には初期値を与える。

【0071】

計算粒子２０を配置した後、バッファ領域内の計算粒子２０が持つ物理量の値を計算する（ステップＳ５）。例えば、バッファ領域Ω_ｓｗｂ、Ω_ｎｓｂにそれぞれ重なる相手側の領域Ω_ｎｓ、Ω_ｓｗの計算粒子２０に付与されている物理量の値に基づいて、バッファ領域Ω_ｓｗｂ、Ω_ｎｓｂの計算粒子２０に付与される物理量の値を計算する（物理量の交換を行う）。

【0072】

区分された複数の領域の各々について支配方程式を解くことにより、計算粒子２０に付与される物理量を計算する（ステップＳ６）。次に、区分された複数の領域の各々について、計算粒子２０の個数の制御を行う（ステップＳ７）。例えば、式（３２）～式（３５）を用いて、流体が流入する側の領域に新たに計算粒子２０を追加する。

【0073】

次に、区分された複数の領域、及びバッファ領域内の各々に含まれる計算粒子２０に付与されている物理量及び位置を更新する（ステップＳ８）。ステップＳ５からステップＳ８までの手順を、解析終了条件が満たされるまで繰り返し実行する（ステップＳ９）。解析終了条件が満たされると、処理部４１は出力部４２に解析結果を出力する（ステップＳ１０）。

【0074】

［優れた効果］
次に、上記実施例の優れた効果について説明する。
上記実施例では、解析空間の一部の領域について、二次元的な浅水方程式を用いて解析を行うため、解析空間の全域について三次元のナビエストークス方程式を用いて解析を行う場合と比べて、計算資源の削減が可能である。また、しぶきや液滴等が発生すると想定される領域については、ナビエストークス方程式を用いて解析を行うため、浅水方程式では解析できない現象を解析することが可能である。

【0075】

異なる支配方程式が適用される２つの領域の間で、計算粒子に付与される物理量の交換を行うことにより、異なる支配方程式が適用される２つの領域の間で物理量の連続性を確保することができる。一方の領域から他方の領域に流れる流体の速度ベクトルに応じて計算粒子の除去及び追加を行うことにより、一方の領域から他方の領域に流れる流体の流れを再現することができる。

【0076】

さらに、区分された複数の領域から相手側の領域に入り込むバッファ領域を定義して物理量の交換を行い、バッファ領域の端面に張られた格子を用いて計算粒子の流入を制御するため、解析空間を任意の形状に区分することが可能である。

【0077】

［実際のシミュレーション結果］
次に、図９～図１０Ｂを参照して、本実施例による方法で実際にシミュレーションを行った結果について説明する。

【0078】

図９は、シミュレーション対象の解析空間Ωを示す模式図である。解析空間Ωは直方体である。直方体の各面に直交する座標軸を持つｘｙｚ直交座標系を定義する。解析空間Ωのｘ方向の寸法Ｌｘ（長さ）を２０ｍとし、ｙ方向の寸法Ｌｙ（高さ）を２ｍとし、ｚ方向（図９の奥行方向）の寸法Ｌｚを１．６ｍとした。ｙ軸の負の向きに重力加速度ｇが作用する。初期状態において、解析空間Ωのｘ方向の中心からｘ軸の負の側、かつ高さＬｙｉが１ｍ以下の領域に水１５が配置されている。水が配置された領域のｘ方向の寸法Ｌｘｉは１０ｍである。水１５は重力によって駆動され、解析空間Ωのｘ軸の正の向きに流れる。

【0079】

解析空間Ωのｘ方向の中心からｘ軸の正の向きに長さＬｘ０だけ移動した位置を境界面とし、境界面よりｘ軸の負側の領域Ω_ｓｗに浅水方程式を適用し、正側の領域Ω_ｎｓにナビエストークス方程式を適用する。長さＬｘ０を１ｍとした。水の物性値には、標準大気圧下での温度２０℃における値を与えた。重力加速度ｇは９．８１ｍ／ｓ^２とした。壁面の境界条件には、すべり無し条件を与えた。計算粒子の直径は、２つの領域Ω_ｓｗ、Ω_ｎｓの両方において０．０５ｍに設定した。初期状態では、水の分布に合わせて、複数の計算粒子を、間隔が計算粒子の直径と等しくなるように、等間隔に配置した。解析の時間刻み幅は０．３７５×１０^－４秒とした。時刻０秒を初期状態として流れ場を時間発展させ、時刻１．５秒の時点のｚ方向の中央における水面の高さのｘ方向の分布を理論解析の結果と比較した。

【0080】

図１０Ａ及び図１０Ｂは、解析結果を示すグラフである。横軸はｘ方向の位置を単位［ｍ］で表し、縦軸は、水面の高さを単位［ｍ］で表す。初期状態のときに水１５が配置されていた領域を破線で表している。初期状態で水１５が配置されていた領域の側面の位置をｘ軸の原点とする。図１０Ａは、解析空間Ωの全体を示すグラフであり、図１０Ｂは、ｘ方向の位置が６．５ｍ以上９．５以下の範囲を拡大して示すグラフである。

【0081】

グラフ中の破線ｂは理論解を示し、実線ａは、本実施例によるシミュレーションによって得られた計算結果を示す。破線ｂの位置を中心として上方向及び下方向に延びる線分のそれぞれの長さは、計算粒子の直径に相当する。

【0082】

ｘ座標が１ｍの位置が、領域Ω_ｓｗとΩ_ｎｓとの境界面に相当し、境界面よりｘ軸の負側の領域が、浅水方程式を適用して解析した領域Ω_ｓｗに相当し、正側の領域が、ナビエストークス方程式を適用して解析した領域Ω_ｎｓに相当する。

【0083】

実線ａで示すシミュレーションによる水面の高さは、領域Ω_ｓｗとΩ_ｎｓとの境界で連続している。破線ｂで示す理論解と、実線ａで示すシミュレーション結果とを比較すると、両者の差は解析空間Ωの全域において計算粒子の直径以下であることがわかる。このシミュレーションにより、水面高さの分布が時間の経過とともに変化する流れ場を、浅水方程式とナビエストークス方程式とを組み合わせた本実施例によるシミュレーション方法を用いて効率的に解析可能であることが確認された。

【0084】

上記実施例では、ＳＰＨ法を用いて解析を行っているが、上記実施例の考え方は、その他の粒子法を用いた解析にも応用可能である。

【0085】

上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

【符号の説明】

【0086】

１０ 傾斜した壁面
１１ 液膜
１２ 液滴
１３ 液膜の表面
１５ 水
２０ 計算粒子
２０Ｓ 液膜の表面を構成する計算粒子
３０ 格子
４０ 入力部
４１ 処理部
４２ 記憶部
４３ 出力部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】