JP2011175327A

JP2011175327A - シミュレーション方法及びプログラム

Info

Publication number: JP2011175327A
Application number: JP2010036999A
Authority: JP
Inventors: Ryota Hirose; 良太広瀬; Daiji Ichijima; 大路市嶋; Yoshitaka Onishi; 良孝大西
Original assignee: Sumitomo Heavy Industries Ltd
Current assignee: Sumitomo Heavy Industries Ltd
Priority date: 2010-02-23
Filing date: 2010-02-23
Publication date: 2011-09-08
Anticipated expiration: 2030-02-23
Also published as: JP5483342B2

Abstract

【課題】計算規模を小さくすることが可能なシミュレーション方法を提供する。
【解決手段】シミュレーション対象を粒子の集まりで表現し、ある時間刻み幅で、粒子に関して数値計算を行い、時間軸上で離散的に分布する時刻ごとの該シミュレーション対象の状態を解析する粒子法におけるシミュレーション方法であって、（ａ）第１の領域、第２の領域、及び第１の領域と第２の領域との境界領域を、それぞれ複数の粒子の集合で定義する工程と、（ｂ）第２の領域の粒子から境界領域の粒子に与えられる熱的寄与を加味して、第１の領域及び境界領域の粒子の温度を求める工程と、（ｃ）工程（ｂ）で求められた境界領域の粒子の温度を加味して、第２の領域の粒子の位置を求める工程とを有するシミュレーション方法を提供する。
【選択図】図１

Description

本発明は、シミュレーション方法及びプログラムに関する。特に、２つの領域のエネルギの授受を、２つの領域の境界に存在する粒子を媒介として行うシミュレーション方法、及びそのシミュレーションをコンピュータで実行するプログラムに関する。

有限要素法(finite elemental method; FEM)と分子動力学法(molecular dynamics method; MD)との連成解析が行われている（たとえば、非特許文献１参照）。しかしながら、ＦＥＭとＭＤとの連成解析や、有限体積法(finite volume method; FVM)とＭＤとの連成解析において、熱伝導による熱流束のやり取りを動的に行った例はない。

ＭＤのみによっても、たとえば温度場と構造変化を同時に解析することは可能である。しかし計算規模が極めて大きくなるため、現実的ではない。

「三次元分子動力学：有限要素法連結手法に関する基礎的検討」，原祥太郎，熊谷知久，泉聡志，酒井信介，計算力学講演会講演論文集２００４（１７），ｐｐ．１３７−１３８

本発明の目的は、新規なシミュレーション方法、及びそのシミュレーションを実行するプログラムを提供することである。

また、本発明の目的は、計算規模を小さくすることが可能なシミュレーション方法、及びそのシミュレーションを実行するプログラムを提供することである。

本発明の一観点によれば、シミュレーション対象を粒子の集まりで表現し、ある時間刻み幅で、粒子に関して数値計算を行い、時間軸上で離散的に分布する時刻ごとの該シミュレーション対象の状態を解析する粒子法におけるシミュレーション方法であって、（ａ）第１の領域、第２の領域、及び前記第１の領域と前記第２の領域との境界領域を、それぞれ複数の粒子の集合で定義する工程と、（ｂ）前記第２の領域の粒子から前記境界領域の粒子に与えられる熱的寄与を加味して、前記第１の領域及び前記境界領域の粒子の温度を求める工程と、（ｃ）前記工程（ｂ）で求められた前記境界領域の粒子の温度を加味して、前記第２の領域の粒子の位置を求める工程とを有するシミュレーション方法が提供される。

また、本発明の他の観点によれば、シミュレーション対象を粒子の集まりで表現し、ある時間刻み幅で、粒子に関して数値計算を行い、時間軸上で離散的に分布する時刻ごとの該シミュレーション対象の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、第１の領域、第２の領域、及び前記第１の領域と前記第２の領域との境界領域の粒子の初期状態を入力する入力手段、前記第２の領域の粒子から前記境界領域の粒子に与えられる熱的寄与を加味して、前記第１の領域及び前記境界領域の粒子の温度を求める手段、前記第１の領域及び前記境界領域の粒子の温度を求める手段で求められた前記境界領域の粒子の温度を加味して、前記第２の領域の粒子の位置を求める手段、少なくとも前記第１の領域の温度、前記第２の領域の粒子の位置のうちの一方を表示する表示手段として機能させるためのプログラムが提供される。

本発明によれば、新規なシミュレーション方法、及びそのシミュレーションを実行するプログラムを提供することができる。

また、本発明によれば、計算規模を小さくすることが可能なシミュレーション方法、及びそのシミュレーションを実行するプログラムを提供することができる。

（Ａ）〜（Ｄ）は、実施例によるシミュレーション方法を説明する図である。実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。（Ａ）〜（Ｄ）は、本願発明者らが行った温度解析シミュレーションについて説明する図である。（Ａ）〜（Ｆ）は、シミュレーション結果を示すグラフである。

以下、シミュレーション対象を粒子の集まりで表現し、ある時間刻み幅Δｔで各粒子に関して数値計算を行い、時間軸上で離散的に分布する時刻ごとのシミュレーション対象の状態を解析する粒子法を用いたシミュレーション方法、及びそのシミュレーションをコンピュータで実行するプログラムの実施例について説明する。

図１（Ａ）〜（Ｄ）を参照して、実施例によるシミュレーション方法を説明する。

図１（Ａ）は、実施例によるシミュレーション方法を示すフローチャートである。実施例によるシミュレーション方法では、まずステップＳ１０１において、シミュレーション対象、たとえばある物体につき、第１の領域、第２の領域、及び両領域の境界領域を、それぞれ複数の粒子の集合で定義する。ここで第１の領域と第２の領域とは、同じ材料で構成される領域でもよいし、異なる材料で構成される領域であってもよい。

図１（Ｂ）を用いてステップＳ１０１を詳説する。第１の領域とは、たとえば本図に示す下側の領域であり、有限体積法（ＦＶＭ）で温度場の計算を行う領域である。第２の領域とは、たとえば本図に示す上側の領域であり、分子動力学法（ＭＤ）を適用して、粒子の挙動を解析する領域である。そして境界領域とは、第１の領域と第２の領域の境界に存在する粒子から構成される領域である。

説明の便宜のため、第１の領域に属する粒子をＦＶＭ粒子、第２の領域に属する粒子をＭＤ粒子、境界領域に属する粒子をＰ粒子と呼称する。ＦＶＭ粒子は、たとえばＦＶＭ解析を行うメッシュの中心点に配置される粒子である。ＭＤ粒子については、ＭＤ解析が行われる。Ｐ粒子は、第１の領域と第２の領域の境界に存在し、ＦＶＭ粒子からの影響と、ＭＤ粒子からの影響の双方を受け、２つの領域のエネルギ授受を媒介する粒子である。なお、Ｐ粒子は、たとえば第１の領域の表面の粒子（ＭＤ粒子と隣接するＦＶＭ粒子）と考えることも可能である。

次に、図１（Ａ）のステップＳ１０２において、粒子の位置、速度、温度など、シミュレーションを行うに当たっての初期状態を決定する。

ステップＳ１０３においては、第２の領域の粒子から境界領域の粒子に与えられる熱的寄与、たとえば熱量を加味して、第１の領域及び境界領域で、ＦＶＭを用い温度場の解析を行う。第２の領域の粒子から境界領域の粒子に与えられる熱的寄与については、ステップＳ１０５を説明する際に触れる。第１の領域と境界領域の粒子には、温度のパラメータをもたせる。温度場解析には、たとえばボロノイ(voronoi)分割と熱伝導方程式を用いる。

図１（Ｃ）を参照して、ボロノイ分割について説明する。ボロノイ分割とは、粒子ｉと、その周囲に存在するすべての粒子とを線分で結び、その線分の垂直二等分面（２次元の場合は垂直二等分線）で領域を分ける分割方法をいう。ボロノイ分割によって形成された多面体をボロノイ多面体、多面体を構成する面をボロノイ面という。本明細書においては、粒子ｉｊ間の距離をｒ_ｉｊ、粒子ｉが含まれるボロノイ多面体をＶ_ｉ、その体積をΔＶ_ｉ、粒子ｉと粒子ｊを結ぶ線分の垂直二等分面の一部で構成されるボロノイ面（粒子ｉとｊとの間にあるボロノイ面）をＳ_ｉｊ、その面積をΔＳ_ｉｊと記述する。ボロノイ分割によって、ボロノイ多面体の寸法や相対位置関係、たとえばΔＳ_ｉｊ、ΔＶ_ｉ、ｒ_ｉｊを得ることができる。

ボロノイ分割で得られるｒ_ｉｊ、ΔＳ_ｉｊ、ΔＶ_ｉ、一つ前の時刻（時間刻み幅Δｔだけ前の時刻）の温度場の情報、及び第２の領域の粒子から境界領域の粒子に与えられる熱的寄与を用いて、現在時刻の温度場解析を行う。解析においては、ボロノイ多面体の各々をＦＶＭの１つのコントロールボリュームとする。

具体的には、たとえば下式（１）で示される熱伝導方程式を解く。

・・（１）

式（１）において、ρは第１の領域及び境界領域を構成する材料の密度、Ｃ_ｖは比熱、Ｋは熱伝導率、Ｔは温度、ｔは時間を示す。

(Qト゛ット)は、単位時間当たりの発熱量を示し、実施例においては、第２の領域の粒子から与えられる単位時間当たりの熱量を表す。

式（１）の両辺を体積積分して、下式（２）を得る。

・・（２）

更に、式（２）の右辺第１項にガウスの定理を適用すると、下式（３）が得られる。

・・（３）

そして、式（３）を離散化すると、ある時刻ｎについて、

・・（４）

を導くことができる。ここで、Ｔ_ｉ ^［ｎ］は、時刻ｎにおける粒子ｉの温度（粒子ｉが属するボロノイ多面体の温度）、Δｔは時間刻み幅（時刻ｎ＋１と時刻ｎとの間の時間間隔）、Ｋ_ｉｊは粒子ｉと粒子ｊとの間の熱伝導率を示す。式（４）により、ｒ_ｉｊ、ΔＳ_ｉｊ、ΔＶ_ｉ、Ｋ_ｉｊ、時刻ｎにおける粒子の温度、及び時刻ｎにおいて第２の領域の粒子から与えられる単位時間当たりの熱量

(Qト゛ット)を用いて、時刻ｎ＋１における粒子の温度を求めることができる。

ステップＳ１０３においては、式（４）を用い、第１の領域及び境界領域の各粒子について、現在時刻の温度を算出する。なお、Ｔ_ｉ ^［０］には、ステップＳ１０２で与えられた初期状態での粒子ｉの温度を用いる。

続いて、ステップＳ１０４において、境界領域の粒子（Ｐ粒子）の温度制御を行いながら、ＭＤを用い、第２の領域及び境界領域で粒子の挙動を解析する。解析によって、一例として第２の領域の粒子の位置が求められる。解析に当たっては、たとえば粒子間ポテンシャルによって第２の領域及び境界領域の各粒子に働く力、及び温度制御によって境界領域の粒子（Ｐ粒子）に働く力を求める。

粒子間ポテンシャルには、たとえば下式（５）で示すＭｏｒｓｅポテンシャルを用いる。

・・（５）

式（５）の両辺を微分することにより、粒子間ポテンシャルによって第２の領域及び境界領域の各粒子に働く力を求めることができる。なお、式（５）において、φ（ｒ）はポテンシャル、ｒは粒子間距離、Ｄはポテンシャルの深さ、Ａはポテンシャルの幅、ｒ_０は安定距離（ポテンシャルの底となる原子間距離）を示す。

Ｐ粒子に対する温度制御は、たとえばＬａｎｇｅｖｉｎ法を用いて行うことが可能である。Ｌａｎｇｅｖｉｎ法はＭＤにおいてフォノンの伝播速度にしたがって熱の授受を行い、さらに温度Ｔの熱浴の中にある粒子の運動を模擬する方法である。温度制御する粒子、すなわちＰ粒子に対して設定温度Ｔ、すなわちステップＳ１０３で求めたＰ粒子の温度となるように、次式（６）で表される標準偏差

・・（６）

をもつランダム力と減衰力とを算出して与える。これが温度制御によってＰ粒子に働く力となる。式（６）において、αは減衰力の減衰係数、ｋ_Ｂはボルツマン定数、Δｔは時間刻み幅である。さらに減衰係数αは、熱の授受がフォノンの伝播速度にしたがうようにするため、

・・（７）

とされる。ここでｍは粒子の質量、ω_Ｄはデバイ振動数である。

粒子間ポテンシャルによって第２の領域及び境界領域の各粒子に働く力、及び温度制御によってＰ粒子に働く力を求めた後、第２の領域及び境界領域の全ての粒子につき運動方程式を解き、各粒子の位置や速度を数値計算で求め、各粒子を移動させる。

また、このとき第２の領域の粒子の温度を速度の分散として求めてもよい。時刻ｎにおける粒子の速度は、たとえば時刻ｎ＋Δｔ／２における速度と、時刻ｎ−Δｔ／２における速度との相加平均として導出する。

このように、ステップＳ１０４においては、境界領域の粒子の温度を加味して、第２の領域及び境界領域の粒子の挙動を計算し、必要に応じてたとえば粒子の移動を行う。

ステップＳ１０４で第２の領域及び境界領域の粒子の挙動を解析した後、ステップＳ１０５において、第２の領域の粒子から境界領域の粒子に与えられる熱的寄与、たとえば熱量を算出する。

図１（Ｄ）を参照し、ステップＳ１０５について説明する。図１（Ｄ）に示す簡易モデルにおいて、粒子ａはＭＤ粒子（第２の領域の粒子）、粒子ｂはＰ粒子（境界領域の粒子）、バルク部分ｋはＦＶＭ粒子（第１の領域の粒子）全体に相当する。本モデルにおいては、たとえばＰ粒子を、第１の領域の表面の粒子（ＭＤ粒子と隣接するＦＶＭ粒子）と考える。

粒子ａに働く粒子ａ、ｂ間の力をｆ_ａ、粒子ａの速度をｖ_ａ、粒子ｂに働く粒子ａ、ｂ間の力をｆ_ｂ、粒子ｂ及びバルク部分ｋの速度をＶとすると、粒子ａが粒子ｂにする仕事の仕事率Ｐは、

・・（８）

と書ける。ここで粒子ｂが粒子ａにＰの仕事率で仕事をされたとき、された仕事のうちの半分が熱となり、半分が内部エネルギの変化に使われるとする。その場合、粒子ｂ内で発生する単位時間当たりの熱量は、

・・（９）

と表すことができる。粒子ａ、粒子ｂがそれぞれＭＤ粒子、Ｐ粒子であることを明瞭にするため、添え字を書き換えると、

・・（１０）

を得る。式（１０）において、ｆ_ＰはＰ粒子に働く、ＭＤ粒子とＰ粒子との間の力であり、ｖ_ＭＤはＭＤ粒子の速度、ＶはＦＶＭ粒子の速度である。そしてＱは、ＭＤ粒子からＰ粒子に与えられる熱量、

(Qト゛ット)は単位時間当たりにＭＤ粒子からＰ粒子に与えられる熱量を表す。

ステップＳ１０５においては、たとえば式（１０）を用いて、第２の領域の粒子（ＭＤ粒子）から境界領域の粒子（Ｐ粒子）に与えられる、単位時間当たりの熱量を算出する。式（１０）から明らかなように、これはたとえばＭＤ粒子がＰ粒子に対してする仕事を加味して計算される。

ステップＳ１０６にてシミュレーションを終了するか否かを決定する。終了を選択した場合は、シミュレーションは完了する。終了しない場合は、時刻を一つ進めて、再びステップＳ１０３に戻り、ステップＳ１０３〜Ｓ１０６の工程を繰り返す。

一つ進んだ時刻のステップＳ１０３において、式（１０）で算出された

(Qト゛ット)は、式（４）の

(Qト゛ット)に代入され、第１の領域及び境界領域の温度場の解析に使用される。

実施例によるシミュレーション方法は、ＦＶＭとＭＤとの連成解析において、ＦＶＭ解析を行う第１の領域と、ＭＤ解析を行う第２の領域との２つの領域の熱エネルギの授受を、両領域の境界に存在する粒子を媒介として行う、新規なシミュレーション方法である。境界粒子を媒介とすることで、第１の領域と第２の領域との連結部における熱エネルギの授受を、自然に再現することができる。実施例によるシミュレーション方法によれば、たとえばＦＶＭだけでは計算が困難な構造変化や相転移を伴う現象のシミュレーションを行うことが可能である。第２の領域の粒子については熱伝導に関し何の仮定も行っていないため、熱膨張、相変化等の影響を反映させることができる。

また、実施例によるシミュレーション方法は、計算規模を小さくすることが可能なシミュレーション方法である。ＭＤだけでは多大な計算時間を要する現象のシミュレーションを、小さい計算規模で行うことができる。たとえば表面にシリコン薄膜が形成されたガラス基板をシミュレーション対象とする場合、ガラス基板の厚さは、シリコン薄膜の厚さに比べて厚く、ＭＤで全体を計算するのは計算コストがかかる。しかし実施例によるシミュレーション方法を用いると、ガラス基板部分をＦＶＭで解析し、シリコン薄膜部分をＭＤで計算することによって、ガラス基板及びシリコン薄膜（シミュレーション対象）全体の温度場と、シリコン薄膜内の構造変化を、小さい計算規模でシミュレーションすることが可能である。

実施例によるシミュレーションは、プログラムによりコンピュータで実行することができる。

図２は、実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。本図に構成を示したシミュレーション装置を使用して、図１（Ａ）にフローチャートで示したシミュレーション方法を実施することができる。

まず、図１（Ａ）のＳ１０１に対応して、第１の領域、第２の領域、及び両領域の境界領域を、それぞれ複数の粒子の集合で定義する。各領域の定義、画定は、キーボードなどの入力装置から入力することで行われる。次に、入力装置から、Ｓ１０２に対応して、粒子の位置、速度、温度などの初期状態が入力される。

なお、第１の領域及び境界領域粒子の温度を計算する式、第２の領域の粒子から境界領域の粒子に与えられる、単位時間当たりの熱量を算出する式、粒子間のポテンシャルを表現する関数等は、これらと同時に、またはこれらに先立って、入力装置から入力することができる。

中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリ中のＦＶＭ・ＭＤ連成解析プログラムを実行する。

ＦＶＭ・ＭＤ連成解析プログラムは、たとえばＦＶＭを用い、第１の領域及び境界領域の温度場解析を行う部分（Ｓ１０３の処理を行う部分）、ＭＤを適用し、第２の領域及び境界領域で粒子の挙動を解析する部分（Ｓ１０４の処理を行う部分）、第２の領域の粒子から境界領域の粒子に与えられる、単位時間当たりの熱量を算出する部分（Ｓ１０５の処理を行う部分）に大別される。

図１（Ａ）を参照して説明したように、第２の領域の粒子から境界領域の粒子に与えられる熱的寄与を加味して、第１の領域及び境界領域で、ＦＶＭを用い温度場の計算を行う。

次に、境界領域の粒子の温度制御を行いながら、ＭＤを用い、第２の領域及び境界領域で粒子の挙動を解析する。解析に当たっては、たとえば粒子間ポテンシャルによって両領域の各粒子に働く力、及び温度制御によって境界領域の粒子に働く力を求める。

続いて、第２の領域の粒子から境界領域の粒子に与えられる単位時間当たりの熱量を算出する。算出された単位時間当たりの熱量は、次の時刻（時間刻み幅Δｔだけ進んだ時刻）のＦＶＭによる温度場解析に用いられる。

シミュレーションの結果は、出力装置、たとえばディスプレイに、時間刻み幅Δｔごとに表示される。たとえば１次元的、２次元的、または３次元的に定められた座標系を用い、少なくとも第１の領域の粒子の温度、第２の領域の粒子の位置のうちのいずれか一方が表示される。境界領域や第２の領域の粒子の温度を表示することもできる。表示は、たとえばグラフ形式で行われてもよい。

本願発明者らは、実施例によるシミュレーション方法の妥当性を検証した。図３（Ａ）〜（Ｄ）を参照し、本願発明者らが行った温度解析シミュレーションについて説明する。

図３（Ａ）にシミュレーション対象となるアルミニウム板を示す。アルミニウム板は、１６．１７９ｎｍの厚さを有する、ＸＹ面内方向に無限に広い板である。このアルミニウム板の下面に、２．００２×１０^１１Ｊ／ｍ^２／ｓの熱エネルギを均等に与える。熱エネルギの付与時間は１０ｐｓである。その後のエネルギ収支はゼロとする。この条件におけるアルミニウム板の、厚さ方向（Ｚ方向）の温度分布の時間変化を調べた。

図３（Ｂ）を参照する。温度解析に当たっては、図３（Ａ）に示すアルミニウム板から、Ｘ方向、Ｙ方向に沿う辺の長さがともに１．６１７ｎｍとなる直方体を厚さ方向に取り出し、この直方体について計算を行った。直方体の高さは、アルミニウム板の厚さと等しく、１６．１７９ｎｍである。熱は、直方体の下面に５．２３５×１０^−７Ｊ／ｓで１０ｐｓ間、付加されることになる。

図３（Ｃ）及び（Ｄ）を参照する。本願発明者らは、ＦＶＭとＭＤとを連成させた実施例によるシミュレーションと、ＦＶＭのみを用いた温度場解析の双方を行い、結果を比較した。

実施例によるシミュレーションにおいては、図３（Ｃ）に示すように、アルミニウム直方体を構成する粒子のうち、厚さ方向の中央に位置する粒子をＰ粒子、Ｐ粒子より下方（Ｚ軸負方向）に位置する粒子をＦＶＭ粒子、Ｐ粒子より上方（Ｚ軸正方向）に位置する粒子をＭＤ粒子とした。ＭＤ粒子相互間、及びＭＤ粒子とＰ粒子との間の粒子間ポテンシャルには、式（５）に示したＭｏｒｓｅポテンシャルを用い、アルミニウムのポテンシャルパラメータとして、式（５）においてＡ＝２．３５×１０^１０ｍ^−１、Ｄ＝１．９２×１０^−２０Ｊ、ｒ_０＝２．８６×１０^−１０ｍを使用した。

なお、図３（Ｄ）に示すように、アルミニウム直方体の温度場をＦＶＭのみで解析する場合、連続体の近似を直接用いるためノイズが少なく、解析値はこの系においては真値と考えることができる。

図４（Ａ）〜（Ｆ）は、シミュレーション結果を示すグラフである。グラフの横軸は、アルミニウム板（アルミニウム直方体）の厚さ方向の位置を、下面（Ｚ＝０）を基準として単位「ｎｍ」で示し、縦軸は、温度を単位「Ｋ」で示す。図４（Ａ）〜（Ｆ）は順に、熱エネルギを付与して、１ｐｓ後、１０ｐｓ後、２０ｐｓ後、５０ｐｓ後、７５ｐｓ後、１００ｐｓ後の温度分布を表す。各グラフにおいては、実施例によるシミュレーションによる解析結果を実線で、ＦＶＭによる解析結果を黒丸を連ねて表示した。

熱エネルギ付与中の温度分布を示す図４（Ａ）、熱エネルギ付与終了直後の温度分布を示す図４（Ｂ）を参照する。実施例によるシミュレーション（実線）においても、ＦＶＭによる解析と同様に、系の温度の上昇を示す分布が得られている。

熱エネルギ付与の終了後の温度分布を示す図４（Ｃ）〜図４（Ｆ）を参照する。どちらの解析手法によっても、時間の経過とともに、系全体の温度が一定値に近づく様子がシミュレートされている。

実施例によるシミュレーション方法を用いた計算の結果は、この場合真値と考えられるＦＶＭによる解析結果と傾向的に一致している。実施例によるシミュレーション方法の妥当性が確認されたといえるであろう。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

温度解析、伝熱解析を伴うシミュレーションに利用可能である。ＦＶＭのみによっては解析が難しい系、一例として物体が流体中に存する系（固体と液体が共存する系）等における温度場解析に適用することができる。また、ＭＤ粒子を冷媒と考えることで、冷凍機の解析へ展開することも可能である。

Ｓ１０１〜Ｓ１０６ステップ

Claims

シミュレーション対象を粒子の集まりで表現し、ある時間刻み幅で、粒子に関して数値計算を行い、時間軸上で離散的に分布する時刻ごとの該シミュレーション対象の状態を解析する粒子法におけるシミュレーション方法であって、
（ａ）第１の領域、第２の領域、及び前記第１の領域と前記第２の領域との境界領域を、それぞれ複数の粒子の集合で定義する工程と、
（ｂ）前記第２の領域の粒子から前記境界領域の粒子に与えられる熱的寄与を加味して、前記第１の領域及び前記境界領域の粒子の温度を求める工程と、
（ｃ）前記工程（ｂ）で求められた前記境界領域の粒子の温度を加味して、前記第２の領域の粒子の位置を求める工程と
を有するシミュレーション方法。
更に、
前記工程（ｃ）の後に、
（ｄ）前記第２の領域の粒子から前記境界領域の粒子に与えられる熱的寄与を算出する工程
を有し、
前記工程（ｂ）において、前記工程（ｄ）で算出された熱的寄与を加味して、前記第１の領域及び前記境界領域の粒子の温度を求め、
前記工程（ｂ）乃至（ｄ）を繰り返す請求項１に記載のシミュレーション方法。
前記工程（ｄ）において、
前記熱的寄与を、前記第２の領域の粒子が、前記境界領域の粒子に対してする仕事を加味して算出する請求項２に記載のシミュレーション方法。
前記工程（ｃ）において、
前記境界領域の粒子が前記工程（ｂ）で求めた温度となるように、前記境界領域の粒子に働く力を算出し、該算出された力を加味して、前記第２の領域の粒子の位置を求める請求項１〜３のいずれか１項に記載のシミュレーション方法。
前記工程（ｃ）において、
更に、粒子間ポテンシャルによって前記第２の領域の粒子に働く力を加味して、前記第２の領域の粒子の位置を求める請求項４に記載のシミュレーション方法。
前記工程（ｂ）において、
有限体積法を用いて、前記第１の領域及び前記境界領域の粒子の温度を求める請求項１〜５のいずれか１項に記載のシミュレーション方法。
前記工程（ｃ）において、更に、
前記工程（ｂ）で求められた前記境界領域の粒子の温度を加味して、前記第２の領域の粒子の速度を求め、前記第２の領域の粒子の温度を、該速度の分散として求める請求項１〜６のいずれか１項に記載のシミュレーション方法。
シミュレーション対象を粒子の集まりで表現し、ある時間刻み幅で、粒子に関して数値計算を行い、時間軸上で離散的に分布する時刻ごとの該シミュレーション対象の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、
第１の領域、第２の領域、及び前記第１の領域と前記第２の領域との境界領域の粒子の初期状態を入力する入力手段、
前記第２の領域の粒子から前記境界領域の粒子に与えられる熱的寄与を加味して、前記第１の領域及び前記境界領域の粒子の温度を求める手段、
前記第１の領域及び前記境界領域の粒子の温度を求める手段で求められた前記境界領域の粒子の温度を加味して、前記第２の領域の粒子の位置を求める手段、
少なくとも前記第１の領域の温度、前記第２の領域の粒子の位置のうちの一方を表示する表示手段
として機能させるためのプログラム。
更に、
前記第２の領域の粒子から前記境界領域の粒子に与えられる熱的寄与を算出する手段を備え、
前記第１の領域及び前記境界領域の粒子の温度を求める手段は、前記熱的寄与を算出する手段によって算出された熱的寄与を加味して、前記第１の領域及び前記境界領域の粒子の温度を求める請求項８に記載のプログラム。
前記熱的寄与を算出する手段は、前記熱的寄与を、前記第２の領域の粒子が、前記境界領域の粒子に対してする仕事を加味して算出する請求項９に記載のプログラム。
前記第２の領域の粒子の位置を求める手段は、前記境界領域の粒子が、前記第１の領域及び前記境界領域の粒子の温度を求める手段によって求められた温度となるように、前記境界領域の粒子に働く力を算出し、該算出された力を加味して、前記第２の領域の粒子の位置を求める請求項８〜１０のいずれか１項に記載のプログラム。
前記第２の領域の粒子の位置を求める手段は、更に、粒子間ポテンシャルによって前記第２の領域の粒子に働く力を加味して、前記第２の領域の粒子の位置を求める請求項１１に記載のプログラム。
前記第１の領域及び前記境界領域の粒子の温度を求める手段は、有限体積法を用いて、前記第１の領域及び前記境界領域の粒子の温度を求める請求項８〜１２のいずれか１項に記載のプログラム。
前記第２の領域の粒子の位置を求める手段は、更に、前記境界領域の粒子の温度を加味して、前記第２の領域の粒子の速度を求め、前記第２の領域の粒子の温度を、該速度の分散として求める請求項８〜１３のいずれか１項に記載のプログラム。