(57)【要約】【課題】固体壁からの熱負荷によって被加熱原子が加熱される加熱現象を、分子動力学法によって従来技術比較して小さい計算コストでかつ高精度でシミュレーションする。【解決手段】加熱シミュレーション装置１０は、加熱原子にてなる固体壁を温度制御することにより、分子動力学法を用いて、複数の面を有する所定の立体形状の計算領域内の被加熱原子を加熱する計算モデルを用いて加熱現象をシミュレーションする。加熱シミュレーション装置１０は、被加熱原子間のポテンシャルをレナード−ジョーンズポテンシャルによりシミュレーションし、立体形状の複数の面のうちの少なくとも１つの面に周期殻境界条件を適用する一方、他の面に周期境界条件を適用するように、分子動力学法を用いて被加熱原子の並進及び速度を計算してその位置及び温度を計算し、加熱原子の作用と被加熱原子の作用を異なる時間ステップ間隔で加熱現象をシミュレーションする。【選択図】図８
2007219831.tif 000002
【特許請求の範囲】
【請求項１】
加熱原子にてなる固体壁を温度制御することにより、分子動力学法を用いて、複数の面を有する所定の立体形状の計算領域内の被加熱原子を加熱する計算モデルを用いて加熱現象をシミュレーションする加熱シミュレーション装置において、
上記被加熱原子間のポテンシャルをレナード−ジョーンズ（Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ）ポテンシャルによりシミュレーションし、上記立体形状の複数の面のうちの少なくとも１つの面に周期殻境界条件を適用する一方、上記立体形状の複数の面のうちの他の面に周期境界条件を適用するように、分子動力学法を用いて上記被加熱原子の並進及び速度を計算してその位置及び温度を計算し、上記加熱原子の作用と上記被加熱原子の作用を異なる時間ステップ間隔で加熱現象をシミュレーションする制御手段を備えたことを特徴とする加熱シミュレーション装置。
【請求項２】
上記制御手段は、上記被加熱原子が１時間ステップ間隔で進行する間に、上記加熱原子が２時間ステップ間隔で進行するように上記加熱現象をシミュレーションすることを特徴とする請求項１記載の加熱シミュレーション装置。
【請求項３】
上記制御手段は、上記加熱シミュレーションの前の定常計算シミュレーションの後に、上記固体壁の温度を上昇させて上記加熱現象をシミュレーションし、
上記定常計算シミュレーションにおいて、上記計算モデルが所定の温度で平衡状態となり、上記被加熱原子を所定の設定温度に保持するように、毎時間ステップ毎に速度スケーリング法を用いて上記被加熱原子の速度に所定の係数を乗算すること特徴とする請求項１又は２記載の加熱シミュレーション装置。
【請求項４】
加熱原子にてなる固体壁を温度制御することにより、分子動力学法を用いて、複数の面を有する所定の立体形状の計算領域内の被加熱原子を加熱する計算モデルを用いて加熱現象をシミュレーションする加熱シミュレーション方法において、
上記被加熱原子間のポテンシャルをレナード−ジョーンズ（Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ）ポテンシャルによりシミュレーションし、上記立体形状の複数の面のうちの少なくとも１つの面に周期殻境界条件を適用する一方、上記立体形状の複数の面のうちの他の面に周期境界条件を適用するように、分子動力学法を用いて上記被加熱原子の並進及び速度を計算してその位置及び温度を計算し、上記加熱原子の作用と上記被加熱原子の作用を異なる時間ステップ間隔で加熱現象をシミュレーションする制御ステップを含むことを特徴とする加熱シミュレーション方法。
【請求項５】
上記制御ステップは、上記被加熱原子が１時間ステップ間隔で進行する間に、上記加熱原子が２時間ステップ間隔で進行するように上記加熱現象をシミュレーションすることを特徴とする請求項４記載の加熱シミュレーション方法。
【請求項６】
上記制御ステップは、上記加熱シミュレーションの前の定常計算シミュレーションの後に、上記固体壁の温度を上昇させて上記加熱現象をシミュレーションし、
上記定常計算シミュレーションにおいて、上記計算モデルが所定の温度で平衡状態となり、上記被加熱原子を所定の設定温度に保持するように、毎時間ステップ毎に速度スケーリング法を用いて上記被加熱原子の速度に所定の係数を乗算すること特徴とする請求項４又は５記載の加熱シミュレーション方法。
【請求項７】
請求項４乃至６のうちのいずれか１つに記載の加熱シミュレーション方法の各ステップを含むことを特徴とする加熱シミュレーションプログラム。
【請求項８】
請求項７記載の加熱シミュレーションプログラムを格納したことを特徴とする、コンピュータにより読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学を用いた分子の加熱シミュレーション装置及び方法、並びに加熱シミュレーションプログラム及びそれを格納した記録媒体に関する。特に、本発明は、固体壁から液体に熱負荷を与えたときの沸騰現象の分子動力学の計算を用いた加熱シミュレーション装置及び方法、並びに加熱シミュレーションプログラム及びそれを格納した記録媒体に関する。
【背景技術】
【０００２】
沸騰現象における初期気泡核の生成過程を把握することは、沸騰開始タイミングや初期ボイド率を正確に予測する上で重要である。近年、分子動力学法に基づくシミュレーションによって、気泡核生成プロセスを原子レベルで解明する試みが行われているが（例えば、非特許文献１参照。）、計算対象となる原子数又は計算ステップ数の制約から、プール沸騰など単純な加熱モデルをシミュレーションする場合においても、適切な境界条件を導入し、計算領域を合理的に設定する必要がある。
【０００３】
従来の分子動力学法によるシミュレーションにおいて、熱負荷を受けた液体分子の沸騰を計算する場合、まず、固体壁を構成する分子を液体分子に接するように配置し、固体壁分子に運動エネルギーを付加することによって液体への熱負荷をシミュレーションする。一方、固体壁に相対する計算領域端には周期境界条件を用いることができないため、等温壁等の境界条件が与えられる。この場合、等温壁における分子の流入や流出がシミュレーションされないため、一般的な沸騰現象にみられるような、固体壁面からバルク流体を加熱させるシミュレーションを行うためには、非現実的な大きさの計算モデルに対する分子動力学計算を実行する必要がある。そこで、熱負荷液体からの沸騰シミュレーションを実現可能な程度の小さな計算領域でシミュレーションするためには、シミュレーションの手法の改良が必要である。
【０００４】
例えば、非特許文献１においては、非常に薄い液膜を固体壁により急加熱し、沸騰させるといった現象を分子動力学法によりシミュレーションしており、その結果、加熱された流体原子の挙動は、マクロな流体挙動に類似した傾向を示している。すなわち、固体壁温度を１１０Ｋから３００Ｋに増加させると、液中に気泡が発生した後、液膜がライデンフロスト現象（又は膜沸騰）のように浮き上がり、やがて全液膜が気相へと蒸散したことを確認している。逆に、固体壁温度を１５０Ｋから１１０Ｋに下げた場合、凝縮により壁面上に液膜が発生することを示している。
【０００５】
また、非特許文献２において、固体壁で流体層を挟んだ体系において、固体壁間距離を増加させることにより圧力を低下させ、減圧沸騰による気泡核生成シミュレーションを行っており、壁面原子と液体原子のポテンシャルを変化させることにより固体壁の濡れ性を変化させている。その結果、壁面が濡れやすい程高い過熱度が必要になり、急激な核生成となることを示している。また、平衡状態での気泡の形状から接触角を測定し、壁−気泡間の関係が液滴の場合と同じになることを示している。なお、非特許文献３において同様のことが開示されている。
【０００６】
さらに、非特許文献４において、物体周りの分子動力学シミュレーションのための外部境界条件として、従来用いられてきた境界条件（周期境界条件や一様流条件）とは異なる周期殻境界条件の提案を行っている。円柱周りの流れに対して、周期殻境界条件、一様流条件及び自由流出条件の３通りの境界条件を適用した解析による加熱シミュレーションを行い、結果を比較している。周期殻境界条件を適用することにより、計算領域を小さく取った場合に生じる密度分布や温度分布の歪みを、他の境界条件を適用した場合よりも削減できることを示している。なお、非特許文献５において同様のことが開示されている。
【０００７】
また、非特許文献５において、温度勾配を有する壁面上での液滴の挙動を分子動力学法により調べており、半円柱状及び半円球状ともに低温度方向への移動が確認され、濡れ性が高く、接触角が小さい場合に移動速さが大きくなることを発見したことを開示している。
【０００８】
さらに、特許文献１においては、温度、速度、形状変化を並列演算処理の電子計算機を用いて解析する場合に、分子数が大きい場合でも計算時間を短くするような半導体製造プロセスシミュレーションなどに用いる分子動力学法の熱流体解析方法を提供するために、計算対象の中に複数の小さな部分領域を設け、部分領域で分子動力学法による詳細な解析を行い、その他の領域で比例的な推算計算を行い計算量を低減し、又は、計算対象のうち連続体として扱うことができない領域で分子動力学法による解析を行い、その他の領域で複数の分子の平均的な時間変化をマクロ解析により計算し、計算量を低減すること」が開示されている。
【０００９】
【特許文献１】特開平１０−０９７５１８号公報。
【非特許文献１】Pan Yi et al., "Molecular dynamics simulations of vaporization of an ultra-thin liquid argon layer on a surface", International Journal of Heat and Mass Transfer, No.45, pp.2087-2100, 2002.
【非特許文献２】Shigeo Maruyama et al., "A molecular dynamics simulations of a bubble nucleation on a solid surface", Heat and Technology, Calore e Technologia, Vol.l8, pp.69-73, 2000.
【非特許文献３】丸山茂夫ほか，“分子動力学法による固体壁面上での気泡核生成”，日本機械学会論文集（Ｂ編），６５巻，６３８号，１９９９年１０月。
【非特許文献４】Akira SATOH, "A new outer boundary condition for molecular dynamics simulations and its application to rarefied gas flow past a cylinder (periodic-shell boundary condition)", Advanced Powder Technology, Vol. 5, No.2, pp.105-117, 1994.
【非特許文献５】佐藤明，“分子動力学シミュレーションのための新しい外部境界条件の開発（周期殻境界条件）”，日本機械学会論文集（Ｂ編），５８巻，５５６号，１９９２年１２月。
【非特許文献６】村田章ほか，“温度勾配を有する壁面上での微小液滴現象の分子動力学シミュレーション”，第４２回日本伝熱シンポジウム公演論文集，ｐｐ．１９９−２００，２００５年６月。
【非特許文献７】神山新一ほか著，“分子動力学シミュレーション”，朝倉書店，分子シミュレーション講座２，ｐｐ．１５−１６，ｐｐ．２９−３０，ｐｐ．４７−４８，ｐｐ．５１−５５，ｐｐ．１３７，１９９７年５月１０日発行。
【非特許文献８】神山新一ほか著，“流体ミクロ・シミュレーション”，朝倉書店，分子シミュレーション講座３，ｐｐ．５２−６１，１９９７年５月１０日発行。
【発明の開示】
【発明が解決しようとする課題】
【００１０】
非特許文献１では、加熱による液体の沸騰現象を解析している。しかし、加熱壁と相対する計算領域の境界に、原子の流入出のない断熱壁（弾性散乱壁）を適用しているため、計算領域をすべて液相で満たすと圧力上昇によって気泡が生じなくなる。そこで、極めて薄い液膜の上に気相を設けた系を扱うこととなり、これにより、非特許文献１の方法による実現象への適用範囲は限定されたものとなっている。
【００１１】
また、非特許文献２及び３では、沸騰のシミュレーションを行っているが、固体壁により流体を挟んだ計算モデルを扱っている。このような対称性のよい計算モデルは比較的容易に生成できるが、得られた結果は、液相の対流などの点において加熱沸騰とは異なるものと考えられる。よって、非特許文献２及び３での計算は、極めて小さな系の減圧沸騰現象をシミュレーションしたものとなっている。
【００１２】
さらに、非特許文献４及び５では、計算領域の内外に比較的広範に広がる流体の取り扱いに対して、新しい境界条件を提案している。しかし、その適用例は円柱周り流に限られており、片側から流体を加熱するような非対称な計算領域への適用について触れられていない。
【００１３】
また、非特許文献６では、特に沸騰現象を扱っているわけではないが、固体壁に相対する境界に断熱壁を用いた解析の例として取り上げている。計算モデルとしては、非特許文献１によく似ている。
【００１４】
さらに、特許文献１では、温度、速度、形状変化を並列演算処理の電子計算機を用いて解析する場合について開示しているが、液体の沸騰現象を取り扱うことはできないという問題点があった。
【００１５】
本発明の目的は以上の問題点を解決し、固体壁からの熱負荷によって被加熱原子が加熱される加熱現象を、分子動力学法によって従来技術比較して小さい計算コストでかつ高精度でシミュレーションすることができる加熱シミュレーション装置及び方法、並びに加熱シミュレーションプログラム及びそれを格納した記録媒体を提供することにある。
【課題を解決するための手段】
【００１６】
第１の発明に係る加熱シミュレーション装置は、加熱原子にてなる固体壁を温度制御することにより、分子動力学法を用いて、複数の面を有する所定の立体形状の計算領域内の被加熱原子を加熱する計算モデルを用いて加熱現象をシミュレーションする加熱シミュレーション装置において、
上記被加熱原子間のポテンシャルをレナード−ジョーンズ（Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ）ポテンシャルによりシミュレーションし、上記立体形状の複数の面のうちの少なくとも１つの面に周期殻境界条件を適用する一方、上記立体形状の複数の面のうちの他の面に周期境界条件を適用するように、分子動力学法を用いて上記被加熱原子の並進及び速度を計算してその位置及び温度を計算し、上記加熱原子の作用と上記被加熱原子の作用を異なる時間ステップ間隔で加熱現象をシミュレーションする制御手段を備えたことを特徴とする。
【００１７】
また、上記加熱シミュレーション装置において、上記制御手段は、上記被加熱原子が１時間ステップ間隔で進行する間に、上記加熱原子が２時間ステップ間隔で進行するように上記加熱現象をシミュレーションすることを特徴とする。
【００１８】
さらに、上記加熱シミュレーション装置において、上記制御手段は、上記加熱シミュレーションの前の定常計算シミュレーションの後に、上記固体壁の温度を上昇させて上記加熱現象をシミュレーションし、
上記定常計算シミュレーションにおいて、上記計算モデルが所定の温度で平衡状態となり、上記被加熱原子を所定の設定温度に保持するように、毎時間ステップ毎に速度スケーリング法を用いて上記被加熱原子の速度に所定の係数を乗算すること特徴とする。
【００１９】
第２の発明に係る加熱シミュレーション方法は、加熱原子にてなる固体壁を温度制御することにより、分子動力学法を用いて、複数の面を有する所定の立体形状の計算領域内の被加熱原子を加熱する計算モデルを用いて加熱現象をシミュレーションする加熱シミュレーション方法において、
上記被加熱原子間のポテンシャルをレナード−ジョーンズ（Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ）ポテンシャルによりシミュレーションし、上記立体形状の複数の面のうちの少なくとも１つの面に周期殻境界条件を適用する一方、上記立体形状の複数の面のうちの他の面に周期境界条件を適用するように、分子動力学法を用いて上記被加熱原子の並進及び速度を計算してその位置及び温度を計算し、上記加熱原子の作用と上記被加熱原子の作用を異なる時間ステップ間隔で加熱現象をシミュレーションする制御ステップを含むことを特徴とする。
【００２０】
また、上記加熱シミュレーション方法において、上記制御ステップは、上記被加熱原子が１時間ステップ間隔で進行する間に、上記加熱原子が２時間ステップ間隔で進行するように上記加熱現象をシミュレーションすることを特徴とする。
【００２１】
さらに、上記加熱シミュレーション方法において、上記制御ステップは、上記加熱シミュレーションの前の定常計算シミュレーションの後に、上記固体壁の温度を上昇させて上記加熱現象をシミュレーションし、
上記定常計算シミュレーションにおいて、上記計算モデルが所定の温度で平衡状態となり、上記被加熱原子を所定の設定温度に保持するように、毎時間ステップ毎に速度スケーリング法を用いて上記被加熱原子の速度に所定の係数を乗算すること特徴とする。
【００２２】
第３の発明に係る加熱シミュレーションプログラムは、上記加熱シミュレーション方法の各ステップを含むことを特徴とする。
【００２３】
第４の発明に係るコンピュータにより読み取り可能な記録媒体は、上記加熱シミュレーションプログラムを格納したことを特徴とする。
【００２４】
本発明では、上記目的を達成するため、固体壁に相対する計算領域端に周期殻境界条件を適用し、殻領域内の液体分子の運動及び配置の情報を計算領域外に複写することにより、計算領域外にも液体が存在するような分子間力のシミュレーションを行う。具体的には、ハーモニック・ポテンシャル原子によって構成される固体壁を計算領域の片側に配置し，それに相対する計算領域端に殻領域を設定する。殻領域の厚さは、液体分子の分子間力計算におけるカットオフ距離程度とする。各計算ステップにおいて、殻領域に存在する液体原子を、位置座標のみ計算領域外にスライドするように複写し、分子間力の積分計算において、複写された液体分子が計算領域内の液体分子に及ぼす寄与も含めるようにする。周期殻境界条件を適用した面における分子の流入及び流出は、殻領域に存在する分子、又はその複写分子が、殻領域端又は計算領域端をまたいで移動した場合に自発的に発生する。なお、壁面の温度上昇は、固体壁と結合したランジュバン方程式に基づくファントム原子に与えるランダム力の時間ステップ間隔を増加させることによりシミュレーションする。
【発明の効果】
【００２５】
従って、本発明に係る加熱シミュレーション装置及び方法によれば、固体壁からの熱負荷によってもたらされる被加熱原子の加熱現象を、従来技術に比較して小規模な計算モデルで計算し、計算コストを大幅に削減できるとともに、従来技術に比較して高精度で加熱現象をシミュレーションでき、相変化冷却などのシミュレーションを行うことが可能となる。
【発明を実施するための最良の形態】
【００２６】
以下、本発明に係る実施形態について図面を参照して説明する。なお、以下の各実施形態において、同様の構成要素については同一の符号を付している。
【００２７】
図１は本発明の一実施形態に係る、分子動力学による高熱負荷液体からの加熱シミュレーション装置１０の構成を示すブロック図であり、図２は図１の加熱シミュレーション装置１０において用いる計算モデルを示す図である。
【００２８】
本発明の実施形態は、固体壁から加熱される被加熱原子群の状態を、分子動力学法を用いて、比較的小さな計算モデルで解析してシミュレーションすることを目的としている。これを実現するために、図２に示すような直交座標系に設定した、例えば直方体から成る計算モデルにおいて、まず４つの側面に周期境界条件を適用し、残る２面のうち底面を固体原子から構成される加熱の固体壁面とし、さらに底面と相対する上面を周期殻境界条件とするような境界条件の組み合わせを提案した。それぞれの境界条件はすでに論文発表されたものであるが、このような境界条件の組み合わせによって、従来の解析手法よりも大幅に少ない記憶容量と計算時間で流体の加熱シミュレーションを実施することが可能となることを特徴としている。
【００２９】
すなわち、本発明に係る実施形態では、固体壁と相対する境界に周期殻境界条件を適用することにより、計算領域の片側から流体が加熱される計算モデルを生成し，固体壁からの熱伝達によって流体が沸騰する現象を分子動力学法により解析して加熱シミュレーションすることができる加熱シミュレーション装置１０について以下に説明する。
【００３０】
図２の計算モデルにおいて、３ｎｍ×３ｎｍ×６ｎｍの直方体内に６８６個の液体アルゴン原子（以下、Ａｒ原子という。）を配置し、当該直方体の一方の面に加熱面となる固体壁を配した。当該固体壁は５層から成り、１層あたり１２８個の白金原子（以下、Ｐｔ原子という。）を第１乃至第３に配置し、第４層に温度制御のためのファントム（ｐｈａｎｔｏｍ）原子、また、第５層に壁面を固定する固定原子を配置した。ここで、Ａｒ原子間及びＰｔ原子間の原子間力は、それぞれ次式に係るレナード−ジョーンズ（Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ）ポテンシャル関数（以下、Ｌ−Ｊポテンシャル関数という。なお、レナード−ジョーンズポテンシャルをＬ−Ｊポテンシャルという。）φＬＪと、ハーモニック（ｈａｒｍｏｎｉｃ）ポテンシャル関数φＨの空間導関数から計算される。
【００３１】
【数１】
2007219831.tif 000003
【数２】
2007219831.tif 000004
【００３２】
ここで、ｒｉｊは原子間距離を表し、Ｌ−Ｊポテンシャル関数の井戸の深さε＝１．６７×１０−２１Ｊ）、ポテンシャルエネルギーが零となる原子間距離（すなわち、実質的な原子直径σ（＝０．３４ｎｍ）、バネ定数ｋ（＝４６．８Ｎ／ｍ）、Ｐｔ原子の最近接原子間距離δＰｔ（＝０．３４ｎｍ）はポテンシャルパラメータである。Ａｒ原子とＰｔ原子との相互作用は、Ｌ−Ｊポテンシャル関数に従うとし、Ｌ−Ｊポテンシャル関数の井戸の深さεを０．８９４×１０−２１Ｊとして計算した。
【００３３】
運動方程式の差分には、ベロシティ・ベルレット法（Ｖｅｌｏｃｉｔｙ  Ｖｅｒｌｅｔ法）（例えば、非特許文献７参照。）を用いた。Ｌ−Ｊポテンシャルを持つＡｒ原子は、時間の差分精度がＰｔ原子と異なるため、Ａｒ原子とＰｔ原子とで異なる時間ステップ間隔を与えるマルチプル・タイムステップ（Ｍｕｌｔｉｐｌｅ  Ｔｉｍｅ  Ｓｔｅｐ）法（例えば、非特許文献１参照。）により運動方程式の差分展開を行った。なお、Ａｒ原子の原子間力の計算において、Ａｒ原子直径の３．８倍のカットオフ距離を適用した。固体壁の温度は、ファントム原子にランダムな力を与えることによって制御される。ファントム原子の運動は３次元のランジュバン（Ｌａｎｇｅｖｉｎ）方程式により記述され、ランダム力は、次式の標準偏差σｒｍｓの正規分布に従うとしている（例えば、非特許文献１及び３参照。）。
【００３４】
【数３】
2007219831.tif 000005
【００３５】
ここで、Ｔは壁面温度、ｋＢはボルツマン定数、Δｔは解析シミュレーションの時間ステップ間隔であり、ξはファントム原子に与える減衰力の係数である。
【００３６】
次いで、図２の計算モデルの境界条件及び解析シミュレーション方法の概要について以下に説明する。
【００３７】
図２の計算モデルは下端の面に固体壁を有し、境界条件は、計算領域の側面（４面）を周期境界条件とし、計算領域上端面にのみ、周期殻境界条件（例えば、非特許文献４及び５参照。）を適用した。計算領域の上方に厚さ１．２９ｎｍの殻領域を設け、時間ステップ間隔毎に殻領域内の原子を計算領域外へ複写する。これより、計算領域上端に存在する原子に対し、あたかも領域外に原子が存在するような原子間力が与えられる。解析シミュレーションは、まず、計算モデルが温度１００Ｋで平衡状態となるように速度スケーリングをかけ、時刻ｔ＝４ｐｓｅｃまで計算した。その後、ファントム原子に与えるランダム力の間隔を０．４ｐｓｅｃから０．１ｐｓｅｃに増加させ、固体壁の壁面温度を上昇させた。ここで、時間ステップ間隔をΔｔ＝０．２×１０−３ｐｓｅｃとし、１０ｐｓｅｃまでの解析シミュレーションを行った。なお、周期殻境界条件の適用により、計算領域上部からのＡｒ原子の流入出が生じるが、加熱直前のＡｒ原子数は６８１個であり、初期状態とほぼ同じであった。なお、本解析シミュレーション（約３．５ＭＰａに相当）におけるＡｒ原子の飽和温度は１４６Ｋである。
【００３８】
次いで、図１を参照して、本実施形態に係る加熱シミュレーション装置１０の構成について以下に説明する。
【００３９】
本実施形態に係る加熱シミュレーション装置１０は、加熱原子にてなる固体壁を温度制御することにより、分子動力学法を用いて、複数の面を有する所定の立体形状の計算領域内の被加熱原子を加熱する計算モデルを用いて加熱現象をシミュレーションする。ここで、加熱シミュレーション装置１０は、上記被加熱原子間のポテンシャルをＬ−Ｊポテンシャルによりシミュレーションし、上記立体形状の複数の面のうちの少なくとも１つの面に周期殻境界条件を適用する一方、上記立体形状の複数の面のうちの他の面に周期境界条件を適用するように、分子動力学法を用いて上記被加熱原子の並進及び速度を計算してその位置及び温度を計算し、上記加熱原子の作用と上記被加熱原子の作用を異なる時間ステップ間隔で加熱現象をシミュレーションすることを特徴とする。
【００４０】
また、上記加熱シミュレーション装置１０は、上記被加熱原子が１時間ステップ間隔で進行する間に、上記加熱原子が２時間ステップ間隔で進行するように上記加熱現象をシミュレーションする。さらに、上記加熱シミュレーション装置１０は、上記加熱シミュレーションの前の定常計算シミュレーションの後に、上記固体壁の温度を上昇させて上記加熱現象をシミュレーションし、上記定常計算シミュレーションにおいて、上記計算モデルが所定の温度で平衡状態となり、上記被加熱原子を所定の設定温度に保持するように、毎時間ステップ毎に速度スケーリング法を用いて上記被加熱原子の速度に所定の係数を乗算する。
【００４１】
本実施形態の加熱シミュレーション装置１０は、図１に示すように、ディジタル計算機を含むように構成され、例えば図６乃至図１０の処理を実行することにより、図２の計算モデルの加熱シミュレーション処理を実行し、そのシミュレーション結果を表示して出力する。加熱シミュレーション装置１０は、
（ａ）当該加熱シミュレーション装置１０の動作及び処理を演算及び制御するコンピュータのＣＰＵ（中央演算処理装置）２０と、
（ｂ）オペレーションプログラムなどの基本プログラム及びそれを実行するために必要なデータを格納するＲＯＭ（読み出し専用メモリ）２１と、
（ｃ）ＣＰＵ２０のワーキングメモリとして動作し、図６の加熱シミュレーション処理と、図７の初期計算処理（ステップＳ１）と、図８の主計算処理（ステップＳ２）と、図９の複写層の生成処理（ステップＳ２６）におけるＡｒ原子の流入出の判別処理（ステップＳ２６（ｂ））と、図１０のＰｔ原子の並進処理（ステップＳ２７）とにおいて必要なパラメータやデータを一時的に格納するＲＡＭ（ランダムアクセスメモリ）２２と、
（ｄ）例えばハードディスクメモリで構成され、入力パラメータのデータやシミュレーション結果のデータなどのデータを格納するデータメモリ２３と、
（ｅ）例えばハードディスクメモリで構成され、ＣＤ−ＲＯＭドライブ装置４５を用いて読みこんだ図６乃至図１０の解析シミュレーション処理プログラムを格納するプログラムメモリ２４と、
（ｆ）所定のデータや指示コマンドを入力するためのキーボード４１に接続され、キーボード４１から入力されたデータや指示コマンドを受信して所定の信号変換などのインターフェース処理を行ってＣＰＵ２０に伝送するキーボードインターフェース３１と、
（ｇ）ＣＲＴディスプレイ４３上で指示コマンドを入力するためのマウス４２に接続され、マウス４２から入力されたデータや指示コマンドを受信して所定の信号変換などのインターフェース処理を行ってＣＰＵ２０に伝送するマウスインターフェース３２と、
（ｈ）ＣＰＵ２０によって処理されたデータや設定指示画面などを表示するＣＲＴディスプレイ４３に接続され、表示すべき画像データをＣＲＴディスプレイ４３用の画像信号に変換してＣＲＴディスプレイ４３に出力して表示するディスプレイインターフェース３３と、
（ｉ）ＣＰＵ２０によって処理されたデータ及び所定のシミュレーション結果などを印字するプリンタ４４に接続され、印字すべき印字データの所定の信号変換などを行ってプリンタ４４に出力して印字するプリンタインターフェース３４と、
（ｊ）図６乃至図１０の解析シミュレーション処理プログラムが記憶されたＣＤ−ＲＯＭ４５ａから解析シミュレーション処理プログラムのプログラムデータを読み出すＣＤ−ＲＯＭドライブ装置４５に接続され、読み出された解析シミュレーション処理プログラムのプログラムデータを所定の信号変換などを行ってプログラムメモリ２４に転送するドライブ装置インターフェース３５とを備える。
ここで、これらの回路２０−２４及び３１−３４はバス３０を介して接続される。
【００４２】
本実施形態に係る解析シミュレーション処理プログラムはＦＯＲＴＲＡＮによって生成され、当該プログラムのフローチャートを図６乃至図１０に示している。本実施形態に係る解析シミュレーション例では、被加熱原子をＡｒ原子とし、加熱原子をＰｔ原子としている。以下、本解析シミュレーション処理プログラムの内容について以下に説明する。なお、図１５は図１の加熱シミュレーション装置１０において用いる各変数の内容、値、単位及び式中の記号を示す図である。
【００４３】
解析シミュレーション処理プログラムの実行時の始めに、計算領域の大きさや原子数などの解析シミュレーション条件についてユーザがキーボード４１やマウス４２を用いて入力し、その条件に従って計算される原子の初期座標や初期速度を各配列変数に代入する。次に、各原子間の距離に基づいて、原子のポテンシャル関数から原子間力を計算し、各原子に作用する力の積分値を記憶する。また、差分化された運動方程式に各原子の座標、速度、力を代入することより、１時間ステップに相当する時間ステップ間隔Δｔ（ＤＴＩＭＥ）後の座標及び速度を計算し、それぞれのデータを格納する配列を上書きする。この一連の計算により、計算領域内の全原子が時間ステップ間隔Δｔだけ並進運動した状態の解を得ることになる。
【００４４】
計算領域の各面に適用した境界条件は次のように扱う。まず、計算モデルの側面に適用した周期境界条件は、全てのＡｒ及びＰｔ原子に対し、原子間力の積分計算及び原子座標の並進計算時に適用される。また、計算領域の上面に設定した周期殻境界条件では、まず計算領域上部に殻領域を設け、その中に存在するＡｒ原子を計算領域の外に複写する。さらに、Ａｒ原子の原子間力の積分値に、複写原子からの相互作用を含めるようにする。殻領域原子の複写は、Ａｒ原子が並進し、かつＡｒ原子の流入出を判断した後の座標値に対して行われ、複写された原子の座標は、Ａｒ原子の座標配列に追加される。
【００４５】
さらに、計算領域の底面に設定した加熱壁面では、Ｐｔ原子の並進運動を計算する。Ａｒ原子の座標配列が更新された後、マルチプル・タイムステップ法（例えば、非特許文献１参照。）により、Ｐｔ原子の並進計算が２回実行され、Ｐｔ原子の座標、速度、原子間力を記憶する配列の値が２回更新される。加熱壁面の温度は、Ｐｔ原子に接続したファントム原子にランダム力を与える時間ステップ間隔で制御する。具体的には、繰り返し回数ＩＲＥＰＴを参照することにより、ランダム力を与える時間ステップ間隔Δｔを非加熱定常状態の時間ステップ間隔ＮＦＯＣ１と、加熱状態の時間ステップ間隔ＮＦＯＣ２とで切り替える。加熱壁面から流体原子への熱の授受は、各原子間の原子間力によって評価される。Ａｒ原子とＰｔ原子間の原子間力はＡｒ原子の座標とＰｔ原子の座標から原子間距離を計算し、その距離をポテンシャル関数に代入することによって計算される。これらの値は、各原子の原子間力の積分値（ＡｒはＦＸ，ＰｔはＷＦＸ）に加算され、次の並進運動に寄与し、結果として運動エネルギーの増減をもたらす。
【００４６】
以上により、１時間ステップ間隔のステップの計算処理が完了する。必要に応じて各原子の座標や速度の値を計算しながら計算が進行し、必要とされる時間ステップだけ繰り返されると計算は終了する。
【００４７】
次いで、初期条件及び原子の初期配置について以下に説明する。
【００４８】
まず、ｘｙｚ座標系を図２に示したように設定する。すなわち、ｘｙｚ座標系の原点を、当該計算モデルの直方体の下端右奧の頂点におく。ここで、解析シミュレーション前に入力する値は、計算モデルのｘ，ｙ方向の一辺の長さ（ＶＳＩＺＥ）、ｚ方向の一辺の長さ（ＺＶＳＩＺＥ）、Ａｒ原子及びＰｔ原子（ファントム原子を含む。）の個数（ＮＵＭ及びＷＮＵＭ）、Ａｒ原子及びＰｔ原子の初期温度（ＴＥＭＰ，ＷＴＥＭＰ）などである。設定圧力におけるＡｒ原子の密度（８４２．２ｋｇ／ｍ３）と計算領域の体積（ＶＳＩＺＥ＊ＶＳＩＺＥ＊ＺＶＳＩＺＥ）からＡｒ原子の原子数を計算し、さらに、Ａｒ原子が計算モデル内に立方格子状に初期配列されるようにＮＵＭ＝６８６個とした。ここで、＊は乗算を表す。Ｐｔ原子は、面心立方格子状に３層配置するため、ＷＮＵＭ＝６４０個（１層当たり１２８個）とした。なお、初期温度はＡｒ原子及びＰｔ原子ともにＴＥＭＰ＝ＷＴＥＭＰ＝１００Ｋとした。
【００４９】
原子の初期配置において、液体Ａｒ原子は、図２の計算モデル内に体心立方格子状に配置する。格子間隔Ａ（＝０．４２９ｎｍ）は計算モデルの体積とＡｒ原子個数から解析シミュレーション処理プログラムにより計算される。さらに、初期配置が計算結果に及ぼす影響を低減するため、Ａｒ原子の座標を格子点から乱数を用いて微小移動させている。すなわち、ｉ番目（ｉ＝１〜ＮＵＭ）のＡｒ原子の移動前の座標を（ｘｉ’，ｙｉ’，ｚｉ’）、移動後の座標を（ｘｉ，ｙｉ，ｚｉ）とすると、次式で表される。
【００５０】
【数４】
2007219831.tif 000006
【００５１】
ここで、Ｒｉ，Ｒｉ＋１，Ｒｉ＋２は［０，１］の任意の乱数（＝ＲＡ）であり、微小移動量Δｌ（ＤＬＴＡ）は格子間隔Ａの０．１倍とした。
【００５２】
解析シミュレーション処理プログラムにおいて、ｘｙｚ座標系の座標値（原子位置）を格納する配列は、（ｘｉ，ｙｉ，ｚｉ）＝（ＲＸ（ＮＵＭ），ＲＹ（ＮＵＭ），ＲＺ（ＮＵＭ））としている。なお、図６乃至図１０の加熱シミュレーション処理のフローチャート内、及び下記の説明では、特別な必要のない限り、３次元の座標配列（ＲＸ（ＮＵＭ），ＲＹ（ＮＵＭ），ＲＺ（ＮＵＭ））を１次元のＲＸ（ＮＵＭ）で代表させて表記し、その他の表記を省略する。ｉ番目のＡｒ原子の初期速度（ｖｉｘ，ｖｉｙ，ｖｉｚ）は、初期温度Ｔ（ＴＥＭＰ）のマクスウェル分布に従うように、ボックス・ミューラー法（例えば、非特許文献７参照。）を用いて次式で表される。
【００５３】
【数５】
2007219831.tif 000007
【００５４】
ここで、ｋＢ（＝６．０２×１０−２３Ｊ）はボルツマン定数、ｍは原子の質量（Ａｒ原子の質量は６．６３×１０−２６ｋｇであり、Ｐｔ原子の質量は３．２４×１０−２４ｋｇ）である。なお、解析シミュレーション処理プログラム内において、Ａｒ原子の速度を格納する変数を、（ｖｉｘ，ｖｉｙ，ｖｉｚ）＝（ＶＸ（ＮＵＭ），ＶＹ（ＮＵＭ），ＶＺ（ＮＵＭ））とした。なお、座標の場合と同様に、特別な必要のない限り、一次元の速度ＶＸで代表させて示すこととする。
【００５５】
Ｐｔ原子は最近接原子間距離ｒ０（＝Ｒ０）を０．２７７ｎｍとして、面心立方格子状に配置した。また、図２に示すように、加熱壁を３層で構成し、その下にファントム原子、さらに下層に固定原子をそれぞれ１層ずつ配置した。Ｐｔ原子に対しては、Ａｒ原子のような格子からの微小移動は特に行っていない。なお、Ｐｔ原子の座標及び速度を代入する変数は、（ＷＲＸ（ＷＮＵＭ），ＷＲＹ（ＷＮＵＭ），ＷＲＺ（ＷＮＵＭ））及び（ＷＶＸ（ＷＮＵＭ），ＷＶＹ（ＷＮＵＭ），ＷＶＺ（ＷＮＵＭ））とし、それぞれＷＲＸ，ＷＶＸで代表して示すこととする。
【００５６】
ハーモニック・ポテンシャル原子における原子間力は、バネの伸縮によって評価されるため、バネの初期長さとバネの先端にどの原子が接続しているかを記憶させることによって、原子間力の演算回数をバネの総数（ＮＵＭＳ）まで低減させることができる。このため、バネの長さと接続原子番号を配列ＲＲ０（１〜ＮＵＭＳ，１〜３）とＫＲＮ（１〜ＮＵＭＳ，１〜２）に記憶させている。
【００５７】
また、Ｐｔ原子よりも下層に配置しているファントム原子は壁面の温度制御を行うために用い、Ｐｔ原子と質量が同一であり、座標（ＷＲＸ）及び速度（ＷＶＸ）を持つ。さらに、ファントム原子よりも下層に配置している固定原子は加熱壁が移動しないために用いているため、質量はＰｔ原子と同一であるものの、速度の値はゼロで、その運動は計算していない。
【００５８】
次いで、境界条件の処理について以下に詳述する。まず、周期境界条件（例えば、非特許文献７参照。）について説明する。
【００５９】
本実施形態では、計算モデルの側面を周期境界条件とし、上面を周期殻境界条件、下面（底面）を加熱壁とする境界条件の組み合わせを提案している（図２参照）。このうち、分子動力学法において最も基本的な境界条件が周期境界条件である。この境界条件は、計算領域内のある原子が側面を通って流出する場合、反対側の側面から同一の原子が流入することから、同様の計算領域が無限に広がる状態をシミュレーションするものである。原子間力の計算では、原子間距離が計算モデルの一辺の長さの半分（以下、一辺の半分長という。）ＶＨＡＬＦ（＝ＶＳＩＺＥ／２）以上離れたか否かを評価し、原子間距離に一辺の半分長ＶＨＡＬＦを加減することによって計算される。例えば、Ｉ番目の原子とＪ番目の原子のｘ方向距離をＸ（＝ＲＸ（Ｉ）−ＲＸ（Ｊ））とすると、ｘ方向距離Ｘが負の一辺の半分長（−ＶＨＡＬＦ）より小さい場合には、ｘ方向距離Ｘの値に一辺の長さＶＳＩＺＥを加え、ｘ方向距離Ｘが一辺の半分長ＶＨＡＬＦより大きい場合には、ｘ方向の距離Ｘから一辺の長さＶＳＩＺＥを減じた値を運動方程式に代入する。さらに、ｙ方向の原子間距離の計算にも同様の計算を行うことにより、計算領域の側面に周期境界条件を適用したこととなる。
【００６０】
なお、加熱原子（ハーモニック・ポテンシャル原子）の場合、バネの両端に接続する原子は常に同じであるため、一辺の長さＶＳＩＺＥを加えるか減ずるか（すなわち、一辺の長さＶＳＩＺＥに乗算する符号データ（１又は−１であり、図１５参照。）を初期計算において配列ＫＲＮ（１〜ＮＵＭＳ，３）と配列ＫＲＮ（１〜ＮＵＭＳ，４）に記憶させることにより、１〜ＮＵＭＳ個のバネ番号を照合するだけで周期境界条件を適用できるようにしている。このような計算は、原子間力の計算処理（図７のステップＳ１４，Ｓ１５、及び図８のステップＳ２８）並びに、Ａｒ原子及びＰｔ原子の並進運動（ステップＳ２５，Ｓ２７）において実行される。
【００６１】
次いで、周期殻境界条件（例えば、非特許文献７参照。）について以下に説明する。計算モデルの上面に適用される周期殻境界条件について述べる。図３は図１の加熱シミュレーション装置１０において実行される周期殻境界条件における原子の処理を示す概念図である。
【００６２】
図３に示すように、計算領域内の上端部に厚さＳＷＺの殻領域を設け、この領域内に存在する原子の位置関係及び速度ベクトルを計算領域外に複写する。複写層を作成することにより、計算領域上面付近の原子間力の計算において、あたかも領域外にも原子が存在するような予測精度を得ることができる。なお、解析シミュレーション例における殻領域の厚さＳＷＺ＝１．２９ｎｍは、原子の相互作用を計算する際のカットオフ距離（後述）と等しくした。計算領域内の原子（例えば原子ｐ）のポテンシャルは、計算領域と複写領域に存在する全ての原子との相互作用から計算される。なお、複写原子は計算領域内に存在する原子に作用する原子間力の精度を向上させることが目的であるため、その運動を次の時間ステップまで追跡する必要はない。よって、複写原子（例えば原子ｑ）に対するポテンシャル計算は行っていない。
【００６３】
計算領域への原子の流入は、殻領域から非殻領域へ移動する原子（原子ａ）が発生した場合、その複写原子（原子ａ’）を殻領域に挿入することでシミュレーションする。また、原子の流出は、原子が殻領域から複写領域へと移動した場合（原子ｂ）、その存在を抹消することでシミュレーションする。なお、それほど多くの頻度で発生するものではないが、例えば、原子ｃが位置ｃ’に移動するような場合、原子ｄを複写した原子ｄ’と過度に重なることがある。このような場合のポテンシャルエネルギーの急上昇を回避するため、本解析シミュレーション処理では、モンテカルロ法（例えば、非特許文献７参照。）によって原子ｃ’の位置を微小移動させている。本解析シミュレーション処理では、複写原子との距離が０．３ｎｍ以下となった場合を過度の重なりとみなし、微小移動の範囲は０．２６ｎｍとしている。このような原子の流入出の判断は、Ａｒ原子の並進後、複写層を作成する前に、毎時間ステップ毎に行っている。また、Ａｒ原子の座標ＲＸ及び速度ＶＸに対する作業メモリ領域をデータメモリ２３内に作成し、流出原子があった場合にも配列の途中でデータが欠落しないよう、Ｉ番目の原子を削除した場合には、Ｉ＋１番目以降の配列値を１つ繰り上げて記憶するような操作を行っている。
【００６４】
次いで、原子間力の計算方法について以下に説明する。本解析シミュレーション処理では、Ａｒ原子間のポテンシャルを次式のＬ−Ｊポテンシャル関数φＬＪを用いてシミュレーションしている。
【００６５】
【数６】
2007219831.tif 000008
【００６６】
ここで、ｒｉｊはｉ番目とｊ番目のＡｒ原子間の距離であり、ｉ番目の原子の位置ベクトルを（ｘｉ，ｙｉ，ｚｉ）とし、ｊ番目の原子の位置ベクトルを（ｘｊ，ｙｊ，ｚｊ）とすると、距離ｒｉｊ（＝ＲＳＱ）は次式で表される。
【００６７】
【数７】
2007219831.tif 000009
【００６８】
εはＬ−Ｊポテンシャル関数の曲線の井戸の深さであり、σはポテンシャルエネルギーが零となる原子間距離、すなわち実質的な原子直径を表す。本実施形態では、Ａｒ原子の一般的な値として、それぞれε（＝ＥＰＳ）＝１．６７×１０−２１（Ｊ）、σ（＝ＳＩＧ）＝３．４０×１０−１０（ｍ）を用いている。
【００６９】
固体壁を構成する原子は振動範囲が非常に小さいため、最近接原子との相互作用のみを考慮した次式のハーモニック・ポテンシャル関数φＨで近似できる。
【００７０】
【数８】
2007219831.tif 000010
【００７１】
ここで、ｋはバネ定数である。Ａｒ原子とＰｔ原子との間のポテンシャルは、Ｌ−Ｊポテンシャル関数（式（３））を用いてシミュレーションされる。ただし、ポテンシャルのパラメータはローレンツ・バーセロ（Ｌｏｒｅｎｔｚ−Ｂｅｒｔｈｅｌｏｔ）則（例えば、非特許文献６参照。）によって与えられ、本実施形態におけるＡｒ−Ｐｔ間のパラメータは、ε（＝ＷＥＰＳ）＝０．８５×１０−２１（Ｊ）、σ（ＷＳＩＧ）＝３．０９×１０−１０（ｍ）を用いた。
【００７２】
各原子間の原子間力は、式（３）及び式（４）で表されるポテンシャルの導関数より導かれる。Ｌ−Ｊ型のポテンシャル関数に対する原子間力ｆＬＪ，ｉｘ（ｒｉｊ），ｆＬＪ，ｉｙ（ｒｉｊ），ｆＬＪ，ｉｚ（ｒｉｊ）（ＦＲＸ，ＦＲＹ，ＦＲＺ）は次式で表される。
【００７３】
【数９】
2007219831.tif 000011
【００７４】
ここで、ｒｉｊ，ｘ，ｒｉｊ，ｙ，ｒｉｊ，ｚはそれぞれ位置ベクトルの成分（ｒｉｊ，ｘ＝ｘｉ−ｘｊ；ｒｉｊ，ｙ＝ｘｉ−ｘｊ；ｒｉｊ，ｚ＝ｘｉ−ｘｊ）である。ハーモニック型のポテンシャルに対する原子間力ｆＨ，ｉｘ，ｆＨ，ｉｙ，ｆＨ，ｉｚは次式で表される。
【００７５】
【数１０】
2007219831.tif 000012
【００７６】
ここで、ｒ０，ｘは初期位置（バネがまったく伸縮しない状態で、配列ＲＲ０（１〜ＮＵＭＳ，１）で表す。）のバネベクトルのｘ方向成分であり、ｒ０，ｙは初期位置（バネがまったく伸縮しない状態で、配列ＲＲ０（１〜ＮＵＭＳ，２）で表す。）のバネベクトルのｙ方向の成分であり、ｒ０，ｚは初期位置（バネがまったく伸縮しない状態で、配列ＲＲ０（１〜ＮＵＭＳ，３）で表す。）のバネベクトルのｚ方向の成分である。本実施形態では、Ａｒ原子同士の原子間力の評価においてカットオフ距離Ｒｃ（ＲＣＵＴ）を適用しており、原子間距離が１．２９ｎｍ以上となる場合は原子間力の計算を行っていない。
【００７７】
なお、後述する運動方程式に代入するべき原子間力は、カットオフ距離内に存在する全ての原子からの力の総和となる。この値を計算するため、解析シミュレーション処理プログラム内で原子間力を記憶する配列を原子１つにつき１つ用意し、ｉとｊの組み合わせを変えて、原子間力（ＦＲＸ，ＦＲＹ，ＦＲＺ）を積分している。これらの原子間力ベクトルの和の配列は、Ａｒ原子が（ＦＸ（ＮＵＭ），ＦＹ（ＮＵＭ），ＦＺ（ＮＵＭ））、Ｐｔ原子が（ＷＦＸ（ＷＮＵＭ），ＷＦＹ（ＷＮＵＭ），ＷＦＺ（ＷＮＵＭ））となる。
【００７８】
次いで、原子の並進及び速度の計算処理について以下に説明する。分子動力学法では、原子を古典力学のニュートンの運動方程式に従う質点であるとみなし、その運動を評価する。これより、原子ｉの運動は、時間ｔに関する次式の微分方程式で表される。
【００７９】
【数１１】
2007219831.tif 000013
【００８０】
ここで、ベクトルＦｉ＝（Ｆｉｘ，Ｆｉｙ，Ｆｉｚ）は式（５）又は式（６）で表される各原子からの原子間力ベクトルの和である。
【００８１】
式（７）は、ベロシティ・ベルレット法（Ｖｅｌｏｃｉｔｙ  Ｖｅｒｌｅｔ法）（例えば、非特許文献７参照。）により次のように離散化され、Ａｒ原子及びＰｔ原子のｎ＋１回目のステップの座標は、ｎ回目のステップの座標、速度及び原子間力から次式のように表される。
【００８２】
【数１２】
2007219831.tif 000014
【００８３】
ここで、上付きのｎは時間ステップを表す。また、ｎ＋１回目のステップの流速は、ｎ及びｎ＋１回目のステップの原子間力を用いて次式のように与えられる。
【００８４】
【数１３】
2007219831.tif 000015
【００８５】
式（９）を計算するため、本実施形態では、１時間ステップ前の原子間力を記憶する配列（ＦＸ２（ＮＵＭ），ＦＹ２（ＮＵＭ），ＦＺ２（ＮＵＭ））及び（ＷＦＸ２（ＷＮＵＭ），ＷＦＹ２（ＷＮＵＭ），ＷＦＺ２（ＷＮＵＭ））を用意し、並進運動の後、これらの配列値を更新するようにしている。
【００８６】
さらに、マルチ・タイムステップ法について以下に説明する。差分化による誤差には局所誤差と累積誤差の二種類がある。局所誤差は１時間ステップの計算過程で生じる差分化に伴う誤差であり、時間ステップ間隔Δｔが小さくなるほど小さくなる。一方、累積誤差はこの局所誤差が全積分区間で累積されたもので、全時間ステップ数（１／Δｔに比例する。）が大きいほどこの誤差は増加する。よって、時間ステップ間隔Δｔは小さければよいというものではなく、ポテンシャル関数から評価される時間スケールτに対して、差分誤差が出ない程度のオーダーに設定する必要がある。
【００８７】
Ｌ−Ｊポテンシャル関数のように２原子間の距離ｒｉｊに対してポテンシャルが（ｒｉｊ／σＡｒ）の関数で表される場合、時間スケールτＡｒは次式で表される。
【００８８】
【数１４】
2007219831.tif 000016
【００８９】
この値は、本解析シミュレーション条件では、時間スケールτＡｒ＝２．１４×１０−１２ｓｅｃとなる。一方、ハーモニック・ポテンシャルの場合、時間スケールτＰｔは次式で表される。
【００９０】
【数１５】
2007219831.tif 000017
【００９１】
この値は、本解析シミュレーション条件を代入すると、τＰｔ＝６．３×１０−１３ｓｅｃとなり、Ｌ−Ｊポテンシャルの半分以下であることがわかる。よって、Ｌ−Ｊ系のポテンシャルを持つＡｒ原子の時間ステップ間隔と、ハーモニック系のポテンシャルを持つＰｔ原子の時間ステップ間隔を同一とすることは、差分化の誤差に配慮する上で好ましくない。よって、Ｌ−Ｊ系ポテンシャル原子とハーモニック・ポテンシャル原子が混在した解析シミュレーションを行う場合には、マルチ・タイムステップ法により、ハーモニック・ポテンシャル原子の作用を２倍とすることが行われている。本実施形態では、周期殻境界条件を適用した解析シミュレーションで初めてＡｒ原子の作用とＰｔ原子の作用を異なる時間ステップ間隔で計算するマルチ・タイムステップ法を用いている。具体的は、Ａｒ原子が１時間ステップ進行する間に（すなわち、図８のステップＳ２５からステップＳ２９の間）、Ｐｔ原子が２時間ステップ進行するようなアルゴリズムとなっている。なお、２時間ステップ進行する間のＰｔ原子の座標及び速度を記憶するため、マルチ・タイムステップ用の配列ＷＲＸ０５，ＷＶＸ０５，ＷＦＸ０５を用意し、ＷＲＸ，ＷＶＸＷＦＸと交互に値を更新している。
【００９２】
次いで、固体壁の温度制御について以下に説明する。固体壁の温度制御は、ファントム原子にランダムな力を与えることによって行う。ファントム原子の運動方程式は、次式のランジュバン方程式（例えば、非特許文献８参照。）により記述される。
【００９３】
【数１６】
2007219831.tif 000018
【００９４】
ここで、ｍＰｈはファントム原子の質量（Ｐｔ原子と同じ質量を有する。）であり、Ｇｉｘ，Ｇｉｙ，Ｇｉｚは外力として与えるランダム力である。α（ＡＬＰ）はファントム原子の減衰係数であり、デバイ周波数ωＤ（角周波数の単位を有する。）を用いて次式で与えられる。
【００９５】
【数１７】
2007219831.tif 000019
【００９６】
Ｐｔ原子の場合、この値はα＝５．１８×１０−１２ｋｇ／ｓｅｃとなる。固体壁温度を一定に保つ計算では、ランダム力によって与えられるエネルギーの期待値が、壁面温度ＴＷにおける原子粘性による減衰で奪われるエネルギーに等しくなるようにする。これより、ランダム力は標準偏差σＦが次式となるとなる正規分布に従うように与えられる。
【００９７】
【数１８】
2007219831.tif 000020
【００９８】
本実施形態では、毎回のランダム力を与える間の時間を温度一定の状態とみなして、Ｐｔ原子の初期温度（ＷＴＥＭＰ）を式（１４）の壁面温度ＴＷに代入し、ランダム力を与える際の標準偏差σＦを求めている。さらに、加熱を開始して壁面温度ＴＷを上昇させる場合には、ランダム力を与える時間ステップ間隔（ＮＦＯＣ）を０．４ｐｓｅｃ（時間ステップ間隔ＮＦＯＣ１＝２０００）から０．１ｐｓｅｃ（時間ステップ間隔ＮＦＯＣ２＝５００）へと短く設定している。なお、ＮＦＯＣ１やＮＦＯＣ２の値は事前に予備計算を行い、系の温度を観察することにより得られた値である。また、本実施形態では、壁面温度ＴＷが１時間ステップ当たり０．００１Ｋずつ上昇するようにした。なお、ランダム力はＰｔ原子の並進計算時において、ファントム原子の原子間力（ＷＦＸ）に与えられる。なお、現状では、６ｐｓｅｃの間に固体壁の平均温度を１００Ｋ程度上昇させる解析シミュレーションを行っているが、時間ステップ間隔ＮＦＯＣを変更することにより、よりゆるやかな温度上昇率で流体を加熱する解析シミュレーションも実施できる。
【００９９】
次いで、温度の計算について以下に説明する（例えば、非特許文献７のｐｐ．６参照。）系内の原子の速度がマクスウェル分布に従ったものと仮定すると、時間ステップ間隔毎の温度Ｔ’とＮ個の原子が存在する系の運動エネルギーとの間には次のような関係が成り立つ。
【０１００】
【数１９】
2007219831.tif 000021
【０１０１】
ここで、ηは自由度を表す定数であり、３次元シミュレーションの場合にはη＝３である。本実施形態では、任意の時間ステップ間隔で、式（１５）で計算されるＡｒ原子及びＰｔ原子の平均温度Ｔ’Ａｒ，Ｔ’Ｐｔを計算して出力する。また、全原子の座標及び速度の値、さらに必要に応じて、ポテンシャルエネルギーの総和などを出力する。
【０１０２】
さらに、速度スケーリング法（例えば、非特許文献７参照。）について以下に説明する。分子動力学法による解析シミュレーションでは、マクスウェル分布に従う初期速度を与えたとしても、系の温度が解析シミュレーション条件に一致しないことがある。このため、ある時間ステップ間隔における系の平均温度が解析シミュレーション条件に一致するように速度スケーリング法を用いた計算を行わなければならない。以下にこの方法について説明する。
【０１０３】
この方法は、毎時間ステップ（時間ステップ間隔）、又は数時間ステップ毎に、設定温度になるように、計算領域内に存在する原子の速度にスケーリングさせたい倍率βをかけて、系の運動エネルギーを一定にするものである。しかし、速度を人為的に変更するため、この方法を行っている間には物理量の統計平均を取ることはできない。系の設定温度（初期温度）をＴとすると、倍率β（＝Ｃ１）は次式で表される。
【０１０４】
【数２０】
2007219831.tif 000022
【０１０５】
本実施形態では、定常計算においてＡｒ原子に対して毎時間ステップに速度スケーリングをかけ、ＮＶＳＣで与えられる時間ステップ数まで定常計算を行っている。なお、加熱が開始された後は速度スケーリングをかけていない。
【０１０６】
さらに、以上の原理に基づいた本実施形態に係る図６乃至図１０の加熱シミュレーション処理の詳細について以下に説明する。
【０１０７】
図６は図１の加熱シミュレーション装置によって実行される加熱シミュレーション処理を示すフローチャートである。図６において、まず、ステップＳ１において、図７の初期計算処理を実行し、ステップＳ２において図８の主計算処理を実行して終了する。なお、以下の処理では、説明の簡単化のために、ｘ座標のみの処理のみを記載する。
【０１０８】
図７は図６のサブルーチンである初期計算処理（ステップＳ１）を示すフローチャートである。なお、変数の記号、内容、値、単位については図１５に示す。
【０１０９】
まず、ステップＳ１１において計算条件の設定処理において変数ＮＲＥＰＴ、ＮＶＳＣ、ＤＴＩＭＥ、ＲＣＵＴ、ＴＥＭＰ、ＶＳＩＺＥ、ＺＶＳＩＺＥ、ＮＵＭ、ＷＮＵＭ、ＮＦＯＣ１、ＮＦＯＣ２を入力した後、ステップＳ１２において以下の原子の初期座標及び初期速度の代入処理を実行する。
（ａ）式（１）を用いて、間隔Ａの体心立方格子の座標（初期座標）と乱数に基づいてＡｒ原子の座標ＲＸ（Ｉ）（Ｉ＝１〜ＮＵＭ）を計算する。
（ｂ）最近接原子間距離がＲ０となる面心立方格子の座標をＷＲＸ（Ｉ）（Ｉ＝１〜ＷＮＵＭ）に代入する。
（ｃ）Ｐｔ原子とＰｔ原子とを結ぶ全ての初期バネ長さをＲＲ０（１〜ＮＵＭＳ，１〜３）に記憶する。なお、バネの長さと、バネに接続されているＰｔ原子の番号は、変数配列ＫＲＮ（１〜ＮＵＭＳ，１）と変数配列ＫＲＮ（１〜ＮＵＭＳ，２）に記憶されている。詳細は図１５参照。
（ｄ）式（２）を用いて、温度ＴＥＭＰ、原子質量ＭＡＳＳ、乱数ＲＡに基づいて、Ａｒ原子の速度ＶＸ（Ｉ）（Ｉ＝１〜ＮＵＭ）を計算する。
（ｅ）式（２）を用いて、温度ＷＴＥＭＰ、原子質量ＷＭＡＳＳ、乱数ＲＡに基づいて、Ｐｔ原子の速度ＷＶＸ（Ｉ）（Ｉ＝１〜ＷＮＵＭ）を計算する。
（ｆ）次いで、次式の処理を実行する。
［数１］
ＶＸ２←ＶＸ
［数２］
ＷＲＸ０５←ＷＲＸ
【０１１０】
なお、当該明細書において、数式がイメージ入力された墨付き括弧の数番号と、数式が文字入力された大括弧の数式番号とを混在して用いており、また、当該明細書において重要な式の数式番号として「式（１）」の形式を用いて数式番号を式の最後部に付与して（付与していない数式も存在する）用いることとする。また、←は右辺の変数値又は変数値の計算値を左辺の変数に代入することを示す。さらに、＊は乗算を示す。
【０１１１】
次いで、ステップＳ１３において複写層の生成処理（上面の周期殻境界条件）を以下のごとく実行する。
【０１１２】
（ａ）次式の処理を実行する。
［数３］
ＩＤＵＭ←ＮＵＭ
（ｂ）複写１回につき次式の処理を実行する。
［数４］
ＩＤＵＭ←ＩＤＵＭ＋１
［数５］
ＲＺ（ＩＤＵＭ）←ＲＺ（ＩＩ）＋ＳＷＺ
ここで、ＩＩは殻領域内に存在する原子の番号である。
（ｃ）複写されるＡｒ原子のｘ座標ＲＸ、ｙ座標ＲＹ、ｚ座標ＲＺ、速度ＶＸ、ＶＹ、ＶＺに基づいて、複写Ａｒ原子の各座標などを次式のようにコピーして計算する。
［数６］
ＲＸ（ＩＤＵＭ）←ＲＸ（ＩＩ）
［数７］
ＲＹ（ＩＤＵＭ）←ＲＹ（ＩＩ）
［数８］
ＶＸ（ＩＤＵＭ）←ＶＸ（ＩＩ）
［数９］
ＶＹ（ＩＤＵＭ）←ＶＹ（ＩＩ）
［数１０］
ＶＺ（ＩＤＵＭ）←ＶＺ（ＩＩ）
【０１１３】
次いで、ステップＳ１４においてＰｔ−Ｐｔ間の原子間力の計算処理において
（ａ）原子間力の積分値を記憶する変数を初期化するために次式の処理を実行する。
［数１１］
ＦＸ←０
［数１２］
ＷＦＸ←０
（ｂ）Ｐｔ原子の座標ＷＲＸとバネの初期長さＲＲ０（Ｋ，１）（Ｋ＝１〜ＮＵＭＳ）に基づいて、次式を用いて、Ｐｔ−Ｐｔ間のバネの伸びＸＳを計算する。
［数１３］
ＸＳ←
｛ＷＲＸ［ＫＲＮ（Ｋ，１）］−ＷＲＸ［ＫＲＮ（Ｋ，２）］｝−ＲＲ０（Ｋ，１）
（ｃ）式（６）を用いて、バネ定数ＫＲ，バネの伸びＸＳに基づいて、次式を用いて、Ｐｔ原子に働く原子間力ＷＦＲＸを計算する。
［数１４］
ＷＦＲＸ←ＫＲ＊ＸＳ
（ｄ）近接原子に対してＷＦＲＸの値を、次式のごとく積分してＰｔ−Ｐｔ間の原子間力ＷＦＸを計算する。
［数１５］
ＷＦＸ［ＫＲＮ（Ｋ，１）］←ＷＦＸ［ＫＲＮ（Ｋ，１）］＋ＷＦＲＸ
［数１６］
ＷＦＸ［ＫＲＮ（Ｋ，２）］←ＷＦＸ［ＫＲＮ（Ｋ，２）］−ＷＦＲＸ
【０１１４】
さらに、ステップＳ１５Ａ，１５Ｂ，１５ＣからなるステップＳ１５においてＡｒ−Ｐｔ及びＡｒ−Ａｒ間の原子間力の計算処理を以下のように実行する。
【０１１５】
ステップＳ１５Ａにおいて以下の処理を実行する。
（ａ）ベクトルの公式を用いて、Ａｒ原子の座標ＲＸ（Ｉ）及びＰｔ原子の座標ＷＲＸ（Ｊ）に基づいて、Ａｒ−Ｐｔ間及びＡｒ−Ａｒ間のｘ方向距離Ｘを計算する。また、同様に、Ａｒ−Ｐｔ間及びＡｒ−Ａｒ間のｙ方向距離Ｙ及びｚ方向距離Ｚを計算する。ここで、ＩはＡｒ原子の番号（Ｉ＝１〜ＩＤＵＭ）、ＪはＰｔ原子の番号（Ｊ＝１〜ＷＮＵＭ）又はＡｒ原子の番号（Ｊ＝１〜ＩＤＵＭ、ただしＩ≠Ｊ）を表す。
（ｂ）ｘ、ｙ、ｚの距離Ｘ、Ｙ、Ｚに基づいて、次式を用いて距離ＲＳＱを計算する。
【０１１６】
【数２１】
2007219831.tif 000023
【０１１７】
ステップＳ１５ＢにおいてＲＳＱ＜ＲＣＵＴであるか否かが判断され、ＹＥＳのときはステップＳ１５Ｃに進む一方、ＮＯのときは元のメインルーチンに戻る。次いで、ステップＳ１５Ｃにおいて以下の処理を実行する。
（ａ）式（５）を用いて、ＷＳＩＧ、ＷＥＰＳ、Ｘに基づいてＡｒ−Ｐｔ間の原子間力ＦＲＸを計算する。
【０１１８】
（ｂ）Ａｒ−Ｐｔ間の個々の原子間力ＦＲＸの値を積分してＡｒ原子及びＰｔ原子の原子間力ベクトルの和ＦＸ及びＷＦＸを計算する。
［数１７］
ＦＸ（Ｉ）←ＦＸ（Ｉ）＋ＦＲＸ
［数１８］
ＷＦＸ（Ｊ）←ＷＦＸ（Ｊ）＋ＦＲＸ
【０１１９】
（ｃ）式（５）を用いて、ＳＩＧ、ＥＰＳ、Ｘに基づいてＡｒ−Ａｒ間に働く原子間力ＦＲＸを計算する。
（ｄ）Ａｒ−Ａｒ間の原子間力ＦＲＸの値を積分してＡｒ原子の原子間力ＦＸを次式を用いて計算する。
［数１９］
ＦＸ（Ｉ）←ＦＸ（Ｉ）＋ＦＲＸ
［数２０］
ＦＸ（Ｊ）←ＦＸ（Ｊ）−ＦＲＸ
【０１２０】
図８は図６のサブルーチンである主計算処理（ステップＳ２）を示すフローチャートである。
【０１２１】
図８において、まず、ステップＳ２１において変数パラメータＩＰＥＰＴに１を代入して初期化する。ステップＳ２２においてステップ計数値パラメータＩＰＥＰＴ＜ＮＶＳＣ（加熱開始ステップ）であるか否かが判断され、ＹＥＳのときは加熱前の定常計算処理なのでステップＳ２３に進む一方、ＮＯのときは加熱開始後の加熱計算処理なのでステップＳ２４に進む。ステップＳ２３においては、定常計算時の処理なので、ランダム力を与えるタイミングの変更なしの処理を実行し、すなわち、次式の処理を行った後、ステップＳ２５に進む。
［数２１］
ＮＦＯＣ←ＮＦＯＣ１（定常計算時の時間ステップ間隔）
一方、ステップＳ２４においては、加熱時の処理なので、ランダム力を与えるタイミングの変更有りの処理を実行し、すなわち、次式の処理を行った後、ステップＳ２５に進む。
［数２２］
ＮＦＯＣ←ＮＦＯＣ２（加熱計算時の時間ステップ間隔）
【０１２２】
ステップＳ２５においてＡｒ原子の並進処理を以下のように行う。すなわち、式（８）を用いて速度ＶＸ、原子間力ＦＸに基づいて、Ａｒ原子の座標ＲＸを初期値としてＡｒ原子の座標ＲＸを計算する。次いで、ステップＳ２６において複写層の生成処理を実行する。すなわち、（ａ）ステップＳ１３と同一の処理を実行し、（ｂ）図９に示す流入出を判別する処理を行う。ここで、Ａｒ原子数ＮＵＭを更新し、前ステップの複写原子を排除し、ＮＵＭの更新に伴ってＶＸ配列も更新する。さらに、ステップＳ２７において、図１０のＰｔ原子の並進処理を実行する。すなわち、式（８）を用いて、速度ＷＶＸ及び原子間力ＷＦＸに基づいて、Ｐｔ原子の座標ＷＲＸを初期値としてＰｔ原子の座標ＷＲＸを計算する。
【０１２３】
次いで、ステップＳ２８において、原子間力の計算処理において、ステップＳ１４、Ｓ１５と同一の処理を実行する。さらに、ステップＳ２９において、以下の処理を実行する。
（ａ）式（９）を用いて、原子間力ＦＸ及びＦＸ２に基づいて、速度ＶＸ２を初期値としてＡｒ原子の速度ＶＸを計算する。
（ｂ）１ステップ前の原子間力ＦＸを記憶するために次式の処理を実行する。
［数２３］
ＦＸ２←ＦＸ
（ｃ）１ステップ前の速度ＶＸを記憶するために次式の処理を実行する。
［数２４］
ＶＸ２←ＶＸ
【０１２４】
次いで、ステップＳ３１においてＩＲＥＰＴ＜ＮＶＳＣ（加熱開始ステップ）か否かが判断され、ＹＥＳのときは加熱前の定常計算処理なのでステップＳ３２に進む一方、ＮＯのときは加熱開始後の加熱計算処理なのでそのままステップＳ３３に進む。ステップＳ３２において、定常計算時の速度スケーリング処理において以下の処理を実行する。
（ａ）式（１６）を用いて、ボルツマン定数ＢＯＬＴＺ、温度ＴＥＭＰ、速度ＶＸに基づいて係数Ｃ１（式の記号は倍率βである。）を計算する。
（ｂ）次式の速度のスケーリング処理を実行する。
［数２５］
ＶＸ←Ｃ１＊ＶＸ
【０１２５】
次いで、ステップＳ３３においてＰｔ原子の並進処理を以下のように実行する。式（９）を用いて、Ｐｔ原子の原子間力ＷＦＸと速度ＷＦ０５に基づいて、速度ＷＶＸ０５を初期値としてＰｔ原子の速度ＷＶＸを計算する。さらに、ステップＳ３４において以下のデータ出力処理を実行する（上述のように、代表してｘ座標のみを示しているが、ｙ座標及びｚ座標に関するデータも出力する。）。
（ａ）座標ＲＸ、ＷＲＸ、速度ＶＸ、Ｐｔ原子の速度ＷＶＸ、原子数ＮＵＭをＣＲＴディスプレイ４３に出力して表示するとともに、プリンタ４４に出力して印字する。
（ｂ）式（１５）を用いて、Ａｒ原子の質量ＳＭＡＳＳ、速度ＶＸ、原子数ＮＵＭ、ボルツマン定数ＢＯＬＴＺに基づいてＡｒ原子の平均温度Ｔ’Ａｒを計算してＣＲＴディスプレイ４３に出力して表示するとともに、プリンタ４４に出力して印字する。
（ｃ）式（１５）を用いて、Ｐｔ原子の質量ＷＭＡＳＳ、速度ＷＶＸ、原子数ＷＮＵＭ、ボルツマン定数ＢＯＬＴＺに基づいて平均温度Ｔ’Ｐｔを計算して出力する。
【０１２６】
さらに、ステップＳ３５においてＩＲＥＰＴ≧ＮＲＥＰＴ（繰り返し回数）であるか否かが判断され、ＹＥＳのときは元のメインルーチンに戻る一方、ＮＯのときはステップＳ３６に進み、ステップＳ３６においてステップ計数値パラメータＩＲＥＰＴが１だけインクリメントされた後、ステップＳ２２に戻る。
【０１２７】
図９は図８のサブルーチンである複写層の生成処理における流入出の判別処理（ステップＳ２６の（ｂ））を示すフローチャートである。図９のＡｒ原子の流入出の判別処理は、ステップＳ４１−Ｓ４２を含む図３におけるＡｒ原子ｂの処理（ステップＳ４０）と、ステップＳ５１−Ｓ５６を含む図３におけるＡｒ原子ａ'の処理（ステップＳ５０）と、計算領域内のＡｒ原子の配列の置き換え処理（ステップＳ５６）とからなる。
【０１２８】
図９において、まず、ステップＳ４１において、各Ａｒ原子の座標の確認し、０．０＜ＲＺ≦ＺＶＳＩＺＥであるか否か、計算領域内の殻領域にあるか否かが判断され、ＹＥＳのときはステップＳ４２に進む一方、ＮＯのときはステップＳ５１に進む。ステップＳ４２において、Ａｒ原子の座標を一時的に次式のように変数配列に代入する。
［数２６］
（ＲＸ２，ＲＹ２，ＲＺ２）←（ＲＸ，ＲＹ，ＲＺ）
【０１２９】
次いで、ステップＳ５１において、殻領域内のＡｒ原子の座標確認を行い、ＲＺ（ＩＩＩ）≦（ＺＶＳＩＺＥ―ＳＷＺ）であるか否かが判断され、すなわちＡｒ原子が殻領域内から非殻領域へ移動したか否かを確認する。ここで、ＩＩＩは１時間ステップ前に殻領域に存在していた原子の番号である。次いで、ステップＳ５２において、次式のように流入Ａｒ原子を追加する。
［数２７］
（ＲＸ２，ＲＹ２，ＲＺ２）←（ＲＸ（ＩＩＩ），ＲＹ（ＩＩＩ），ＲＺ（ＩＩＩ）＋ＳＷＺ）
【０１３０】
さらに、ステップＳ５３において、ステップＳ５２で計算した流入Ａｒ原子の座標とその他のＡｒ原子の座標ＲＸ２に基づいて、その間の距離ＲＳＱを計算し、ステップＳ５４において距離ＲＳＱ≧０．３ｎｍ（所定の範囲内）であるか否かを判断する。ステップＳ５４においてＹＥＳのときはステップＳ５５に進む一方、ＮＯのときはステップＳ５６に進む。ステップＳ５５において、モンテカルロ法を用いてＡｒ原子を０〜０．２６ｎｍの範囲の移動距離を、Ａｒ原子のｘ，ｙ，ｚ座標に加算する。さらに、ステップＳ５６において計算領域内のＡｒ原子の配列を次式のように置き換えた後、元のルーチンに戻る。
［数２８］
（ＲＸ，ＲＹ，ＲＺ）←（ＲＸ２，ＲＹ２，ＲＺ２）
【０１３１】
図１０は図８のサブルーチンであるＰｔ原子の並進処理（ステップ２７）を示すフローチャートである。
【０１３２】
図１０において、まず、ステップＳ６１においてＰｔ原子の並進処理（１回目）を以下のように実行する。
（ａ）式（８）を用いて、ステップＳ２５の処理と同様の計算処理を実行する。すなわち、式（８）を用いて速度ＷＶＸ及び原子間力ＷＦＸに基づいてＰｔ原子の座標ＷＲＸ０５を初期値として、次式のごとくＰｔ原子の座標ＷＲＸを計算する。
［数２９］
ＷＲＸ
←ＷＲＸ０５＋ＷＤＴＩＭＥ＊ＷＶＸ＋ＷＤＴＳＱ＊ＷＦＸ０５/（２＊ＷＭＡＳＳ）
ここで、
［数３０］
ＷＤＴＩＭＥ←ＤＴＩＭＥ＊０．５
［数３１］
ＷＤＴＳＱ←ＷＤＴＩＭＥ＊ＷＤＴＩＭＥ
（ｂ）次いで、Ｐｔ原子の原子間力ＷＦＸ値の保存のために次式を実行する。
［数３２］
ＷＦＸ０５←ＷＦＸ
【０１３３】
次いで、ステップＳ６２においてＰｔ原子のみの力の計算処理（１回目）を実行する。すなわち、原子間力ＷＦＸに関して、ステップＳ１４と同一の計算処理を行う。
【０１３４】
そして、ステップＳ６３において、Ｐｔ原子のみの速度の計算処理（１回目）を実行する。すなわち、式（９）を用いて、Ｐｔ原子の原子間力ＷＦＸと１時間ステップ前の原子間力ＷＦＸ０５に基づいてステップＳ３３の処理と同様にＰｔ原子の速度ＷＶＸ０５を次式のように計算する。
［数３３］
ＷＶＸ０５←ＷＶＸ＋ＷＤＴＩＭＥ＊（ＷＦＸ＋ＷＦＸ０５）／（２＊ＷＭＡＳＳ）
【０１３５】
さらに、ステップＳ６４において、Ｐｔ原子の並進処理（２回目）を１回目の並進処理と同様に実行する。すなわち、式（８）を用いて、Ｐｔ原子の座標ＷＲＸ０５、速度ＷＶＸ０５及び原子間力ＷＦＸに基づいて、Ｐｔ原子の座標ＷＲＸを以下のように計算する。
［数３４］
ＷＦＸ０５←ＷＦＸ
［数３５］
ＷＦＸ←０
［数３６］
ＷＲＸ０５←
ＷＲＸ＋ＷＤＴＩＭＥ＊ＷＶＸ０５＋ＷＤＴＳＱ＊ＷＦＸ０５／（２＊ＷＭＡＳＳ）
【実施例１】
【０１３６】
図１の加熱シミュレーション装置１０を用いた加熱シミュレーションとその結果を以下に説明する。本実施例に係る加熱シミュレーションでは、初期配置の影響を排除するため、計算モデルが温度１００Ｋで平衡状態となるよう、時刻ｔ＝４ｐｓｅｃ（ＮＶＳＣ＝２００００）まで定常計算を行った。Ａｒ原子の時間ステップ間隔はΔｔ＝０．２×１０−３ｐｓｅｃとした。定常計算時には、Ａｒ温度を設定値に保つために毎時間ステップ毎に速度スケーリング法による制御を行った。そして、定常計算の後（ＩＲＥＰＴ＞ＮＶＳＣ）、壁面温度を上昇させ、加熱をシミュレーションした。加熱時には、Ａｒ原子の速度スケーリング処理を行っていない。
【０１３７】
式（１６）を用いて速度ＶＸと原子数ＮＵＭから計算したＡｒ原子の平均温度Ｔ’Ａｒは、Ｐｔ原子の固体壁の平均温度Ｔ’Ｐｔの上昇に伴って増加した。これにより、固体壁の壁面からの熱伝達が達成されていることが確認された。また、Ａｒ原子の温度が飽和温度に達すると、計算領域内に原子密度の低い領域が増加し始め、沸騰による気泡が生成される。これらの結果から本方法による沸騰解析シミュレーションの有効性が確認された。
【０１３８】
さらに、シミュレーション計算結果について以下に説明する。図１１は図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であるＡｒ原子と固体壁の平均温度の時間的変動を示すグラフである。図１１から明らかなように、加熱開始後（時刻ｔ＞４ｐｓｅｃ）、固体壁の平均温度が増加するにつれて、液体Ａｒ原子の温度も増加した。これより、本解析シミュレーションによって壁面からの熱伝達がシミュレーションされていることを確認した。
【０１３９】
図１２は図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であって、時刻ｔ＝１０ｐｓｅｃのときのＡｒ分子及びＰｔ分子の分布を示す斜視図である。また、図１３は図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であって、時刻ｔ＝１０ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図である。さらに、図１４は図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であって、図１４（ａ）は時刻ｔ＝４ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図であり、図１４（ｂ）は時刻ｔ＝６ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図であり、図１４（ｃ）は時刻ｔ＝８ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図であり、図１４（ｄ）は時刻ｔ＝１０ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図である。
【０１４０】
図１２において、時刻ｔ＝１０ｐｓｒｃにおける、計算領域中心断面内の原子分布を示す。ここでは壁面から３ｎｍまでの分布を示すが、壁面近傍にＡｒ原子密度が局所的に減少した領域が発生していることがわかる。原子数密度の減じた領域を可視化するため、計算領域内に０．２ｎｍ間隔の検査格子を配置し、格子点から半径０．４５ｎｍの空間内にＡｒ原子が存在しないような検査格子のみを白丸でプロットした。その結果を図１３に示す。壁面近傍に原子数密度の低い領域の集合が生じており、気泡が生成されていることを示している。
【０１４１】
図１３と同様のプロットを各時間に対して行った結果を図１４に示す。加熱開始時（時刻ｔ＝４ｐｓｅｃ）において、原子数密度の低い領域は計算領域内に散在しており、各領域の体積もそれほど大きくない。しかし、加熱時間が増加するにつれて、壁面近傍における気泡が徐々に結合し、連続した大きな気泡へと成長する。これより、本解析シミュレーション方法によって壁面からの熱伝達により流体が加熱され、沸騰へと至る過程がシミュレーションされることが示された。
【０１４２】
以上説明したように、本実施例においては、計算領域の上端に周期殻境界条件を適用し、下端の固体壁から流体を加熱する分子動力学シミュレーションを行った。壁面の加熱を開始すると、熱伝達によって流体温度が増加し、壁面近傍において沸騰による気泡核が生成・成長する過程が観察された。
【０１４３】
＜本発明の適用範囲＞
さらに、本発明に係る解析シミュレーション方法の適用範囲について以下に説明する。本発明に係る解析シミュレーション方法では、被加熱原子にＬ−Ｊポテンシャル原子である液体Ａｒ、加熱原子にハーモニック・ポテンシャル原子である固体Ｐｔを設定して、加熱シミュレーションを行っている。しかしながら、本発明で提案した境界条件の組み合わせ方法は、原子の種類に依存しないため、被加熱原子や、加熱原子を別の種類の物質とした解析シミュレーションを実施することは可能である。例えば、Ｌ−Ｊポテンシャルは１原子分子からなる物質（希ガスなど）の原子間力を良好にシミュレーションすると言われているが、短距離の反発力と長距離引力を併せ持つため、式（３）及び式（４）におけるパラメータε、σ及び原子の質量ｍの値を変化させることによって、二酸化炭素や水などの運動もシミュレーションすることが可能である。これらの値の例を表１に示す。なお、被加熱原子がＬ−Ｊポテンシャル原子である限り、シミュレーション物質の変更に対して、アルゴリズムを変更する必要はない。
【０１４４】
【表１】
2007219831.tif 000024
【０１４５】
また、ハーモニック・ポテンシャル原子についても、式（６）におけるパラメータであるバネ係数ｋと、質量ｍの値を変化させることにより、アルゴリズムの変更なしに別の原子の運動をシミュレーションすることが可能である。
【０１４６】
次に被加熱物質の相についてであるが、被加熱原子に適用される周期殻境界条件及び周期境界条件のアルゴリズムは原子の相に依存しない。よって、被加熱原子が、Ｌ−Ｊポテンシャル原子である限り、気体、液体、固体のどのような状態でも解析シミュレーションが可能である。なお、対象物質の状態は、密度（すなわち、ユーザによって代入される体積あたりの原子数）によって設定される。ただし、Ｌ−Ｊポテンシャルによって固体をシミュレーションすることはまれである。一方、加熱物質の相については、固体に限定される。これは、ファントム原子に固定原子を連結するようなモデルが、液体及び気体の解析シミュレーションに適用できないからである。よって、本解析シミュレーション手法は、固体の加熱物質から、気体、液体、（固体）の被加熱物質を加熱する解析シミュレーションに応用しうる。
【０１４７】
ポテンシャルの種類については、被加熱原子はＬ−Ｊポテンシャルに限定される。ハーモニック・ポテンシャルでは、結合する相手が変化しない（すなわち、バネは伸縮するが切れることがない）ので、被加熱原子にハーモニック・ポテンシャルを想定した解析シミュレーションを行っても、周期殻境界条件を適用する利点がなくなってしまう。つまり、被加熱原子を、近接原子間の力のみを考慮するハーモニック・ポテンシャルでシミュレーションする場合は、何層か積層した後に固定原子を置いておけば十分である。逆に、加熱原子のポテンシャルは、Ｌ−Ｊ原子でも可能である（すなわち、本発明では、ハーモニック・ポテンシャルに限定する必要もない）。この場合、マルチ・タイムステップは不要となるので、図８のステップＳ２７の処理を、ステップＳ２５の処理に置き換える必要がある。ただし、加熱原子のポテンシャルがＬ−Ｊ型となるのは、固体物質の性質がＬ−Ｊポテンシャルによってシミュレーションされる場合に限られ、このような実例はあまりないようである。
【０１４８】
解析シミュレーションモデルについては、被加熱原子及び加熱原子の初期配置を変更することにより、様々な固体形状から流体原子を加熱する解析シミュレーションが可能である。周期殻境界条件は解析シミュレーションモデルの上部のみならず、複数面又は全面に設定することができる。この場合の複写層の生成では、殻領域厚さ（ＳＷＺ）の加算（又は減算）対象となる原子位置ベクトルの成分（すなわちｘ，ｙ，ｚ方向成分のいずれか）をプログラムにおいて変更する必要がある。また、周期境界条件の適用方向も、周期殻境界条件と重複しないように設定する必要がある。
【０１４９】
なお、周期境界条件の性質から、周期殻境界条件に相対する端面の境界は、周期殻境界条件か又は固体壁が適用されなければならない。このようなプログラムの変更により、例えば、加熱原子を円柱状（図４参照。）や、球状に配置した解析シミュレーションを行うことができ、流体中におかれた微小円柱や微小球などからの加熱シミュレーションを行うことが可能である。図４は本発明の変形例に係る円柱内の加熱解析の一例を示す平面図である。ただし、固定原子及びファントム原子と、被加熱原子との相互作用は考慮していないため、これらの仮想原子が直接、被加熱原子に接するような解析シミュレーションは行えない。
【０１５０】
図５（ａ）は本発明の加熱シミュレーション方法で解析可能な場合の状態を示す斜視図であり、図５（ｂ）は本発明の加熱シミュレーション方法で解析不可能な場合の状態を示す斜視図である。また、上記の内容を表２にまとめる。
【０１５１】
【表２】
2007219831.tif 000025
【０１５２】
最後に、本解析シミュレーションを２次元の解析シミュレーションに適用することは可能である。ただし、座標ベクトルや速度ベクトルの次元を全て２次元のものに変更し、原子間力の計算や運動方程式などの計算も全て２次元成分に対して行うように変更しなくてはならない。これらの変更は可能であるが、原子運動の自由度が実現象と異なるため、解析シミュレーションの予測精度は劣化する。
【０１５３】
＜変形例＞
以上の実施形態又は実施例においては、図６乃至図１０の加熱シミュレーション処理プログラムデータをＣＤ−ＲＯＭ４５ａに格納して実行するときにプログラムメモリ２４にロードして実行しているが、本発明はこれに限らず、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ、ＭＯなどの光ディスク又は光磁気ディスクの記録媒体、もしくは、フロッピーディスクなどの磁気ディスクの記録媒体など種々の記録媒体に格納してもよい。これらの記録媒体は，コンピュータで読み取り可能な記録媒体である。また、図６乃至図１０の加熱シミュレーション処理プログラムのデータを予めプログラムメモリ２４に格納して当該処理を実行してもよい。
【０１５４】
＜本願発明と非特許文献との相違点＞
非特許文献１では、加熱による液体の沸騰現象を解析している。しかし、加熱壁と相対する計算領域の境界に、原子の流入出のない断熱壁（弾性散乱壁）を適用しているため、計算領域をすべて液相で満たすと圧力上昇によって気泡が生じなくなる。そこで、極めて薄い液膜の上に気相を設けた系を扱うこととなり、これにより、非特許文献１の実現象への適用範囲は限定されたものとなっている。これに対して、本発明は、非特許文献１で断熱壁とした境界に周期殻境界条件を適用することにより、流体の流入出を自然にシミュレーションし、連続的に広がった液体を壁面から加熱する状態を解析シミュレーションできるようにした。これにより、一般的な加熱沸騰モデルを、比較的小さな計算領域の分子動力学シミュレーションによって扱えるようになった。
【０１５５】
非特許文献２及び３では、沸騰のシミュレーションを行っているが、固体壁により流体を挟んだ計算モデルを扱っている。このような対称性のよい計算モデルは比較的容易に生成できるが、得られた結果は、液相の対流などの点において加熱沸騰とは異なるものと考えられる。よって、非特許文献２及び３の計算は、極めて小さな系の減圧沸騰現象をシミュレーションしたものとなっている。これに対して、本発明は、非特許文献２及び３で扱われている固体壁を加熱面として用いているが、それと相対する境界に周期殻境界条件を用いた点で非特許文献２及び３と相違する。また、固体壁の取り扱いにおいても、一定温度ではなく、時間経過とともに壁面温度を上昇させている。
【０１５６】
非特許文献４及び５では、計算領域の内外に比較的広範に広がる流体の取り扱いに対して、新しい境界条件を提案している。しかし、その適用例は円柱周り流に限られており、片側から流体を加熱するような非対称な計算領域への適用について触れられていない。これに対して、本発明は、加熱される液相に周期殻境界条件を適用し、非特許文献４及び５と同様に計算領域外にも液体が広がる状態をシミュレーションする効果を得ている。しかし、それと相対する面に固体壁を配置し、加熱による沸騰シミュレーションに本境界条件を適用した点で非特許文献４及び５と相違する。
【０１５７】
非特許文献６では、特に沸騰現象を扱っているわけではないが、固体壁に相対する境界に断熱壁を用いた解析シミュレーション例として取り上げた。計算モデルとしては、非特許文献１によく似ている。本発明との直接的な関連性は小さい。
【産業上の利用可能性】
【０１５８】
以上詳述したように、本発明に係る加熱シミュレーション装置及び方法によれば、固体壁からの熱負荷によってもたらされる被加熱原子の加熱現象を、従来技術に比較して小規模な計算モデルで計算し、計算コストを大幅に削減できるとともに、従来技術に比較して高精度で加熱現象をシミュレーションでき、相変化冷却などのシミュレーションを行うことが可能となる。
【０１５９】
半導体部品の高密度実装に伴い、半導体パッケージやプリント基板を冷却する手法として、液体冷却や沸騰を利用した相変化冷却が利用されている。数年後には、電子部品からの発熱量は数十Ｗ／ｃｍ２に達すると予想されており、このような高熱負荷に対応するため、マイクロオーダーの流路に液体を流動させる冷却システムや、ヒートシンク内部で液体を沸騰させる冷却システムの製品化が進められている。また、フロリナートやノベックなどのフッ素系液体を半導体パッケージに直接触れさせ、沸騰させることによって除熱する方式も検討されている。これらの産業では、微小流路内の液体沸騰や、特殊な濡れ性を持つ液体の沸騰現象をシミュレーションする方法が望まれており、本発明で提供する分子動力学法の計算方法が利用される可能性が高い。また、このような動向に伴い、分子動力学法による流体シミュレーションの需要は増大する方向にあり、将来的に分子動力学法による汎用熱流体解析コードが販売される場合、広範な沸騰現象を模擬するモデルとして、本発明が組み込まれる可能性がある。なお、原子炉事故時においても、数ＭＷ／ｍ２の高熱負荷によって冷却水が沸騰し、減速材の体積割合が減少することによって自発的な出力抑制が生ずることが指摘されており、このような現象の把握においても、本発明に係る加熱シミュレーション方法が利用される可能性が高い。
【図面の簡単な説明】
【０１６０】
【図１】本発明の一実施形態に係る、分子動力学による高熱負荷液体からの加熱シミュレーション装置１０の構成を示すブロック図である。
【図２】図１の加熱シミュレーション装置１０において用いる計算モデルを示す図である。
【図３】図１の加熱シミュレーション装置１０において実行される周期殻境界条件における原子の処理を示す概念図である。
【図４】本発明の変形例に係る円柱内の加熱解析の一例を示す平面図である。
【図５】（ａ）は本発明の加熱シミュレーション方法で解析可能な場合の状態を示す斜視図であり、（ｂ）は本発明の加熱シミュレーション方法で解析不可能な場合の状態を示す斜視図である。
【図６】図１の加熱シミュレーション装置によって実行される加熱シミュレーション処理を示すフローチャートである。
【図７】図６のサブルーチンである初期計算処理（ステップＳ１）を示すフローチャートである。
【図８】図６のサブルーチンである主計算処理（ステップＳ２）を示すフローチャートである。
【図９】図８のサブルーチンである複写層の生成処理における流入出の判別処理（ステップＳ２６の（ｂ））を示すフローチャートである。
【図１０】図８のサブルーチンであるＰｔ原子の並進処理（ステップＳ２７）を示すフローチャートである。
【図１１】図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であるＡｒ原子と固体壁の平均温度の時間的変動を示すグラフである。
【図１２】図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であって、時刻ｔ＝１０ｐｓｅｃのときのＡｒ分子及びＰｔ分子の分布を示す斜視図である。
【図１３】図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であって、時刻ｔ＝１０ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図である。
【図１４】図１の加熱シミュレーション装置１０によるシミュレーション結果の一例であって、（ａ）は時刻ｔ＝４ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図であり、（ｂ）は時刻ｔ＝６ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図であり、（ｃ）は時刻ｔ＝８ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図であり、（ｄ）は時刻ｔ＝１０ｐｓｅｃのときの壁面近傍領域の気泡分布を示す斜視図である。
【図１５】図１の加熱シミュレーション装置１０において用いる各変数の内容、値、単位及び式中の記号を示す図である。
【符号の説明】
【０１６１】
１０…加熱シミュレーション装置、
２０…ＣＰＵ、
２１…ＲＯＭ、
２２…ＲＡＭ、
２３…データメモリ、
２４…プログラムメモリ、
３０…バス、
３１…キーボードインターフェース、
３２…マウスインターフェース、
３３…ディスプレイインターフェース、
３４…プリンタインターフェース、
３５…ドライブ装置インターフェース、
４１…キーボード、
４２…マウス、
４３…ＣＲＴディスプレイ、
４４…プリンタ、
４５…ＣＤ−ＲＯＭドライブ装置、
４５ａ…ＣＤ−ＲＯＭ。
【図１】
2007219831.tif 000026
【図３】
2007219831.tif 000027
【図４】
2007219831.tif 000028
【図５】
2007219831.tif 000029
【図６】
2007219831.tif 000030
【図７】
2007219831.tif 000031
【図８】
2007219831.tif 000032
【図９】
2007219831.tif 000033
【図１０】
2007219831.tif 000034
【図１１】
2007219831.tif 000035
【図１５】
2007219831.tif 000036
【図２】
2007219831.tif 000037
【図１２】
2007219831.tif 000038
【図１３】
2007219831.tif 000039
【図１４】
2007219831.tif 000040