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

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

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

(11)【公開番号】P2024025512

(43)【公開日】2024-02-26

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

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20240216BHJP

【ＦＩ】

G16Z99/00

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

(21)【出願番号】P 2022129013

(22)【出願日】2022-08-12

(71)【出願人】

【識別番号】000002107

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

(74)【代理人】

【識別番号】100105887

【弁理士】

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

(72)【発明者】

【氏名】廣田 真人

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

5L049

【Ｆターム（参考）】

5L049DD02

(57)【要約】

【課題】粒子法を用いた解析において、解析精度を高く維持し、かつ計算量を削減することが可能なシミュレーション方法を提供する。
【解決手段】流体を複数の計算粒子で表し、複数の計算粒子のそれぞれに物理量に付与し、解析空間に複数の計算粒子を配置し、複数の計算粒子のそれぞれに付与された物理量及び複数の計算粒子の位置を、支配方程式を解いて時間発展させる。複数の計算粒子で表される流れ場の状態に応じて、複数の計算粒子の大きさに依存する空間解像度の過不足を、複数の計算粒子のそれぞれについて評価し、空間解像度の過不足の評価の結果に応じて、空間解像度を調整する。
【選択図】図２

【特許請求の範囲】

【請求項1】

流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量に付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、支配方程式を解いて時間発展させるシミュレーション方法であって、
前記複数の計算粒子で表される流れ場の状態に応じて、前記複数の計算粒子の大きさに依存する空間解像度の過不足を、前記複数の計算粒子のそれぞれについて評価し、
空間解像度の過不足の評価の結果に応じて、空間解像度を調整するシミュレーション方法。

【請求項2】

解析対象の流体は壁面に接触しており、
空間解像度の過不足の評価は、前記複数の計算粒子のそれぞれの位置において、流れ場のレイノルズ数、または前記壁面からの距離に基づいて行い、
空間解像度の調整は、前記複数の計算粒子を分割または合体させることにより行う請求項１に記載のシミュレーション方法。

【請求項3】

前記複数の計算粒子を分割または合体させることにより生じる前記複数の計算粒子の分布の歪みを低減させる方向に作用する人工移流項を前記支配方程式に付与し、
時間発展の際に、前記人工移流項が付与された前記支配方程式を解く請求項２に記載のシミュレーション方法。

【請求項4】

前記支配方程式はラグランジュ微分項を含んでおり、ラグランジュ微分項を、前記人工移流項を加味したものに修正して前記支配方程式を解くことにより、前記人工移流項を付与することによって生じる前記複数の計算粒子ごとの物理量の歪みを修正する方向に物理量を変化させる請求項３に記載のシミュレーション方法。

【請求項5】

粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて、粒子法を用いて複数の粒子のそれぞれの位置を時間発展させる処理部と
を備え、
前記処理部は、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量に付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、支配方程式を解いて時間発展させる手順と、
前記複数の計算粒子で表される流れ場の状態に応じて、前記複数の計算粒子の大きさに依存する空間解像度の過不足を、前記複数の計算粒子のそれぞれについて評価する手順と、
空間解像度の過不足の評価の結果に応じて、空間解像度を調整する手順と
を実行するシミュレーション装置。

【請求項6】

粒子法を用いて流体の流れを解析するシミュレーションをコンピュータに実行させるプログラムであって、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量に付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、支配方程式を解いて時間発展させる手順と、
前記複数の計算粒子で表される流れ場の状態に応じて、前記複数の計算粒子の大きさに依存する空間解像度の過不足を、前記複数の計算粒子のそれぞれについて評価する手順と、
空間解像度の過不足の評価の結果に応じて、空間解像度を調整する手順と
を実行させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

流体の運動を粒子系の運動として近似して、流体の挙動を解析するシミュレーション方法が知られている。このシミュレーション方法は、粒子法と呼ばれる。粒子法では、流体が複数の粒子で表され、流体の流れ場の中に配置される壁境界も一般的に複数の粒子で表される。

【先行技術文献】

【非特許文献】

【0003】

【非特許文献1】原田隆宏、越塚誠一、「ＳＰＨにおける壁境界計算手法の改良」、情報処理学会論文誌、Ｖｏｌ．４８、Ｎｏ．４、１８３８～１８４６頁、２００７年

【発明の概要】

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

【0004】

粒子法による流体解析においては、空間解像度が粒子の直径によって決まる。粒子の直径は、解析しようとする流体運動の長さスケールに対して十分小さく設定しなければならない。一般に、流体運動の長さスケールは場所によって異なる。例えば、流体内に配置される物体の壁から遠い領域では、粒子の直径を大きくしても十分精度の高い解析を行うことができる。ところが、壁の近傍では、流体運動の長さスケールは境界層の厚さで代表されるため、壁近傍の流体の運動を精度よく解析するためには、粒子の直径を境界層の厚さよりも十分小さくしなければならない。

【0005】

解析空間全体において高精度の解析を行うためには、壁から遠い領域においても、境界層の厚さに応じて粒子の直径を小さくしなければならない。このため、解析空間に配置する粒子数が多くなり、その結果、計算量が増大してしまう。

【0006】

本発明の目的は、粒子法を用いた解析において、解析精度を高く維持し、かつ計算量を削減することが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。

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

【0007】

本発明の一観点によると、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量に付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、支配方程式を解いて時間発展させるシミュレーション方法であって、
前記複数の計算粒子で表される流れ場の状態に応じて、前記複数の計算粒子の大きさに依存する空間解像度の過不足を、前記複数の計算粒子のそれぞれについて評価し、
空間解像度の過不足の評価の結果に応じて、空間解像度を調整するシミュレーション方法が提供される。

【0008】

本発明の他の観点によると、
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて、粒子法を用いて複数の粒子のそれぞれの位置を時間発展させる処理部と
を備え、
前記処理部は、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量に付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、支配方程式を解いて時間発展させる手順と、
前記複数の計算粒子で表される流れ場の状態に応じて、前記複数の計算粒子の大きさに依存する空間解像度の過不足を、前記複数の計算粒子のそれぞれについて評価する手順と、
空間解像度の過不足の評価の結果に応じて、空間解像度を調整する手順と
を実行するシミュレーション装置が提供される。

【0009】

本発明のさらに他の観点によると、
粒子法を用いて流体の流れを解析するシミュレーションをコンピュータに実行させるプログラムであって、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量に付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、支配方程式を解いて時間発展させる手順と、
前記複数の計算粒子で表される流れ場の状態に応じて、前記複数の計算粒子の大きさに依存する空間解像度の過不足を、前記複数の計算粒子のそれぞれについて評価する手順と、
空間解像度の過不足の評価の結果に応じて、空間解像度を調整する手順と
を実行させるプログラムが提供される。

【発明の効果】

【0010】

空間解像度の過不足の評価の結果に応じて、空間解像度を調整することにより、解析空間全体において、適切な空間解像度で解析を行うことができる。

【図面の簡単な説明】

【0011】

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

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

【図3】図３は、解析空間、及び解析空間内に配置された複数の計算粒子を示す模式図である。

【図4】図４は、空間解像度を評価する方法の手順を示すフローチャートである。

【図5】図５は、壁面及び評価対象の計算粒子の位置関係を表す模式図である。

【図6】図６は、ステップＳＡ６（図２）で実行する空間解像度の調整方法の手順を示すフローチャートである。

【図7】図７Ａは、分割前の計算粒子の模式図であり、図７Ｂは、１つの計算粒子を分割して生成された２個の計算粒子の模式図である。

【図8】図８Ａは、１つの計算粒子を分割して生成された４個の計算粒子の模式図であり、図８Ｂは、１つの計算粒子を分割して生成された８個の計算粒子の模式図である。

【図9】図９Ａは、複数の計算粒子の分布の一例を示す模式図であり、図９Ｂは、２つの計算粒子を合体させた後の計算粒子の分布の一例を示す模式図である。

【図10】図１０Ａは、計算粒子を分割する前の複数の計算粒子の分布を示す模式図であり、図１０Ｂは、１つの計算粒子を分割した後の複数の計算粒子の分布を示す模式図であり、図１０Ｃは、図１０Ｂに示した状態から一定のタイムステップが経過した時点の複数の計算粒子の分布を模式図である。

【図11】図１１Ａは、２つの計算粒子を合体させる前の複数の計算粒子の分布を示す模式図であり、図１１Ｂは、２つの計算粒子を合体させた後の複数の計算粒子の分布を示す模式図であり、図１１Ｃは、図１１Ｂに示した状態から一定のタイムステップが経過した時点の複数の計算粒子の分布を模式図である。

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

【図13】図１３は、実際にシミュレーションを行った解析モデルの模式図である。

【図14】図１４Ａ及び図１４Ｂは、それぞれ計算開始時点及び計算終了時点における計算粒子の分布を示す模式図である。

【図15】図１５は、計算終了時点の速度分布及び初期状態の速度分布を示すグラフである。

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

【0012】

本発明の一実施例について説明する前に、一般的なＳＰＨ（Smoothed Particle Hydrodynamics）法について簡単に説明する。ＳＰＨ法は、連続体を複数の計算粒子の集合を用いて離散化し、流れ場の時間発展を計算するラグランジュ型の数値計算手法である。ＳＰＨ法では、粒子１個の質量分布として、連続かつ微分可能なカーネル関数を用いるとともに、カーネル関数の微分係数の重ね合わせを空間微分とする離散化モデルを用いる。

【0013】

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

【数1】

本明細書において、下付きの添え字ｓ及びｍは、計算粒子のそれぞれに付与された通し番号を表し、下付きの添え字ｓ、ｍが付された物理量は、それぞれｓ番目及びｍ番目の計算粒子に付された物理量を意味する。

【0014】

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

【数2】

【0015】

ｈは、平滑化長さであり、カーネル幅ともいわれる。カーネル幅ｈは、複数の粒子を同一の質量及び物性密度の球体に置き換えた等価球の直径に等しい。本明細書において、カーネル幅ｈを粒子の直径という場合がある。

【0016】

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

【数3】

【0017】

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

【数4】

【数5】

【0018】

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

【数6】

なお、式（２）の関数Ｗ（ｑ）は１変数の関数であるため、その勾配∇Ｗはスカラー量である。

【0019】

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

【数7】

【0020】

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

【数8】

【0021】

［流れの支配方程式］
次に、圧縮性流れを解析対象とする場合の支配方程式について説明する。高レイノルズ数の圧縮性流れを、限られた計算資源で解析する手段として、例えばReynolds-Averaged Navier-Stokes（ＲＡＮＳ）解析を用いることができる。これは、流れの支配方程式をレイノルズ分解によって平均場と変動場とに分解し、平均場に対する変動場の寄与をモデル化し、平均場の支配方程式だけを数値的に解くことで、少ない計算資源で高レイノルズ数流れの平均場を解析する技術である。

【0022】

ＲＡＮＳ解析の定式化方法として、例えばShear Stress Transport（ＳＳＴ）モデルを用いることができる。ＳＳＴモデルについては、F. R. MENTER, “Two-equation eddy-viscosity turbulence models for engineering applications”, AIAA Journal Vol.32 No.8, (1994), pp.1598-1605, _eprint: https://doi.org/10.2514/3.12149に説明されている。

【0023】

ＳＳＴモデルを用いた圧縮性流れの支配方程式は、例えば次のように記述することができる。

【数9】

【数10】

【数11】

【数12】

【数13】

ここで、ｔは時間、ｘは空間座標、ρは密度、ｖは速度、ｐは圧力、ｕは内部エネルギ、及びＴは温度を表す。下付きの添え字ｉ、ｊ、ｋは、ベクトル量またはテンソル量の成分を表す。ｋ及びωは、ＳＳＴモデルの変数であり、ｋは乱流の運動エネルギ、ωは乱流の運動エネルギの散逸時間スケールに相当する量である。ｃ_ｐは低圧比熱であり、Ｐｒはプラントル数を表す。

【0024】

σ_ｉｊは粘性応力テンソルを表し、次の式で与えられる。

【数14】

ここで、Ｓ_ｉｊは歪み速度テンソルを表す。μは粘性係数であり、例えばサザーランドの式を用いて算出することができる。δ_ｉｊはクロネッカーのδを表す。τ_ｉｊはレイノルズ応力項であり、次式による表現が多く用いられる。

【数15】

【0025】

ここで、μ_ｔは渦粘性係数である。圧力ｐは、例えば下記の理想気体の状態方程式を用いて算出することができる。

【数16】

ここで、Ｒは気体状数である。

【0026】

式（９）～式（１３）のその他の変数は、ＳＳＴモデルで用いられる定数である。ＳＳＴモデルは一般的に広く用いられていることから、これらの定数については、ここでは説明を省略する。

【0027】

式（９）～式（１３）を、式（４）を用いて離散化して数値的に解くことにより、ＳＰＨ法を用いて、高レイノルズ数の圧縮性流れの解析を行うことができる。

【0028】

ＳＰＨ法によって離散化された支配方程式は、以下のように記述される。

【数17】

【数18】

【数19】

【数20】

【数21】

【数22】

ここで、δｖ_ｉは計算粒子の空間分布の不均一性を均すための人工移流項である。

【0029】

人工移流項δｖ_ｉとして、例えば以下の式を用いることができる。

【数23】

【数24】

ここで、Ｓ及びκは、人工移流項の強さを示すパラメータである。例えば、Ｓ＝０．４、κ＝４と設定することができる。ｃ_０は代表音速である。代表音速ｃ_０として、よどみ点条件における音速を用いることができる。

【0030】

［空間解像度の過不足］
次に、粒子法を用いて流れ場を解析する場合の空間解像度の過不足について説明する。ＳＰＨ法をはじめとする粒子法を用いた流体解析においては、解析の空間解像度は計算粒子の直径（例えば、カーネル幅）によって決まる。一般に、流体運動の長さスケールは場所によって異なる。例えば、物体の壁面近傍では流体運動の長さスケールは、境界層厚さで代表され、相対的に高い空間解像度が必要とされる。

【0031】

一方で、壁面からの距離が遠い領域では、流体運動の長さスケールは、境界層厚さより大きいことが多い。このため、壁面からの距離が遠い領域では、壁面近傍の領域と比べて、低い空間解像度で解析を行うことができる。

【0032】

高い空間解像度で解析を行う必要がある領域に、直径の大きな計算粒子を配置すると、解析の空間解像度が不足してしまう。逆に、低い空間解像度で解析を行うことが可能な領域に、直径の小さな計算粒子を配置すると、解析の空間解像度が過剰になってしまう。以下に説明する実施例では、解析の空間解像度の過不足を低減させ、必要十分な空間解像度で解析を行うことができる。

【0033】

［実施例によるシミュレーション手順］
図１及び図２を参照して、本実施例によるシミュレーション装置及びシミュレーション方法について説明する。

【0034】

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

【0035】

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

【0036】

図２は、本実施例によるシミュレーション方法の手順を示すフローチャートである。
まず、入力部３０に入力されたシミュレーション条件を処理部３１が取得する（ステップＳＡ１）。シミュレーション条件には、例えば、シミュレーション対象の流体を定義する情報、境界条件、初期条件等が含まれる。流体を定義する情報には、流体の密度、粘度等の物性値が含まれる。境界条件には、解析空間の形状、大きさを指定する情報が含まれる。初期条件には、流体の初期物性値を指定する方法、計算粒子のカーネル幅の初期値、及び複数の計算粒子を解析空間に配置するための情報が含まれる。

【0037】

シミュレーション条件には、さらに、空間解像度の過不足の評価指標、計算粒子の分割方法、最大カーネル幅を指定するための情報が含まれる。これらの情報は、以下に説明する本実施例のシミュレーションにおいて参照される。

【0038】

次に、処理部３１は、取得したシミュレーション条件に基づいて解析空間に複数の計算粒子を配置する（ステップＳＡ２）。

【0039】

図３は、解析空間１０、及び解析空間１０内に配置された複数の計算粒子２０を示す模式図である。一例として、解析空間１０は直方体の形状を有しており、解析空間１０の各面に、滑り無し条件、周期境界条件等の条件が課される。

【0040】

次に、複数の計算粒子２０のそれぞれについて式（１７）～式（２２）に示した支配方程式を数値的に解くことにより、計算粒子２０のそれぞれに付与される物理量を計算し、計算粒子２０のそれぞれの物理量及び位置を更新する（ステップＳＡ３）。

【0041】

次に、空間解像度の過不足を評価するか否かを判定する（ステップＳＡ４）。例えば、処理部３１は、空間解像度の過不足の評価を行った直近のタイムステップから一定のタイムステップ数が経過した時点で、空間解像度の過不足を評価すると判定する。

【0042】

処理部３１は、空間解像度の過不足を評価すると判定した場合は、空間解像度の過不足の評価を行う（ステップＳＡ５）。例えば、解析に必要とされる空間解像度に比べて計算粒子２０の直径が大きすぎる場合は、空間解像度が不足していると判定する。逆に、解析に必要とされる空間解像度に比べて計算粒子２０の直径が小さすぎる場合は、空間解像度が過剰であると判定する。

【0043】

処理部３１は、空間解像度の過不足を評価した後、評価結果に基づいて空間解像度を調整する（ステップＳＡ６）。例えば、空間解像度が不足していると判定された領域においては、計算粒子２０の直径を現在の直径より小さくすることにより、空間解像度を高める。このとき、計算粒子２０の直径を小さくしたことに伴い、計算粒子２０の個数を増加させる。空間解像度が過剰であると判定された領域においては、計算粒子２０の直径を現在の直径より大きくすることにより、空間解像度を低下させる。このとき、計算粒子２０の直径を大きくしたことに伴い、計算粒子２０の個数を減少させる。

【0044】

ステップＳＡ４で、空間解像度の過不足を評価しないと判定した場合は、ステップＳＡ５及びステップＳＡ６を実行しない。

【0045】

解析終了条件が満たされるまで、ステップＳＡ３以降の処理を繰り返し実行する（ステップＳＡ７）。解析終了条件は、例えば実行するタイムステップ数等で規定される。解析を終了したら、解析結果を出力部３２に出力する（ステップＳＡ８）。

【0046】

［空間解像度の過不足の評価］
次に、図４及び図５を参照して、ステップＳＡ５（図２）で実行する空間解像度の評価方法の一例について説明する。

【0047】

図４は、空間解像度を評価する方法の手順を示すフローチャートである。この評価は、複数の計算粒子２０のそれぞれについて実行する。

【0048】

まず、評価対象の計算粒子２０が壁面に近接しているか否かを判定する（ステップＳＢ１）。ここで、図３に示した解析空間１０の壁面のうち、周期境界条件を適用する壁面については、判定対象の壁面に含めない。例えば、滑り無し条件を適用している壁面について判定を行う。

【0049】

図５は、壁面１１及び評価対象の計算粒子２０ｓの位置関係を表す模式図である。壁面１１の表面をｘｚ面とし、壁面１１の法線方向をｙ方向とするｘｙｚ直交座標系を定義する。計算粒子２０ｓの直径をｈ_ｓと標記し、壁面１１から計算粒子２０ｓの中心までの距離をｙ_ｓと標記する。一例として、ｙ_ｓ＜ｈ_ｓを満たす計算粒子２０ｓを、壁面１１の近傍に位置する粒子と判定する。流体の速度ベクトルの壁面１１に平行な成分をｖ_ｈと標記する。壁面１１の近傍においては、壁面１１から遠ざかるにしたがって速度ｖ_ｈが大きくなる。

【0050】

壁面に近接していない計算粒子２０ｓについて、計算粒子２０ｓの中心位置におけるせん断強さγドットによって定義されるレイノルズ数Ｒｅ_γｓを、以下の式を用いて計算する（ステップＳＢ２）。

【数25】

ここで、せん断強さγドットは、以下の式を用いて計算することができる。

【数26】

ここで、ｖ_ｘ、ｖ_ｙ、ｖ_ｚは、それぞれ計算粒子２０の速度ベクトルのｘ成分、ｙ成分、及びｚ成分を表す。

【0051】

処理部３１は、計算粒子２０ｓの位置におけるレイノルズ数Ｒｅ_γｓを計算した後、レイノルズ数Ｒｅ_γｓが１０より大きいか否かを判定する（ステップＳＢ４）。レイノルズ数Ｒｅ_γｓが１０より大きい場合、計算粒子２０ｓの位置における空間解像度は不足であると判定する（ステップＳＢ６）。レイノルズ数Ｒｅ_γｓが１０以下のとき、処理部３１は、レイノルズ数Ｒｅ_γｓが２．５より小さいか否かを判定する（ステップＳＢ８）。レイノルズ数Ｒｅ_γｓが２．５より小さい場合、計算粒子２０ｓの位置における空間解像度は過剰であると判定する（ステップＳＢ１０）。本明細書において、計算粒子の位置における空間解像度を、単に「計算粒子の空間解像度」という場合がある。

【0052】

計算粒子２０ｓが壁面１１に近接している場合、処理部３１は壁面１１から計算粒子２０ｓの中心までの無次元距離ｙ_ｓ ^＋を、以下の式により計算する（ステップＳＢ３）。

【数27】

ここで、ｖ_τｓは壁面摩擦速度であり、以下の式で表される。

【数28】

【数29】

ここで、τ_ｗは、壁面１１に作用する摩擦応力を表す。なお、壁面１１の近傍に位置する計算粒子２０ｓに対しては、近似的に、計算粒子２０ｓの位置における速度勾配を用いて壁面摩擦応力τ_ｗｓを以下の式で計算することができる。

【数30】

【0053】

計算粒子２０ｓについて無次元距離ｙ_ｓ ^＋が求まると、処理部３１は、無次元距離ｙ_ｓ ^＋が３より大きいか否かを判定する（ステップＳＢ５）。無次元距離ｙ_ｓ ^＋が３より大きい場合、計算粒子２０ｓの位置の空間解像度は不足であると判定する（ステップＳＢ７）。無次元距離ｙ_ｓ ^＋が３以下の場合、処理部３１は、無次元距離ｙ_ｓ ^＋が１より小さいか否かを判定する（ステップＳＢ９）。無次元距離ｙ_ｓ ^＋が１より小さい場合、処理部３１は、計算粒子２０ｓの空間解像度は過剰であると判定する（ステップＳＢ１１）。

【0054】

すべての計算粒子２０について、空間解像度の過不足を判定する上記ステップを繰り返す（ステップＳＢ１２）。

【0055】

［空間解像度の調整方法］
次に、図６～図９Ｂを参照して、ステップＳＡ６（図２）で実行する空間解像度の調整方法の一例について説明する。

【0056】

図６は、ステップＳＡ６（図２）で実行する空間解像度の調整方法の手順を示すフローチャートである。この調整は、複数の計算粒子２０のそれぞれについて実行する。

【0057】

処理部３１は、ステップＳＡ５で行った評価結果に基づいて、着目する計算粒子２０の空間解像度が過剰か否かを判定する（ステップＳＣ１）。空間解像度が過剰ではないと判定された場合、計算粒子２０の空間解像度が不足しているか否かを判定する（ステップＳＣ５）。計算粒子２０の空間解像度が不足している場合、計算粒子２０を分割し（ステップＳＣ６）、カーネル幅がより小さい複数の計算粒子を生成する。

【0058】

１つの計算粒子２０を複数の計算粒子２０に分割する場合には、計算粒子２０の質量や体積等の物理量が保存される必要がある。また、分割後の複数の計算粒子２０の個数及び配置には任意性がある。

【0059】

図７Ａは、分割前の計算粒子２０ａの模式図である。図７Ｂは、１つの計算粒子２０ａを分割して生成された２個の計算粒子２０ａ１、２０ａ２の模式図である。図８Ａは、１つの計算粒子２０ａを分割して生成された４個の計算粒子２０ａ１～２０ａ４の模式図である。図８Ｂは、１つの計算粒子２０ａを分割して生成された８個の計算粒子２０ａ１～２０ａ８の模式図である。

【0060】

まず、図７Ｂに示すように、２個の計算粒子２０ａ１、２０ａ２に分割する場合について説明する。分割前の計算粒子２０ａの位置をｘ_１、カーネル幅をｈ_１、質量をｍ_１、体積をＶ_１と標記する。分割後の２つの計算粒子２０ａ１、２０ａ２の各々のカーネル幅をｈ_２、質量をｍ_２、体積をＶ_２と標記する。分割後の２つの計算粒子２０ａ１、２０ａ２の位置を、それぞれｘ_２１、ｘ_２２と標記する。ここで、Ｖ_１＝ｈ_１ ^３、Ｖ_２＝ｈ_２ ^３である。

【0061】

分割後の計算粒子２０ａ１、２０ａ２のそれぞれの質量ｍ_２、体積Ｖ_２、カーネル幅ｈ_２は、以下の式で計算することができる。

【数31】

【0062】

分割後の２つの計算粒子２０ａ１、２０ａ２の位置ｘ_２１、ｘ_２２は、例えば以下の式で与えることができる。

【数32】

ここで、ベクトルｉ、ｊ、ｋは、それぞれ解析空間１０（図３）に定義されたｘｙｚ直交座標系のｘ方向、ｙ方向、ｚ方向の単位ベクトルを表す。

【0063】

次に、図８Ａに示すように、４個の計算粒子２０ａ１～２０ａ４に分割する場合について説明する。分割後の４個の計算粒子２０ａ１～２０ａ４の各々のカーネル幅をｈ_２、質量をｍ_２、体積をＶ_２と標記する。分割後の４個の計算粒子２０ａ１～２０ａ４の位置を、それぞれｘ_２１～ｘ_２４と標記する。

【0064】

分割後の計算粒子２０ａ１～２０ａ４のそれぞれの質量ｍ_２、体積Ｖ_２、カーネル幅ｈ_２は、以下の式で計算することができる。

【数33】

【0065】

分割後の４個の計算粒子２０ａ１～２０ａ４の位置ｘ_２１～ｘ_２４は、例えば以下の式で与えることができる。

【数34】

【0066】

次に、図８Ｂに示すように、８個の計算粒子２０ａ１～２０ａ８に分割する場合について説明する。分割後の８個の計算粒子２０ａ１～２０ａ８の各々のカーネル幅をｈ_２、質量をｍ_２、体積をＶ_２と標記する。分割後の８個の計算粒子２０ａ１～２０ａ８の位置を、それぞれｘ_２１～ｘ_２８と標記する。

【0067】

分割後の計算粒子２０ａ１～２０ａ８のそれぞれの質量ｍ_２、体積Ｖ_２、カーネル幅ｈ_２は、以下の式で計算することができる。

【数35】

【0068】

分割後の８個の計算粒子２０ａ１～２０ａ８の位置ｘ_２１～ｘ_２８は、例えば以下の式で与えることができる。

【数36】

【0069】

分割後のｓ番目の計算粒子２０ｓの任意の物理量Ａ_２ｓは、例えば分割前の計算粒子２０の物理量Ａ_１の勾配を用いて以下の式で与えることができる。

【数37】

ここで、物理量Ａ_２ｓ、Ａ_１は、式（１８）～式（２２）の支配方程式の左辺の物理量に相当する。

【0070】

なお、解析空間が二次元である場合も、同様の手法で分割後の計算粒子２０の位置、及び計算粒子２０の物理量を求めることができる。

【0071】

ステップＳＣ１（図６）で空間解像度が過剰であると判定された場合、最近接ペアを構成する相手粒子があるか否かを判定する（ステップＳＣ２）。「最近接ペア」とは、互いに最近接の関係を有する計算粒子２０のペアを意味する。

【0072】

図９Ａを参照して、最近接ペアについて具体的に説明する。図９Ａは、複数の計算粒子２０の分布の一例を示す模式図である。図９Ａに示した矢印は、計算粒子２０のそれぞれから最近接の計算粒子２０に向かう。例えば、計算粒子２０ａの最近接粒子は、計算粒子２０ｂであり、計算粒子２０ｂの最近接粒子は計算粒子２０ａである。したがって、２つの計算粒子２０ａ、２０ｂは、最近接ペアを構成する。

【0073】

計算粒子２０ｃの最近接粒子は計算粒子２０ｄである。ところが、計算粒子２０ｄの最近接粒子は、計算粒子２０ｅであり、計算粒子２０ｃではない。したがって、計算粒子２０ｃは、他の計算粒子２０と最近接ペアを構成しない。計算粒子２０ｅの最近接粒子は計算粒子２０ｄである。このため、計算粒子２０ｄと２０ｅとが、最近接ペアを構成している。

【0074】

ステップＳＣ２（図６）において最近接ペアを構成する相手の計算粒子２０が有ると判定された場合、最近接ペアを構成する相手粒子の空間解像度が過剰か否かを判定する（ステップＳＣ３）。相手粒子の空間解像度が過剰であると判定された場合は、最近接ペアを構成する２つの計算粒子２０を合体させ（ステップＳＣ４）、２つの計算粒子２０のいずれのカーネル幅より大きいカーネル幅を持つ１つの大きな計算粒子を生成する。すなわち、最近接ペアを構成する２つの計算粒子２０のそれぞれの空間解像度が共に過剰である場合に、２つの計算粒子２０を合体させる。

【0075】

図９Ａに示した計算粒子２０のうち、計算粒子２０ａ、２０ｂ、２０ｃ、２０ｄは、空間解像度が過剰と判定されたものであり、他の計算粒子２０は、空間解像度が過剰と判定されていない。図９Ａにおいて、空間解像度が過剰と判定された計算粒子２０ａ、２０ｂ、２０ｃ、２０ｄにハッチングを付している。

【0076】

最近接ペアを構成する計算粒子２０ａと２０ｂとは、共に空間解像度が過剰と判定されたものであるため、計算粒子２０ａと２０ｂとを合体させる。最近接ペアを構成する計算粒子２０ｄと２０ｅとは、一方の計算粒子２０ｅの空間解像度が過剰ではないため、計算粒子２０ｄと２０ｅとは合体させない。

【0077】

図９Ｂは、計算粒子２０ａと２０ｂとを合体させた後の計算粒子２０の分布の一例を示す模式図である。計算粒子２０ａと２０ｂとが合体されて、大きな計算粒子２０ａｂが生成されている。図９Ｂにおいて、合体後の計算粒子２０ａｂに、相対的に淡いハッチングを付している。

【0078】

次に、２つの計算粒子２０を合体させる方法について説明する。
合体前の一方の計算粒子２０ａの位置をｘ_１１、カーネル幅をｈ_１１、質量をｍ_１１、体積をＶ_１１と標記し、他方の計算粒子２０ｂの位置をｘ_１２、カーネル幅をｈ_１２、質量をｍ_１２、体積をＶ_１２と標記する。合体後の計算粒子２０ａｂ（図９Ｂ）の質量ｍ_２、カーネル幅ｈ_２は、以下の式を用いて計算することができる。

【数38】

【0079】

合体後の計算粒子２０ａｂの位置ｘ_２は、合体前の２つの計算粒子２０ａ、２０ｂの位置ｘ_１１、ｘ_１２を、質量に基づく重み付き平均値を用いて、以下の式で決めることができる。

【数39】

【0080】

合体後の計算粒子２０ａｂの任意の物理量Ａ_２も同様に、合体前の計算粒子２０ａ、２０ｂの物理量Ａ_１１、Ａ_１２を用いて以下の式で計算することができる。

【数40】

または、合体前の計算粒子２０ａ、２０ｂの物理量Ａ_１１、Ａ_１２と、これらの物理量の勾配とを用いて、例えば以下の式で計算してもよい。

【数41】

【0081】

すべての計算粒子２０について、空間解像度の調整の処理が終了するまで、ステップＳＣ１からの手順を繰り返す（ステップＳＣ７）。

【0082】

次に、本実施例の優れた効果について説明する。
本実施例では、図２に示したステップＳＡ５で空間解像度の過不足を評価し、ステップＳＡ６で空間解像度を調整するため、解析空間内の全域で、適切な空間解像度で解析を行うことができる。例えば、壁面の近傍では高い空間解像度で解析が行われることになり、流れ場に関する十分な情報を得ることができる。また、壁面から遠い領域では、相対的に低い空間解像度で解析が行われることになり、計算負荷を軽減することができる。

【0083】

解析空間内の場所によって異なる空間解像度で解析を行う手法として、解析空間を、空間解像度の異なる複数の領域に分割し、分割された領域ごとに、大きさの異なる計算粒子を配置する手法が考えられる。ただし、この手法では、分割された領域の間で物理量を交換する処理を行う必要があり、この処理の計算負荷が大きい。分割数を増やすと、計算負荷がさらに大きくなってしまう。また、複雑な流れ場においては、分割された領域を跨って計算粒子が移動する際に、質量保存則に不整合が生じ、解析が数値的に不安定になる場合がある。さらに、解析空間の分割はオペレータの手作業で行われるため、オペレータの技量によっては、不適切な空間解像度の設定がなされてしまう可能性もある。

【0084】

これに対して本実施例では、解析空間を複数の領域に分割しないため、上述のような課題が生じない。

【0085】

本実施例では、図４に示したステップＳＢ４及びステップＳＢ８での判定の基準となる閾値を、１０及び２．５に設定し、ステップＳＢ５及びステップＳＢ９での判定の基準となる閾値を、３及び１に設定しているが、これらの閾値として他の値を用いてもよい。ただし、空間解像度が過剰と判定する基準となる閾値を、空間解像度が不足であると判定する基準となる閾値より小さくし、両者の間に差を設けることが好ましい。このように閾値を設定することにより、計算粒子２０が分割と合体とを頻繁に繰り返すことが防止される。

【0086】

また、少ないタイムステップ数の間に、計算粒子２０が分割や合体を繰り返すと、解析空間１０（図３）を満たす複数の計算粒子２０の分布に不整合が生じ、計算が数値的に不安定になる場合がある。これを回避するために、空間解像度を調整する処理（ステップＳＡ６）を、ある程度長い時間間隔で実行することが好ましい。この時間間隔に相当するタイムステップ数をＮ_ｒと標記し、最大カーネル幅をｈ_ｍａｘと標記し、流体の代表音速をｃ_０と標記し、計算の時間刻み幅をΔｔと標記したとき、タイムステップ数Ｎ_ｒを以下のように設定するとよい。

【数42】

ここで、最大カーネル幅ｈ_ｍａｘ及びαは任意のパラメータであり、最大カーネル幅ｈ_ｍａｘは、例えば流れ場の代表長さの１／１０以上１／３０以下程度の値に設定することができ、αは例えば１０に設定することができる。

【0087】

［計算粒子の分割、合体に伴う誤差の低減］
次に、図１０Ａ～図１１Ｃを参照して、計算粒子２０を分割、合体させたことにより生じ得る誤差を低減させる方法について説明する。

【0088】

図１０Ａは、計算粒子２０ａを分割する前の複数の計算粒子２０の分布を示す模式図である。図１０Ａにおいて、分割すべき計算粒子２０ａにハッチングを付している。分割前の状態では、計算粒子２０ａを含めて複数の計算粒子２０がほぼ均一に分布している。

【0089】

図１０Ｂは、計算粒子２０ａ（図１０Ａ）を分割した後の複数の計算粒子２０の分布を示す模式図である。１個の計算粒子２０ａ（図１０Ａ）が２個の計算粒子２０ａ１と２０ａ２とに分割されている。図１０Ｂにおいて、分割された２個の計算粒子２０ａ１、２０ａ２にハッチングを付している。

【0090】

１つの計算粒子２０ａ（図１０Ａ）を２つの計算粒子２０ａ１、２０ａ２に分割したことにより、２つの計算粒子２０ａ１、２０ａ２を含む複数の計算粒子２０の分布に偏りが生じる。

【0091】

図１１Ａは、２つの計算粒子２０ａ及び２０ｂを合体させる前の複数の計算粒子２０の分布を示す模式図である。図１１Ａにおいて、合体前の２つの計算粒子２０ａ、２０ｂにハッチングを付している。合体前の状態では、計算粒子２０ａ、２０ｂを含めて複数の計算粒子２０がほぼ均一に分布している。

【0092】

図１１Ｂは、計算粒子２０ａ及び２０ｂ（図１１Ａ）を合体させた後の複数の計算粒子２０の分布を示す模式図である。２つの計算粒子２０ａ、２０ｂ（図１１Ａ）が合体されて１つの計算粒子２０ａｂが生成されている。図１１Ｂにおいて、合体により生成された計算粒子２０ａｂにハッチングを付している。

【0093】

２つの計算粒子２０ａ、２０ｂ（図１１Ａ）を合体させて１つの計算粒子２０ａｂを生成したことにより、合体後の計算粒子２０ａｂを含む複数の計算粒子２０の分布に偏りが生じる。

【0094】

図１０Ｂ及び図１１Ｂに示したように、計算粒子２０の分布に偏りが生じると、式（１７）に付加している人工移流項δｖ_ｉの作用により、計算粒子２０のそれぞれに、分布の偏りを解消する方向に人工的な速度が付加される。図１０Ｂ及び図１１Ｂに示された矢印は、人工的に付加される速度の方向を示す。

【0095】

図１０Ｃは、図１０Ｂに示した状態から一定のタイムステップが経過した時点の複数の計算粒子２０の分布を示す模式図である。図１１Ｃは、図１１Ｂに示した状態から一定のタイムステップが経過した時点の複数の計算粒子２０の分布を示す模式図である。式（１７）の人工移流項δｖ_ｉの作用によって、複数の計算粒子２０の分布の偏りが解消されている。

【0096】

この人工移流項によってもたらされる計算粒子２０の運動は、流体運動とは無関係のものである。このため、計算粒子２０の分布の偏りが解消されることによって、流れ場が元々持っていた物理量の分布に歪みが生じてしまう。以下、この歪みを抑制する方法について説明する。

【0097】

式（１７）～式（２２）に示した支配方程式はラグランジュ微分項を含んでいる。以下に説明する方法では、このラグランジュ微分項を、人工移流項を加味したものに修正することにより、物理量の歪みを低減させる。

【0098】

式（１７）の人工移流項δｖ_ｉによる非物理的な運動に伴う物理量の歪みを抑制するために、人工移流項を加味したラグランジュ微分ｄＡ／ｄｔを以下の式で定義する。なお、ＤＡ／Ｄｔは、人工移流項を加味していない通常のラグランジュ微分である。

【数43】

式（４３）は、観測者が式（１７）の人工移流項δｖ_ｉの方向に移動したときに観測される物理量Ａの変化を抑制する方向に、物理量Ａを変化させることを意味する。

【0099】

式（４３）の右辺第２項は、任意のスカラー量Ａに対して以下のように書き換えることができる。

【数44】

【0100】

式（４３）の右辺第２項は、任意のベクトル量Ａに対して以下のように書き換えることができる。

【数45】

【0101】

式（４３）を用いて、式（１７）～式（２２）を書き換えると、ＳＰＨ法で解くべき支配方程式は、以下のように記述される。

【数46】

【数47】

【数48】

【数49】

【数50】

【数51】

【0102】

式（４７）～式（５１）は、式（１８）～式（２２）に含まれるラグランジュ微分項を、人工移流項を加味したものに修正した支配方程式である。式（４６）～式（５１）に示した支配方程式を解くことにより、人工移流項を付与することによって生じる物理量の歪みが修正される方向に物理量が変化する。これにより、計算粒子２０の分割及び合体と人工移流項の作用とによって発生する物理量の歪みを、低減させることができる。

【0103】

［空間解像度の過不足を評価する時期］
次に、図１２を参照して、空間解像度の過不足を評価する時期を好適化した他の実施例によるシミュレーション方法について説明する。図１２は、本実施例によるシミュレーション方法の手順を示すフローチャートである。以下、図２に示したフローチャートとの差に着目して説明する。

【0104】

粒子法を用いた流体解析手法においては、計算粒子２０のそれぞれについて、その計算粒子２０の近傍に存在する他の計算粒子２０の番号が保存された近接粒子リストを用いた解析が広く用いられている。この近接粒子リストを更新する処理は、一定の間隔で実行される。図１２に示すように、タイムステップごとに、近接粒子リストを更新するか否かを判定する（ステップＳＡ１０）。近接粒子リストを更新するタイムステップにおいて、近接粒子リストを更新し（ステップＳＡ１１）、その他のタイムステップでは、近接粒子リストの更新は行わない。

【0105】

図２に示した方法では、タイムステップごとに、空間解像度の過不足を評価するか否かを判定する（ステップＳＡ４）。これに対して図１２に示した実施例では、近接粒子リストを更新するタイムステップのときにのみ、近接粒子リストを更新する前に、空間解像度の過不足を評価するか否かを判定する（ステップＳＡ４）。

【0106】

近接粒子リストを更新する間隔は、一般的に、式（４２）で計算されるタイムステップ数Ｎ_ｒより短い。ステップＳＡ４において、直近に空間解像度の過不足の評価を行った時点からタイムステップ数Ｎ_ｒ以上経過した場合に、空間解像度の過不足の評価を行うと判定すればよい。

【0107】

計算粒子２０の分割や合体を行った後は、近接粒子リストを更新する必要がある。したがって、近接粒子リストを更新するタイミングで計算粒子２０の分割や合体を行うことにより、計算効率を高めることができる。

【0108】

［実際に行ったシミュレーションの結果］
次に、図１３～図１５を参照して、図１２に示したシミュレーション方法を用いてシミュレーションを行った結果について説明する。このシミュレーションには、式（４６）～式（５１）の支配方程式を用い、平行平板間の流れの解析を行った。

【0109】

図１３は、解析モデルの模式図である。解析空間１０を直方体の形状とした。流れの方向（図１３の左右方向）の寸法を０．０２５ｍとし、平行な一対の壁面１１の間隔（図１３の高さ方向の寸法）を０．１ｍとし、流れに直交する幅方向の寸法を０．０２５ｍとした。流れ方向及び幅方向の境界条件として、周期境界条件を課した。壁面１１の境界条件として、滑り無し条件を与えた。

【0110】

計算粒子２０のカーネル幅の初期値を０．００５ｍとし、複数の計算粒子２０を解析空間１０に等間隔に配置した。作動流体として空気を想定し、温度の初期値２９３Ｋ、圧力の初期値１０１．３ｋＰａを、解析空間１０内の全ての計算粒子２０に与えた。速度の初期条件として、平均流速が１７０ｍ／ｓの放物線分布になるように、壁面１１からの距離に応じて、計算粒子２０のそれぞれに速度の初期値を与えた。粘性係数はサザーランドの法則を用いて算出し、その１００倍の値を与えた。解析中は、流れ場の平均速度が１７０ｍ／ｓになるように流れ方向に圧力勾配を与えた。

【0111】

解析の時間刻み幅を７．２８×１０^－７ｓとした。時刻０ｓを初期状態とし、流れ場を時間発展させ、時刻０．１ｓにおける速度分布の値を出力した。また、比較対象として、解析に要する最大の空間解像度に合わせて、すべての計算粒子２０のカーネル幅を０．０００５ｍに固定した解析も行った。

【0112】

図１４Ａ及び図１４Ｂは、それぞれ計算開始時点及び計算終了時点における計算粒子２０の分布を示す模式図である。計算開始時点は、カーネル幅が０．００５ｍの複数の計算粒子２０が等間隔に整列している。計算終了時点では、壁面１１に近い領域ほど、細かい計算粒子２０が分布していることがわかる。これは、壁面１１の近傍では速度勾配が大きく、より高い空間解像度が必要なため、計算中に計算粒子２０が分割されてその大きさが小さくなっていったためである。計算終了時点では、計算粒子２０は約１２，０００個であった。なお、計算粒子２０のカーネル幅を０．０００５ｍに固定した比較例においては、計算粒子２０は約５００，０００個であった。

【0113】

図１５は、計算終了時点の速度分布及び初期状態の速度分布を示すグラフである。図１５の横軸は流速を単位［ｍ／ｓ］で表し、縦軸は、一方の壁面１１からの距離を単位［ｍ］で表す。図１５中の細い実線は初期状態の速度分布を示し、太い実線及び破線は、それぞれ本実施例及び比較例によるシミュレーション方法で求めた計算終了時点における速度分布を示す。

【0114】

本実施例によるシミュレーション方法で求めた速度分布は、比較例によるシミュレーション方法で求めた速度分布とほぼ同じであることがわかる。すなわち、本実施例では、比較例と比べて計算粒子２０の個数を約１／４０に低減させるとともに、流れ場の全域に亘って十分高い解像度で解析を行うことができることが確認された。

【0115】

上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。上記実施例では、粒子法としてＳＰＨ法を用いているが、その他の粒子法、例えばＭＰＳ法等を用いた流体解析にも、上記実施例による手法を適用することが可能である。

【符号の説明】

【0116】

１０ 解析空間
１１ 壁面
２０、２０ａ、２０ａ１～２０ａ８、２０ａｂ、２０ｂ、２０ｃ、２０ｄ、２０ｅ、２０ｓ 計算粒子
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】