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

Abstract

【課題】 高精度なシミュレーションを実現する。
【解決手段】 物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法において、（ａ）第１の時刻における複数の粒子の各々に作用する力に基づいて、第１の時刻の次の第２の時刻における、複数の粒子の位置を求める。（ｂ）第２の時刻における複数の粒子の位置をもとに、ボロノイ多面体を構成する。（ｃ）第１の時刻における複数の粒子の位置における物体の温度、及び、工程（ｂ）で構成されたボロノイ多面体の寸法と相対位置関係に基づいて、第２の時刻の複数の粒子の位置における物体の温度を求める。
【選択図】 図１

Description

本発明は、シミュレーション方法及びプログラムに関する。特に、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法を用いたシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムに関する。

コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。

また、繰り込み群の手法を応用した分子動力学法（繰り込み群分子動力学法）に基づく分子シミュレーションの発明が開示されている（たとえば、特許文献１参照）。

分子動力学法や繰り込み群分子動力学法を用いたシミュレーションにおいては、通常、温度は全体の粒子速度の分散から求められ、また熱は粒子の振動（格子振動）のみによって伝達されると考える。

しかしながら、銅などの良導体では電子が熱伝導に重要な役割を果たしているため、伝熱における自由電子の寄与を考慮しない上記シミュレーションにあっては、熱解析を正しく行うことは困難である。

特開２００６−２８５８６６号公報

本発明の目的は、高精度なシミュレーションを実現する、新規なシミュレーション方法を提供することである。

また、高精度なシミュレーションを実現する、新規なシミュレーション方法を実行させるプログラムを提供することである。

本発明の一観点によれば、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、（ａ）第１の時刻における複数の粒子の各々に作用する力に基づいて、前記第１の時刻の次の第２の時刻における、前記複数の粒子の位置を求める工程と、（ｂ）前記第２の時刻における前記複数の粒子の位置をもとに、ボロノイ多面体を構成する工程と、（ｃ）前記第１の時刻における前記複数の粒子の位置における前記物体の温度、及び、前記工程（ｂ）で構成されたボロノイ多面体の寸法と相対位置関係に基づいて、前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める工程とを有するシミュレーション方法が提供される。

また、本発明の他の観点によると、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、複数の粒子の初期状態を入力する手段、第１の時刻における複数の粒子の各々に作用する力に基づいて、前記第１の時刻の次の第２の時刻における、前記複数の粒子の位置を求める手段、前記第２の時刻における前記複数の粒子の位置をもとに、ボロノイ多面体を構成する手段、前記第１の時刻の前記複数の粒子の位置における前記物体の温度、及び、前記ボロノイ多面体の寸法と相対位置関係に基づいて、前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める手段、前記第２の時刻の前記複数の粒子の位置における前記物体の温度を表示する表示手段として機能させるためのプログラムが提供される。

本発明によれば、高精度なシミュレーションを実現する、新規なシミュレーション方法を提供することができる。

また、高精度なシミュレーションを実現する、新規なシミュレーション方法を実行させるプログラムを提供することができる。

（Ａ）は、実施例によるシミュレーション方法を示すフローチャートであり、（Ｂ）は、ボロノイ解析について説明するための図である。 実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。 （Ａ）〜（Ｃ）は、第１の計算系による検証について説明するための図である。 （Ａ）及び（Ｂ）は、第２の計算系による検証について説明するための図である。 （Ａ）及び（Ｂ）は、第３の計算系による検証について説明するための図である。 （Ａ）〜（Ｅ）は、第４の計算系による検証について説明するための図である。

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

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

図１（Ａ）は、実施例によるシミュレーション方法を示すフローチャートである。

実施例によるシミュレーション方法では、まずステップＳ１０１において、粒子の位置、速度、温度（粒子の位置における物体の温度）など、シミュレーションを行うに当たっての初期状態を決定する。次に、ステップＳ１０２において、ボロノイ(Voronoi)解析（ボロノイ分割）を行う。

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

再び、図１（Ａ）を参照する。ステップＳ１０２のボロノイ解析によって、ボロノイ多面体の寸法や相対位置関係、たとえばΔＳ_ｉｊ、ΔＶ_ｉ、ｒ_ｉｊを得ることができる。

ボロノイ解析後、ステップＳ１０３において、粒子間のポテンシャルや、重力場、電磁場等の外場が、各粒子に及ぼす力を、全ての粒子について導出する。更に、摩擦力等を考慮して計算することもできる。

続いて、ステップＳ１０４において、有限体積法(Finite Volume Method; FVM)を用い、温度場の解析を行う。温度場の解析においては、ステップＳ１０２のボロノイ解析で得られるｒ_ｉｊ、ΔＳ_ｉｊ、ΔＶ_ｉ、及び、一つ前の時刻（時間刻み幅Δｔだけ前の時刻）の温度場の情報を用いる。また、ボロノイ多面体の各々をＦＶＭの１つのコントロールボリュームにして解析する。

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

・・（１）

式（１）において、ρは密度、Ｃ_ｖは比熱、Ｔは温度、ｔは時間、Ｋは熱伝導率、Ｑは単位体積当たりの発熱量を示す。

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

・・（２）

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

・・（３）

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

・・（４）

を導くことができる。ここで、Ｔ_ｉ ^［ｎ］は、時刻ｎにおける粒子ｉの温度（粒子ｉが属するボロノイ多面体の温度）、Δｔは時間刻み幅（時刻ｎ＋１と時刻ｎとの間の時間間隔）、Ｋ_ｉｊは粒子ｉと粒子ｊとの間の熱伝導率を示す。式（４）により、ステップＳ１０２で得られたｒ_ｉｊ、ΔＳ_ｉｊ、ΔＶ_ｉ、及び、Ｋ_ｉｊと時刻ｎにおける粒子の温度、それから発熱がある場合には、時刻ｎ＋１における単位体積当たりの発熱量Ｑを用いて、時刻ｎ＋１における粒子の温度を求めることができる。

ステップＳ１０４においては、式（４）を用いて、ある時刻のすべての粒子の温度を算出する。なお、Ｔ_ｉ ^［０］には、ステップＳ１０１で与えられた初期状態での粒子ｉの温度を用いる。

なお、式（４）において、発熱量Ｑがゼロの場合には、ステップＳ１０３とＳ１０４の順序は前後してもよい。

ステップＳ１０５において、ステップＳ１０３で導出された、粒子ごとに働く力を用いて、全ての粒子につき運動方程式を解き、時間Δｔ後の各粒子の位置を数値計算で求め、各粒子を移動させる。

ステップＳ１０５で全ての粒子についての移動が終わったら、ステップＳ１０６にてシミュレーションを終了するか否かを決定する。終了を選択した場合はステップＳ１０７に進み、シミュレーションは完了する。終了しない場合は再びステップＳ１０２に戻り、ステップＳ１０２〜Ｓ１０６の工程を繰り返す。

なお、前述のように、Ｔ_ｉ ^［０］には、ステップＳ１０１で与えられた初期状態での粒子ｉの温度を用いるため、初期状態決定直後のボロノイ解析（ステップＳ１０２）は省くことが可能である。

更に、たとえば初期状態決定（Ｓ１０１）後、Ｓ１０３、Ｓ１０５、Ｓ１０２、Ｓ１０４、Ｓ１０６の順に工程を実施し、Ｓ１０６で終了しないと判定した場合は、Ｓ１０３に戻るシミュレーションとすることもできる。

実施例によるシミュレーション方法は、温度のパラメータをもった粒子を、粒子法（分子動力学法）を用いて移動させ、そのたびに移動後の粒子配置に基づいてボロノイ解析を実施し、ＦＶＭで温度場の解析を行う。粒子に直接温度のパラメータをもたせ、ＦＶＭ解析を行うことで、実際の比熱、熱伝導率を反映した解析が可能となる。実施例によるシミュレーション方法によれば、温度に関する高精度のシミュレーションを行うことができる。また、粒子の移動が生じるたびに、ボロノイ解析を実施するため、物体が変形する場合は勿論、たとえば流体のように、時間的に形状が著しく変化する対象を解析することも可能である。

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

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

まず、キーボードなどの入力装置から、図１（Ａ）のＳ１０１に対応して、粒子の位置、速度、温度などの粒子の初期状態が入力される。また、これと同時に、またはこれに先立って、キーボードなどの入力装置から、粒子間のポテンシャルや外場等を表現する関数を入力し、設定することができる。

中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリ中の粒子移動・温度解析プログラムを実行する。

粒子移動・温度解析プログラムは、たとえばボロノイ解析を行う部分（図１（Ａ）のＳ１０２で示される部分）、構成されたボロノイ多面体を用いて、ＦＶＭで温度場解析を実施する部分（図１（Ａ）のＳ１０４で示される部分）、粒子間のポテンシャルや外場等が各粒子に及ぼす力を、全ての粒子について導出し、粒子ごとに働く力を用いて、全ての粒子につき運動方程式を解き、時間Δｔ後の各粒子の位置を数値計算で求め、各粒子を移動させる部分（図１（Ａ）のＳ１０３及びＳ１０５で示される部分）に大別される。

図１（Ａ）を参照して説明したように、ある時刻における粒子位置に関する情報に基づいて、ボロノイ解析が行われる。

次に、粒子間ポテンシャルや外場等が各粒子に及ぼす力を、全ての粒子について導出する。

続いて、ボロノイ多面体を用い、ＦＶＭで温度場解析を実施する。

そして全ての粒子につき運動方程式を解き、時間Δｔ後の各粒子の位置を算出し、算出された新たな粒子位置に基づいてボロノイ解析を行い、時間Δｔ後の各粒子の温度を求める。

温度場解析の結果は、出力装置、たとえばディスプレイに、時間刻み幅Δｔごとに表示される。表示は、たとえば２次元的、または３次元的に定められた座標系において、粒子の位置と温度とを対応づけて行われる。たとえば、粒子位置に、色や濃淡で温度を示す形式で行われる。

本願発明者らは、実施例によるシミュレーション方法の妥当性を、以下に示すように４つの計算系で検証した。

図３（Ａ）〜（Ｃ）を参照して、第１の計算系による検証について説明する。

図３（Ａ）及び（Ｂ）に、２次元での温度場解析の計算系である、第１の計算系を示す。第１の計算系においては、まず、図３（Ａ）に示す長さ０．０５［ｍ］の短冊状領域の長さ方向（ｘ方向）に、下式（５）で表される初期（ｔ＝０）温度Ｔを与えた。

・・（５）

この場合において、短冊状領域の長さ方向の一方の端部（図３（Ａ）においては左端）をｘ＝０とし、他方の端部をｘ＞０の領域（ｘ＝０．０５）に配置した。なお、式（５）において、Ｌは温度分布をもたせる領域の長さであり、ここではＬ＝０．０５［ｍ］である。

図３（Ｂ）に、ｘ方向に沿う初期温度分布をグラフで示す。そして第１の計算系においては、時間が経過しても、短冊状領域の長さ方向の両端が常に０［Ｋ］となるように、吸熱を行わせることとした。

第１の計算系については、時刻ｔにおける温度分布の理論解が、下式（６）で与えられる。

・・（６）

ここでａは熱拡散率である。

図３（Ｃ）は、実施例によるシミュレーション方法で得られた温度分布と、式（６）で表される理論解とを比較するグラフである。グラフの横軸は、短冊状領域における、ｘ＝０を基点としたときのｘ方向に沿う長さ（ｘ座標）を単位「ｍ」で示し、グラフの縦軸は、温度を単位「Ｋ」で示す。比較は、２．０［ｓｅｃ］及び４．０［ｓｅｃ］経過時の温度分布について行った。「理論」の矢印で理論解を表し、「ＲＭＤ」の矢印で、実施例によるシミュレーション方法で得られた温度分布を表す。

実施例によるシミュレーション方法を用いた計算結果は、理論値とほぼ一致しており、誤差は０．５［％］以内であることがわかる。

図４（Ａ）及び（Ｂ）を参照して、第２の計算系による検証について説明する。第２の計算系による検証は、第１のそれと同様の検証を３次元モデルについて行うものである。

図４（Ａ）に、３次元での温度場解析の計算系である、第２の計算系を示す。第２の計算系には、一方向に長い（長さ０．０２［ｍ］）直方体をシミュレーション対象として用いた。直方体はＳＵＳで構成され、その密度は７８００［ｋｇ／ｍ^３］、比熱は５０．０［Ｊ／ｋｇ／Ｋ］、熱伝導率は２１．０［Ｗ／ｍ／Ｋ］、そして熱拡散率は５．３８×１０^−５［ｍ^２／ｓ］である。この対象の長さ方向をｘ方向として、式（５）で表される初期温度を与え、長さ方向の両端（ｘ＝０、０．０２［ｍ］）が常に０［Ｋ］となるように、吸熱を行わせた。

図４（Ｂ）は、実施例によるシミュレーション方法で得られた温度分布と、温度分布の理論値とを比較するグラフである。グラフの横軸は、直方体の長さ方向の一方端を基点としたときの、長さ方向の位置を単位「ｍ」で示し、グラフの縦軸は、温度を単位「Ｋ」で示す。比較は、０．２［ｓｅｃ］、０．５［ｓｅｃ］、及び１．０［ｓｅｃ］経過時の温度分布について行った。

実施例によるシミュレーション方法を用いた計算結果のグラフと、理論値によるグラフとには、画然として区別がつかないほどの一致が見られる。

図５（Ａ）及び（Ｂ）を参照して、第３の計算系による検証について説明する。第３の計算系による検証は、図５（Ａ）に示すように、原子（粒子）のサイズを大きくして、第２のそれと同様の検証を行うものである。原子（粒子）のサイズ以外の条件は、第２の計算系による検証と等しい。

図５（Ｂ）に示すように、第３の計算系による検証においても、実施例によるシミュレーション方法を用いた計算結果のグラフと、理論値によるグラフとには、画然として区別がつかないほどの一致が見られる。原子が大きい場合（メッシュが粗な場合）でも、原子が小さい場合（メッシュが密な場合）と同程度の精度で解析を行えることが確認された。

図６（Ａ）〜（Ｅ）を参照して、第４の計算系による検証について説明する。

図６（Ａ）に、第４の計算系を示す。第４の計算系は、荷重をかけた状態の円筒を、下地のブロックと擦りながら回転させる計算系である。回転方向は、図６（Ａ）において、反時計回りとした。

この場合において、円筒と下地ブロックの接触部分から発生する摩擦によるエネルギは、下式（７）で表される。

・・（７）

ここで、Ｆは摩擦力、ｖ_ｒは接触部粒子の相対速度、ΔＶは粒子の占める体積である。第４の計算系を用いて行う実施例によるシミュレーション方法においては、摩擦力Ｆを求め、式（７）を式（４）に代入することで、温度場を解析することになる。なお、摩擦力Ｆ、接触部粒子の相対速度ｖ_ｒ、及び粒子の占める体積ΔＶは、図１（Ａ）を参照して行った説明においては、すべて時刻ｎ＋１におけるものである。

図６（Ｂ）〜（Ｅ）に解析結果を示す。図６（Ｂ）〜（Ｅ）においては、温度の高低を表示の濃淡と対応づけて表している。図面に濃く表示されている領域が、温度の高い部分である。なお、シミュレーション結果は、時間を追って、図６（Ｂ）、（Ｃ）、（Ｄ）、（Ｅ）の順に表示してある。円筒と下地ブロックの接触部分で発生した熱が、徐々に伝達される様子が再現されている。

以上、第１〜第４の計算系による検証から、実施例によるシミュレーション方法の妥当性が確認された。

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

温度解析、伝熱解析を伴うシミュレーション一般に利用可能である。殊に、場の方程式の数値解法と、繰り込み群分子動力学法の連成解析に好適に応用可能であろう。

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

Claims (9)

1. 物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、
   （ａ）第１の時刻における複数の粒子の各々に作用する力に基づいて、前記第１の時刻の次の第２の時刻における、前記複数の粒子の位置を求める工程と、
   （ｂ）前記第２の時刻における前記複数の粒子の位置をもとに、ボロノイ多面体を構成する工程と、
   （ｃ）前記第１の時刻における前記複数の粒子の位置における前記物体の温度、及び、前記工程（ｂ）で構成されたボロノイ多面体の寸法と相対位置関係に基づいて、前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める工程と
   を有するシミュレーション方法。
2. 前記工程（ｃ）において、更に、前記第２の時刻に発生する熱量に基づいて、前記複数の粒子の位置における前記物体の温度を求める請求項１に記載のシミュレーション方法。
3. 前記熱量は前記第２の時刻における摩擦力により発生し、前記工程（ｃ）は、前記第２の時刻における摩擦力を求める工程を含む請求項２に記載のシミュレーション方法。
4. 前記工程（ｃ）において、前記ボロノイ多面体の各々を、有限体積法の１つのコントロールボリュームにして、前記複数の粒子の位置における前記物体の温度を求める請求項１〜３のいずれか１項に記載のシミュレーション方法。
5. 更に、前記工程（ｃ）の後に、
   （ｄ）ｋ＝２とする工程と、
   （ｅ）第ｋの時刻における複数の粒子の各々に作用する力に基づいて、前記第ｋの時刻の次の第ｋ＋１の時刻における、前記複数の粒子の位置を求める工程と、
   （ｆ）前記第ｋ＋１の時刻における前記複数の粒子の位置をもとに、ボロノイ多面体を構成する工程と、
   （ｇ）前記第ｋの時刻における前記複数の粒子の位置における前記物体の温度、及び、前記工程（ｆ）で構成されたボロノイ多面体の寸法と相対位置関係に基づいて、前記第ｋ＋１の時刻の前記複数の粒子の位置における前記物体の温度を求める工程と、
   （ｈ）ｎを２以上の自然数とするとき、ｋ＝ｎとなるまで、ｋの値を１ずつ増やして前記工程（ｅ）〜（ｇ）を繰り返す工程と
   を有する請求項１〜４のいずれか１項に記載のシミュレーション方法。
6. 物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、
   複数の粒子の初期状態を入力する手段、
   第１の時刻における複数の粒子の各々に作用する力に基づいて、前記第１の時刻の次の第２の時刻における、前記複数の粒子の位置を求める手段、
   前記第２の時刻における前記複数の粒子の位置をもとに、ボロノイ多面体を構成する手段、
   前記第１の時刻の前記複数の粒子の位置における前記物体の温度、及び、前記ボロノイ多面体の寸法と相対位置関係に基づいて、前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める手段、
   前記第２の時刻の前記複数の粒子の位置における前記物体の温度を表示する表示手段
   として機能させるためのプログラム。
7. 前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める手段は、更に、前記第２の時刻に発生する熱量に基づいて、前記複数の粒子の位置における前記物体の温度を求める請求項６に記載のプログラム。
8. 前記熱量は前記第２の時刻における摩擦力により発生し、前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める手段は、前記第２の時刻における摩擦力を求める請求項７に記載のプログラム。
9. 前記第２の時刻の前記複数の粒子の位置における前記物体の温度を求める手段は、前記ボロノイ多面体の各々を、有限体積法の１つのコントロールボリュームにして、前記複数の粒子の位置における前記物体の温度を求める請求項６〜８のいずれか１項に記載のプログラム。