(57)【要約】      （修正有）【課題】流れがある空間における粒子成長過程を短時間に精度良く模擬することが可能なシミュレーション方法等を提供する。【解決手段】制御部は、予め、空間における流体の流速分布及び粒子の濃度分布を計算し、空間内を流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により流体の流速分布に基づいて微小領域の軌跡を計算する。次に、燃焼計算の結果データを入力し（ステップ１１）、粒子の濃度分布に基づいて流線上における粒子の生成位置を示す時間の関数を決定し（ステップ１２）、流線上における粒子の生成位置を示す時間の関数に基づいて、粒子成長計算を行う（ステップ１３）。終了条件を満たす場合（ステップ１５のＹｅｓ）、結果データファイルを出力し（ステップ１６）、処理を終了する。【選択図】図４
2014035614.tif 000002
【特許請求の範囲】
【請求項１】
流れがある空間における粒子成長過程のシミュレーション方法であって、
少なくとも、前記空間における流体の流速分布、温度分布及び粒子の濃度分布を計算する流体計算ステップと、
前記空間内を前記流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により前記流体の流速分布に基づいて前記微小領域の軌跡を計算する軌跡計算ステップと、
少なくとも、前記温度分布に基づいて前記微小領域の軌跡上における温度を示す時間の関数及び前記粒子の濃度分布に基づいて前記微小領域の軌跡上における粒子の生成位置を示す時間の関数を決定し、少なくとも前記微小領域の軌跡上における温度を示す時間の関数及び前記微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算ステップと、
を含むことを特徴とするシミュレーション方法。
【請求項２】
前記粒子計算ステップは、計算対象の前記微小領域に含まれる粒子の密度を、計算対象の前記微小領域の軌跡と隣接する前記微小領域の軌跡との間隔、及び計算対象の前記微小領域の移動速度を用いて補正する
ことを特徴とする請求項１に記載のシミュレーション方法。
【請求項３】
前記粒子計算ステップは、
前記微小領域の軌跡に垂直な面上、又は前記流れの主流方向に垂直な面上であって、計算対象の前記微小領域の軌跡と隣接する前記微小領域の軌跡との間隔から定まる所定の範囲の面積を算出し、
前記所定の範囲の面積と計算対象の前記微小領域の移動速度の積を補正係数とし、
計算対象の前記微小領域に含まれる粒子の実数を前記補正係数で割った値を、計算対象の前記微小領域に含まれる粒子の密度とする
ことを特徴とする請求項２に記載のシミュレーション方法。
【請求項４】
前記流体計算ステップは、第１計算時間刻みによって定義される計算時刻ごとに計算を行い、
前記粒子計算ステップでは、前記第１計算時刻刻みと異なる第２計算時間刻みによって定義される計算時刻ごとに計算を行う
ことで、それぞれ前記流体計算ステップの計算負荷と前記粒子計算ステップでの計算負荷を独立に設定することで、計算負荷を低減することを特徴とする請求項１乃至請求項３のいずれかに記載のシミュレーション方法。
【請求項５】
コンピュータに、
流れがある空間における粒子成長過程のシミュレーション方法の各ステップとして、
少なくとも、前記空間における流体の流速分布、温度分布及び粒子の濃度分布を計算する流体計算ステップと、
前記空間内を前記流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により前記流体の流速分布に基づいて前記微小領域の軌跡を計算する軌跡計算ステップと、
少なくとも、前記温度分布に基づいて前記微小領域の軌跡上における温度を示す時間の関数及び前記粒子の濃度分布に基づいて前記微小領域の軌跡上における粒子の生成位置を示す時間の関数を決定し、少なくとも前記微小領域の軌跡上における温度を示す時間の関数及び前記微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算ステップと、
を実行させるためのプログラム。
【請求項６】
流れがある空間における粒子成長過程のシミュレーション装置であって、
少なくとも、前記空間における流体の流速分布、温度分布及び粒子の濃度分布を計算する流体計算手段と、
前記空間内を前記流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により前記流体の流速分布に基づいて前記微小領域の軌跡を計算する軌跡計算手段と、
少なくとも、前記温度分布に基づいて前記微小領域の軌跡上における温度を示す時間の関数及び前記粒子の濃度分布に基づいて前記微小領域の軌跡上における粒子の生成位置を示す時間の関数を決定し、少なくとも前記微小領域の軌跡上における温度を示す時間の関数及び前記微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算手段と、
を備えることを特徴とするシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流れがある空間における粒子成長過程のラグランジェ記述を用いたシミュレーション方法等に関するものである。
【背景技術】
【０００２】
光ファイバの代表的な製造方法として、ＶＡＤ法（Vapor phase
axial deposition method、気相軸付け法）が知られている。ＶＡＤ法は、酸素−水素バーナーの火炎中に石英の原料ガス（四塩化珪素）を導入してシリカ粒子を生成し、その集合体を形成する方法である。
【０００３】
また、リチウムイオン電池の正極活物質として、リン酸鉄リチウムやケイ酸鉄リチウム等が注目されている。これらの材料はそのままでは導電性が低く、結晶内部のリチウムイオン拡散速度が遅いため、１０〜１００ｎｍの一次粒子径を有する微粒子にする必要がある。このようなナノサイズの微粒子の製造方法として原料液をバーナーに導入し、バーナーの先端に取付けたノズルから噴霧して火炎中で燃焼させることによって、微粒子を製造する方法がある。このような粒子は衝突を繰り返すことで、粒子の塊（凝集粒子）として成長していく。また凝集粒子を構成する個々の粒子（一次粒子）は、高温下では焼結により周囲の一次粒子と融合し、より大きな一次粒子に成長することもある。
【０００４】
このように、粒子の成長は複雑な過程である。このような現象を精度よく計算するには、衝突頻度を正確に計算する必要がある。衝突頻度は粒子の大きさと粒子の周囲の気体の平均自由行程の比によって、式が異なる。粒子の大きさが平均自由行程より小さい領域は自由分子領域と呼ばれ、衝突頻度は温度の関数である。一方、粒子の大きさが平均自由行程より大きい領域は連続領域と呼ばれ、衝突頻度の計算には温度と平均自由行程が必要となる。粒子の径が平均自由行程と同じくらいの領域は遷移領域と呼ばれ、この領域でも衝突頻度の計算には温度と平均自由行程が必要となる。
【０００５】
ＶＡＤなどの火炎中における粒子成長過程を模擬するためには、粒子成長計算と燃焼計算を連成させる必要がある。連成方法としては、オイラー記述に基づく方法とラグランジュ記述に基づく方法の２つが知られている。
【０００６】
例えば、川の流れの流速を計測する場合を考える。オイラー記述に基づく方法は、川の様々な場所に速度センサを固定しておき、複数個の速度センサの値をモニタリングする方法と捉えることができる。一方、ラグランジュ記述に基づく方法は、川に速度センサ付きの「浮き」を浮かべて、流れに沿って移動する「浮き」に付けられた速度センサの値をモニタリングする方法と捉えることができる。両者は、数学的には等価である。火炎中における粒子成長過程の計算を含む一般的な流体の計算では、オイラー記述に基づく方法が用いられており、ラグランジュ記述に基づく方法はあまり用いられない。
【０００７】
特許文献１には、自熔炉（銅精錬プロセスにおいて、銅精鉱を燃焼して銅成分の純度を上げるための設備）内の精鉱粒子同士の衝突及び成長現象を考慮した燃焼現象を予測する方法が開示されている。特許文献１に記載の方法は、まず、ガス相と、精鉱粒子の粒径に従って分類された複数の粒子相とからなるオイラー法の多相流モデルを作成する。次に、オイラー法の多相流モデルを解いて燃焼解析を行い、それぞれの粒子相の少なくとも体積分率を計算する。次に、複数の粒子相のうちの２つの粒子相で構成される粒子相の組み合わせの全てについて、それぞれの組み合わせ毎に、当該組み合わせ内のそれぞれの粒子相の燃焼解析の出力を利用して、粒径の大きい方の粒子相であるα相と粒径の小さい方の粒子相であるβ相との間の衝突確率ηαを求め、β相の燃焼解析によって計算された体積分率の一部又は全部を、求められた衝突確率ηαに従って、α相の燃焼解析によって計算された体積分率に加える。
【０００８】
また、非特許文献１には、閉じた空間の中での粒子の成長過程の計算方法が開示されている。非特許文献１に記載の計算方法では、粒子が多数衝突してできた凝集粒子の形状を考慮するために、体積と表面積をいくつかの階級に分けて、階級間の総体積の変化で粒子の成長を計算する。ここでの粒子の総体積は、その階級を代表する粒子の半径から計算される粒子の体積にその階級に属する粒子の数を乗じたもので定義される。また、非特許文献２には、オイラー記述に基づく方法を用いて、非特許文献１記載の粒子の成長計算を火炎中で行う計算方法が開示されている。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開２００６−７７２８８号公報
【非特許文献】
【００１０】
【非特許文献１】Yun Xiong and Sotiris E.Pratsinis, "Formation of agglomerate particles by coagulation andsintering-Part I. A two-dimensional solution of the population balance equation",Journal of Aerosol Science, 24, p283-, 1993
【００１１】
【非特許文献２】B.W. Lee, and J.I. Jeong, and J.Y.Hwang, and M. Choi, and S.H. Chung, "Analysis of growth of non-sphericalsilica particles in a counterflow diffusion flame considering chemicalreactions, coagulation and coalescence", Aerosol Science 32, p165-, 2001
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、特許文献１、非特許文献２を含めて、オイラー記述に基づく方法では、粒子の大きさの階級の数と同じ数の方程式を解く必要があり、計算負荷が高いという問題がある。例えば、非特許文献１に従って粒子の体積を３５の階級に分けた場合、表面積の階級も考慮すると、体積と表面積で６３０の階級ができる。この６３０個の階級について、それぞれ方程式を解く必要があるため、非常に計算負荷が高く、膨大な計算時間を要する。
【００１３】
また、ＶＡＤ法や噴霧燃焼法などを模擬するための系では、火炎の温度や流れの計算を必要とする領域が、粒子の成長の計算に関与する領域よりも広いことが一般的である。一方、粒子の成長、特に反応直後の粒子の成長は、火炎の温度や流れの変化と比較して非常に速いため、計算時間刻みを小さくする必要がある。従って、火炎の温度や流れの計算と粒子の成長の計算を同時に行う場合、精度良く計算するためには、計算領域を広くするとともに、計算時間刻みを小さくする必要がある。計算領域を広くすることと、計算時間刻みを小さくすることは、共に計算負荷を増大させるため、オイラー記述に基づく方法では、精度良く計算しようとすると、計算負荷が飛躍的に増大してしまうことになる。
【００１４】
本発明は、前述した問題点に鑑みてなされたものであり、その目的とすることは、流れがある空間における粒子成長過程を短時間に精度良く模擬することが可能なシミュレーション方法等を提供することである。
【課題を解決するための手段】
【００１５】
前述した目的を達成するための第１の発明は、流れがある空間における粒子成長過程のシミュレーション方法であって、少なくとも、前記空間における流体の流速分布、温度分布及び粒子の濃度分布を計算する流体計算ステップと、前記空間内を前記流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により前記流体の流速分布に基づいて前記微小領域の軌跡を計算する軌跡計算ステップと、少なくとも、前記温度分布に基づいて前記微小領域の軌跡上における温度を示す時間の関数及び前記粒子の濃度分布に基づいて前記微小領域の軌跡上における粒子の生成位置を示す時間の関数を決定し、少なくとも前記微小領域の軌跡上における温度を示す時間の関数及び前記微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算ステップと、を含むことを特徴とするシミュレーション方法である。第１の発明によって、流れがある空間における粒子成長過程を短時間に精度良く模擬することができる。
【００１６】
第１の発明における前記粒子計算ステップは、計算対象の前記微小領域に含まれる粒子の密度を、計算対象の前記微小領域の軌跡と隣接する前記微小領域の軌跡との間隔、及び計算対象の前記微小領域の移動速度を用いて補正するようにしても良い。これによって、流れの影響を受けて、計算対象の微小領域の軌跡と隣接する微小領域の軌跡との間隔、及び計算対象の微小領域の移動速度が大きく変化する場合も、粒子成長過程を精度良く模擬することができる。
【００１７】
また、第１の発明における前記粒子計算ステップは、前記微小領域の軌跡に垂直な面上、又は前記流れの主流方向に垂直な面上であって、計算対象の前記微小領域の軌跡と隣接する前記微小領域の軌跡との間隔から定まる所定の範囲の面積を算出し、前記所定の範囲の面積と計算対象の前記微小領域の移動速度の積を補正係数とし、計算対象の前記微小領域に含まれる粒子の実数を前記補正係数で割った値を、計算対象の前記微小領域に含まれる粒子の密度とするようにしても良い。これによって、精度良く正確な粒子数を算出することができる。また、流れの主流方向に垂直な面で所定の範囲の面積を算出する場合には、計算負荷を少なくすることができる。
【００１８】
前記流体計算ステップは、第１計算時間刻みによって定義される計算時刻ごとに計算を行い、前記粒子計算ステップでは、前記第１計算時刻刻みと異なる第２計算時間刻みによって定義される計算時刻ごとに計算を行うようにしても良い。これによって、計算領域と計算時間刻みを独立に考えることが可能となる。特に、流体計算では、広い計算領域によって精度を確保し、大きい計算時間刻みによって計算負荷を減少させることが可能となる。また、粒子計算では、小さい計算時間刻みによって精度を確保し、微小領域内のみを計算することによって計算負荷を大幅に減少させることが可能となる。このように、流体計算と粒子計算とで、計算時間刻みと計算領域の設定を、それぞれ独立して適切に設定することで、計算負荷の減少と計算精度の向上を実現することができる。
【００１９】
第２の発明は、コンピュータに、流れがある空間における粒子成長過程のシミュレーション方法の各ステップとして、少なくとも、前記空間における流体の流速分布、温度分布及び粒子の濃度分布を計算する流体計算ステップと、前記空間内を前記流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により前記流体の流速分布に基づいて前記微小領域の軌跡を計算する軌跡計算ステップと、少なくとも、前記温度分布に基づいて前記微小領域の軌跡上における温度を示す時間の関数及び前記粒子の濃度分布に基づいて前記微小領域の軌跡上における粒子の生成位置を示す時間の関数を決定し、少なくとも前記微小領域の軌跡上における温度を示す時間の関数及び前記微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算ステップと、を実行させるためのプログラムである。第２の発明のプログラムを汎用のコンピュータにインストールすることによって、第１の発明のシミュレーション方法を実行することができる。
【００２０】
第３の発明は、流れがある空間における粒子成長過程のシミュレーション装置であって、少なくとも、前記空間における流体の流速分布、温度分布及び粒子の濃度分布を計算する流体計算手段と、前記空間内を前記流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により前記流体の流速分布に基づいて前記微小領域の軌跡を計算する軌跡計算手段と、少なくとも、前記温度分布に基づいて前記微小領域の軌跡上における温度を示す時間の関数及び前記粒子の濃度分布に基づいて前記微小領域の軌跡上における粒子の生成位置を示す時間の関数を決定し、少なくとも前記微小領域の軌跡上における温度を示す時間の関数及び前記微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算手段と、を備えることを特徴とするシミュレーション装置である。第３の発明によって、流れがある空間における粒子成長過程を短時間に精度良く模擬することができる。
【発明の効果】
【００２１】
本発明により、流れがある空間における粒子成長過程を短時間に精度良く模擬することが可能なシミュレーション方法等を提供することができる。
【図面の簡単な説明】
【００２２】
【図１】シミュレーション装置１のハードウエア構成図
【図２】シミュレーション装置１のソフトウエア構成図
【図３】熱流体計算処理の流れを示すフローチャート
【図４】粒子成長計算処理の流れを示すフローチャート
【図５】粒子密度補正処理を説明する図
【図６】シリカ粒子の質量流量の比較図
【発明を実施するための形態】
【００２３】
以下図面に基づいて、本発明の実施の形態を詳細に説明する。最初に、図１、図２を参照しながら、シミュレーション装置１の概要を説明する。
【００２４】
図１は、シミュレーション装置１のハードウエア構成図である。尚、図１のハードウエア構成は一例であり、用途、目的に応じて様々な構成を採ることが可能である。また、シミュレーション装置１は、１台の筐体から構成されても良いし、複数台の筐体から構成されても良い。
【００２５】
シミュレーション装置１は、制御部１１、記憶部１２、メディア入出力部１３、通信制御部１４、入力部１５、表示部１６、周辺機器Ｉ／Ｆ部１７等が、バス１８を介して接続される。
【００２６】
制御部１１は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等で構成される。
【００２７】
ＣＰＵは、記憶部１２、ＲＯＭ、記録媒体等に格納されるプログラムをＲＡＭ上のワークメモリ領域に呼び出して実行し、バス１８を介して接続された各装置を駆動制御し、シミュレーション装置１が行う後述する処理を実現する。ＲＯＭは、不揮発性メモリであり、コンピュータのブートプログラムやＢＩＯＳ等のプログラム、データ等を恒久的に保持している。ＲＡＭは、揮発性メモリであり、記憶部１２、ＲＯＭ、記録媒体等からロードしたプログラム、データ等を一時的に保持するとともに、制御部１１が各種処理を行うために使用するワークエリアを備える。
【００２８】
記憶部１２は、ＨＤＤ（ハードディスクドライブ）であり、制御部１１が実行するプログラム、プログラム実行に必要なデータ、ＯＳ（オペレーティングシステム）等が格納される。プログラムに関しては、ＯＳ（オペレーティングシステム）に相当する制御プログラムや、後述する処理をコンピュータに実行させるためのアプリケーションプログラムが格納されている。これらの各プログラムコードは、制御部１１により必要に応じて読み出されてＲＡＭに移され、ＣＰＵに読み出されて各種の手段として実行される。
【００２９】
メディア入出力部１３（ドライブ装置）は、データの入出力を行い、例えば、ＣＤドライブ（−ＲＯＭ、−Ｒ、−ＲＷ等）、ＤＶＤドライブ（−ＲＯＭ、−Ｒ、−ＲＷ等）等のメディア入出力装置を有する。通信制御部１４は、通信制御装置、通信ポート等を有し、コンピュータとネットワーク間の通信を媒介する通信インタフェースであり、ネットワークを介して、他のコンピュータ間との通信制御を行う。ネットワークは、有線、無線を問わない。
【００３０】
入力部１５は、データの入力を行い、例えば、キーボード、マウス等のポインティングデバイス、テンキー等の入力装置を有する。入力部１５を介して、コンピュータに対して、操作指示、動作指示、データ入力等を行うことができる。表示部１６は、液晶パネル等のディスプレイ装置、ディスプレイ装置と連携してコンピュータのビデオ機能を実現するための論理回路等（ビデオアダプタ等）を有する。
【００３１】
周辺機器Ｉ／Ｆ（インタフェース）部１７は、コンピュータに周辺機器を接続させるためのポートであり、周辺機器Ｉ／Ｆ部１７を介してコンピュータは周辺機器とのデータの送受信を行う。周辺機器Ｉ／Ｆ部１７は、ＵＳＢやＩＥＥＥ１３９４やＲＳ−２３２Ｃ等で構成されており、通常複数の周辺機器Ｉ／Ｆを有する。周辺機器との接続形態は有線、無線を問わない。バス１８は、各装置間の制御信号、データ信号等の授受を媒介する経路である。
【００３２】
図２は、シミュレーション装置１のソフトウエア構成図である。シミュレーション装置１は、例えば、市販のシミュレーションソフトがインストールされることによって、プリプロセッサ２１、ソルバー２２、ポストプロセッサ２３等の機能を備える。
【００３３】
プリプロセッサ２１は、後述するシミュレーション処理を行うための前処理を行い、条件データ２４を生成する機能である。条件データ２４は、シミュレーション処理に必要な計算条件、初期値等である。
【００３４】
ソルバー２２は、条件データ２４に基づいて、シミュレーションを行い、結果データ２５を生成する機能である。ソルバー２２は、時間ステップごとに、解析対象に応じて定式化された方程式を解く。方程式の定式化とソルバー２２の処理の詳細は後述する。
【００３５】
ポストプロセッサ２３は、結果データ２５を表示部１６等に出力する機能である。ポストプロセッサ２３は、変形図、コンター図、ベクトル図、グラフ、アニメーションなど様々な表示を行う。また、ポストプロセッサ２３は、結果データ２５の時間ステップごとの解をリスト形式などで表示する。
【００３６】
次に、図３〜図６を参照しながら、シミュレーション装置１の詳細を説明する。以下では、シミュレーション装置１が実行するシミュレーション処理について、火炎中での粒子成長計算を例に挙げて説明する。
【００３７】
本発明の実施の形態では、従来技術のように燃焼計算と粒子成長計算を同時に行うのではなく、予め燃焼計算を行い、燃焼計算による結果を適用して粒子成長計算を行うことによって、燃焼計算と粒子成長計算を連成させる。このような連成方法を実現するために、本発明の実施の形態では、従来技術のようなオイラー記述に基づく方法ではなく、ラグランジュ記述に基づく方法を採用する。特に、本発明の実施の形態における粒子成長計算では、空間内を流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により微小領域内にある粒子の物理量を計算する。
【００３８】
このように燃焼計算と粒子成長計算を連成させることによって、各計算において計算領域と計算時間刻みを独立に考えることが可能となる。特に、燃焼計算では、広い計算領域によって精度を確保し、大きい計算時間刻みによって計算負荷を減少させることが可能となる。また、粒子成長計算では、小さい計算時間刻みによって精度を確保し、微小領域内のみを計算することによって計算負荷を大幅に減少させることが可能となる。従って、本発明の実施の形態では、流れがある空間における粒子成長過程を短時間に精度良く模擬することが可能となると言える。
【００３９】
図３は、燃焼計算処理の流れを示すフローチャートである。燃焼計算処理では、空間（火炎内）における流体（気体）の流速分布及び粒子の濃度分布を計算する。
【００４０】
シミュレーション装置１の制御部１１は、記憶部１２に記憶されている条件データファイルから、条件データを入力する（ステップ１）。燃焼計算処理における条件データは、例えば、火炎の属性データ（計算モデル、バーナーガスの流量条件等）や、粒子の属性データ（モノマーの質量等）等である。ここで、モノマーとは、反応で生成される最小の粒子を意味する。一方、凝集粒子を構成するそれぞれの粒子は、一次粒子と呼ばれる。モノマーは相互に衝突し凝集粒子を形成し、凝集粒子間およびモノマーと凝集粒子の衝突が起き凝集粒子が成長して行く。焼結が起こると凝集粒子を構成するモノマー同士の融合が進み、凝集粒子を構成する粒子の径は大きくなる。したがって、モノマーや凝集粒子の衝突と焼結で形成された一次粒子の大きさは、様々である。そこで、粒子成長計算には、モノマーを最小単位の粒子として用いて後述する計算を実行する。
【００４１】
制御部１１は、メディア入出力部１３を介して、可搬記憶メディアに記憶されている条件データファイルから、条件データを入力しても良い。また、制御部１１は、通信制御部１４を介して、他のコンピュータの記憶部１２に記憶されている条件データファイルから、条件データを入力しても良い。また、制御部１１は、入力部１５を介してユーザによって入力されるデータを、条件データとして入力しても良い。
【００４２】
次に、制御部１１は、火炎の燃焼計算を行う（ステップＳ２）。制御部１１は、火炎内の流体（気体）の流速分布、火炎内の温度分布、火炎内の化学種の濃度分布を計算する。火炎の燃焼計算に用いられる方程式の一例は、次の通りである。
【００４３】
・連続の式（質量保存則）
【数１】
2014035614.tif 000003
ここで、Ｖ：速度ベクトル、ρ：密度、である。
【００４４】
・運動量の保存式（ナビィエ・ストークス方程式）
【数２】
2014035614.tif 000004
ここで、ｕｉ：ｘｉ方向の速度成分、ｉｉ：ｘｉ方向の単位ベクトルである。添え字のｉが取り得る値は、ｉ＝１、２、３（ｘ、ｙ、ｚ方向を示す。）である。
【００４５】
・エネルギーの方程式
【数３】
2014035614.tif 000005
ここで、Ｅ：エネルギー、ｋ：熱伝導度、Ｓ：反応によるエネルギー生成率、である。
【００４６】
・化学種ｉの保存式
【数４】
2014035614.tif 000006
ここで、ｍｉ：化学種ｉの質量分率、Ｄｉ：拡散係数、Ｒｉ：質量生成率、である。
【００４７】
制御部１１は、上述の連続の式、運動量の保存式、エネルギーの保存式、計算対象の化学種に対する化学種ｉの保存式を解くことで、火炎内の流体（気体）の流速分布、火炎内の温度分布、火炎内の化学種の濃度分布を計算する。
【００４８】
上述の各式は、本発明を理解し易いように、本発明の技術的思想を逸脱しない範囲で単純化したものである。ステップ２における火炎の燃焼計算は、市販の流体解析ソフトを用いて実行することが可能であるため、詳細な説明は省略する。
【００４９】
次に、制御部１１は、第１計算時間刻みによって定義される計算時刻をインクリメントする（ステップ３）。「第１計算時間刻み」は、燃焼計算用の計算時間刻みである。本発明の実施の形態では、後述する粒子成長計算処理用の計算時間刻みを「第２計算時間刻み」と表記し、両者を区別する。
【００５０】
尚、ステップ３は必ずしも必須ではない。例えば、定常状態における燃焼計算を行う場合、ステップ２の計算は１回であるから、計算時刻のインクリメントは必要ない。一方、燃焼開始時刻から過渡期を経て、定常状態までの一連の燃焼計算を行う場合、制御部１１は、燃焼開始時刻を「０」とし、ステップ２の計算ごとに計算時刻をインクリメントしていく。
【００５１】
次に、制御部１１は、終了条件を満たすか否かを確認する（ステップ４）。終了条件を満たさない場合（ステップ４のＮｏ）、制御部１１は、ステップ２から繰り返す。終了条件を満たす場合（ステップ４のＹｅｓ）、制御部１１は、ステップ５に進む。
【００５２】
定常状態における燃焼計算を行う場合、終了条件は「常に終了」である。また、燃焼開始時刻からａ秒後までの一連の燃焼計算を行う場合、終了条件は「計算時刻ｔ＞ａ」である。また、燃焼開始時刻から、指標ｂの値が閾値Ｂを超えるまでの一連の燃焼計算を行う場合、終了条件は「指標ｂの値＞Ｂ」である。このように、終了条件は、ユーザが任意に決定することができる。
【００５３】
次に、制御部１１は、空間内を流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により流体（気体）の流速分布に基づいて、微小領域の軌跡を計算する（ステップ５）。微小領域は、火炎内の環境を精度良く模擬するために必要な数及び大きさを定義する。そして、制御部１１は、微小領域ごとに、微小領域の軌跡を計算する。以下では、微小領域の軌跡を「流線」と表記することもある。
【００５４】
次に、制御部１１は、結果データを出力する（ステップ６）。制御部１１は、例えば、記憶部１２に結果データファイルを出力する。結果データファイルは、例えば、微小領域ごとに（流線ごとに）生成される。結果データファイルには、少なくとも、時刻ごとに、「微小領域の位置」、「微小領域内の温度」、「微小領域の移動速度」、「微小領域の断面積」、「微小領域内の粒子の濃度分布」が含まれる。また、「微少領域内の気体の密度」、「微少領域内の気体の分子量」もしくは「微少領域内の気体の平均分子量」、「微少領域内の気体の粘度」などが含まれていてもよい。「微小領域の位置」の時系列データが、微小領域の軌跡に相当する。反応によってモノマーが生成されるような場合、「微小領域内の粒子の濃度分布」は、モノマーの生成速度（単位時間、単に体積当たりの発生質量）の分布となる。その他、結果データファイルには、流体の物理量や粒子の物理量が含まれても良い。
【００５５】
後述する粒子成長計算処理では、ラグランジュ記述に基づく方法を用いるため、計算ステップごとに、微小領域内の温度、微小領域内の粒子の濃度分布等を与える必要がある。前述の通り、図３に示す燃焼計算処理によってこれらの値を算出するので、シミュレーション装置１は、後述する粒子成長計算処理を実行することが可能となる。
【００５６】
図４は、粒子成長計算処理の流れを示すフローチャートである。シミュレーション装置１の制御部１１は、記憶部１２に記憶されている燃焼計算の結果データファイルから、燃焼計算の結果データを入力する（ステップ１１）。
【００５７】
次に、制御部１１は、温度分布に基づいて流線上における温度を示す時間の関数を決定し、粒子の濃度分布（モノマーの生成速度の分布）に基づいて流線上における粒子の生成位置を示す時間の関数を決定する（ステップ１２）。
【００５８】
次に、制御部１１は、流線上における粒子の生成位置を示す時間の関数及び流線上における温度を示す時間の関数に基づいて、粒子成長計算を行う（ステップ１３）。粒子成長計算に必要な物理量は、流線上における粒子の生成位置（＝生成速度が実質的に０ではない領域）と、温度分布である。但し、衝突過程のモデルによっては、粒子のバックグラウンドとなる気体分子の平均自由行程が必要な場合もある。尚、気体の平均自由行程は、燃焼計算（流体計算）から直接求めることができないため、制御部１１は、燃焼計算（流体計算）とは別に、気体の粘度等から推定する。具体的には、別途、気体分子の平均自由行程の推定プログラムを用いて計算しても良いし、燃焼計算（流体計算）のポスト処理として実行しても良い。粒子成長計算に用いられる方程式の一例は、次の通りである。各方程式は、ｓ個のモノマーからなる粒子の個数の変化を表す式である。ここでは、説明を簡単にするために、焼結による粒子の融合を考慮しない式を掲げる。
【００５９】
・閉空間において衝突だけを考える場合（ラグランジュ記述）
【数５】
2014035614.tif 000007
ここで、ｔ：衝突が始まってからの時間、ｎｓ：ｓ個のモノマーからなる粒子の数密度、βｒ、ｓ：ｒ個のモノマーからなる粒子とｓ個のモノマーからなる粒子の衝突頻度関数、である。式（５）の右辺第１項は、足し合わせてモノマーの数がｓ個になる粒子ペアの衝突によりｓ個のモノマーからなる凝集粒子の数密度が増加することを表している。式（５）の右辺第２項は、ｓ個のモノマーからなる凝集粒子が他の粒子又は凝集粒子と衝突して数密度が減少することを表している。
【００６０】
・閉空間において衝突と拡散を考える場合（ラグランジュ記述）
【数６】
2014035614.tif 000008
ここで、Ｄｓ：粒子の拡散係数、である。添え字のｉが取り得る値は、ｉ＝１、２、３（ｘ、ｙ、ｚ方向を示す。）である。式（６）の右辺第３項は、拡散を意味している（拡散項）。
【００６１】
粒子同士の衝突だけを考える場合、制御部１１は、式（５）の方程式を解き、粒子成長計算を行う。また、粒子同士の衝突に加えて、粒子の拡散を考える場合、制御部１１は、式（６）の方程式を解き、粒子成長計算を行う。参考として、以下にオイラー記述の方程式を記載しておく。
【００６２】
・空間に固定された、ある微小領域での粒子数の変化を記述する式（オイラー記述）
【数７】
2014035614.tif 000009
ここで、ｕｉ：空間内を流れる流体の速度、である。式（７）の左辺第２項は、計算対象微小領域に隣接する微小領域から流入する粒子の数を表す項である。この項は、対流項又は移流項と呼ばれている。
【００６３】
尚、ここでは、式（６）、式（７）の記述の一部は、以下のアインシュタインの総和規約の記述形式を取るものとする。
【数８】
2014035614.tif 000010
具体的には、式（６）の右辺第３項、式（７）の左辺第２項及び右辺第３項は、式（８）のアインシュタインの総和規約に従って記述している。
【００６４】
式（７）に示すように、オイラー記述に基づく方法では、粒子の移動に関係する項が方程式に陽に存在するので、方程式が複雑になり、計算負荷が大きくなる。一方、式（５）及び（６）に示すように、ラグランジュ記述に基づく方法では、移動に関係する項が方程式に陽に存在しないので、方程式が簡単になり、計算負荷が小さくなる。焼結による粒子の融合を考慮した場合も、ラグランジュ記述に基づく方法の場合に、方程式が簡単になり、計算負荷が小さくなることは同様である。
【００６５】
図４の説明に戻る。次に、制御部１１は、第２計算時間刻みによって定義される計算時刻をインクリメントする（ステップ１４）。前述の通り、「第２計算時間刻み」は、粒子成長計算処理用の計算時間刻みである。
【００６６】
次に、制御部１１は、終了条件を満たすか否かを確認する（ステップ１５）。終了条件は、燃焼計算処理と同様、ユーザが任意に決定することができる。終了条件を満たさない場合（ステップ１５のＮｏ）、制御部１１は、ステップ１３から繰り返す。終了条件を満たす場合（ステップ１５のＹｅｓ）、制御部１１は、結果データファイルを出力し（ステップ１６）、処理を終了する。
【００６７】
粒子成長計算処理の結果データファイルには、少なくとも、「粒子の物理量」が含まれる。「粒子の物理量」の具体例としては、例えば、粒子数、粒子のバルク密度、モノマーの質量、表面積、体積、区分領域ごとの粒子数や粒子の総体積流量等である。
【００６８】
以上、図３及び図４に示すように、シミュレーション装置１は、少なくとも空間における流体の流速分布及び粒子の濃度分布、気体の温度分布を計算する流体計算ステップと、空間内を流れに沿って移動する微小領域を定義し、ラグランジュ記述に基づく方法により流体の流速分布に基づいて微小領域の軌跡を計算する軌跡計算ステップと、少なくとも、粒子の濃度分布に基づいて微小領域の軌跡上における粒子の生成位置を示す時間の関数及び温度分布に基づいて微小領域の軌跡上における温度を示す時間の関数を決定し、微小領域の軌跡上における粒子の生成位置を示す時間の関数に基づいて粒子の物理量を計算する粒子計算ステップと、を含むシミュレーション方法を実行する。
【００６９】
これによって、シミュレーション装置１は、流れがある空間における粒子成長過程を短時間に精度良く模擬することができる。具体的には、従来技術では、火炎の計算領域全体に対して粒子の成長を計算する必要があったが、本発明の実施の形態におけるシミュレーション方法では、粒子が生成する領域を通過する流線だけを計算対象にすれば、計算精度を確保できるので、計算時間を大幅に短縮することができる。また、本発明の実施の形態におけるシミュレーション方法では、燃焼計算と粒子成長計算が独立であるため、双方を計算精度が確保される最適な計算時間刻みで実行できるので、この点においても計算時間を短縮することができると言える。
【００７０】
但し、本発明の実施の形態におけるシミュレーション方法が成立するためには、以下の条件を満たす必要がある。
（条件１）粒子の拡散が無視できる。
（条件２）粒子が流体の運動にすぐに応答する。
（条件３）粒子の運動が流体の流れに影響を与えない。
【００７１】
まず、（条件１）について詳細に説明する。粒子の拡散が無視できない場合、微小領域内で粒子数が保存せず、計算が成り立たない。そこで、本発明の実施の形態におけるシミュレーション方法を適用する際には、粒子間の衝突による粒子数の変化と拡散による粒子数の変化との比較検証を行うことが望ましい。比較検証の結果、衝突による粒子数の変化に比べて、拡散による粒子数の変化が小さければ（例えば、１％未満等）、拡散は実質的に無視できると判断し、本発明の実施の形態におけるシミュレーション方法を適用しても問題がないと言える。尚、本発明者は、少なくともＶＡＤ法の火炎について比較検証を行い、拡散が実質的に無視できることを確認済である。
【００７２】
次に、（条件２）について詳細に説明する。気体中を運動する粒子は気体との速度差に比例した抵抗力を受ける。この抵抗力により、十分な時間の経過の後、粒子は気体と同じ速度で運動するようになる。静止している粒子が気体の速度の約６３％に到達するまでの時間を応答時間τvと呼ぶ。応答時間終端速度τvは
【数９】
2014035614.tif 000011
と書くことができる。ここでｍｐは粒子の質量、ｄｐは粒子の直径、ρｐは粒子の密度、μｃは気体の粘度である。上式に、ナノ粒子としては大きめの直径であるｄｐ＝１０−６ｍ、石英のバルク密度に近い値ρｐ＝２０００ｋｇ／ｍ３、火炎内の高温領域のおおよその粘度μｃ＝１０−４Ｐａ・ｓを代入すると、τv＝１０−６秒となる。これは火炎内の流速（数ｍ／ｓ）と比較すると非常に短い時間である。また、応答時間は直径の２乗に比例するので、粒子径が小さいほど短くなる。したがって、ナノｍサイズの粒子の応答時間はこれより短い。したがって、火炎内ではナノ粒子は流れにすぐに応答することがわかる。
【００７３】
以上より、（条件２）は、ナノ粒子であれば一般的に成り立つ条件であり、ＶＡＤ法や噴霧燃焼法などを含む多くの系において成り立つ。
【００７４】
次いで、（条件３）について詳細に説明する。また、粒子の運動が流体の運動に影響を与えるかを判断する指標として、交換パラメータ
【数１０】
2014035614.tif 000012
が用いられる。この交換パラメータが1よりもはるかに小さければ、粒子の運動が流体の運動に与える影響が小さいと見なすことができる。ここで、ρｄ、ρｃ（式（１０）ではρの記号に上線が付されているが、文章中では使用できない記号であるため、文章中では省略する。）は粒子と気体の嵩密度（単位体積の空間内にある粒子と気体の質量）、Ｓｔはストークス数でτｖ／τｆと書くことができる、τｆは流れの特性時間である。流れの特性時間の決め方には任意性があるが、考えている流れ場の代表的な長さを流れ場の平均的な流速または最高流速で割った値を使用する場合が一般的である。このとき流れ場の代表長さはバーナーの直径などを取ることが多い。ナノ粒子の場合、前段落で述べたようにτvは火炎の大きさや流速から考えられる流れの特性時間に比べて非常に小さいので、Ｓｔは０と考えることができる、したがって、ρｄ／ρｃ＜＜１が流体の運動に粒子の運動が与える影響の大小を判断する基準となる。
【００７５】
ρｄ／ρｃは粒子と気体の嵩密度の比であるので、粒子となる分子の量とそれ以外の分子の量から計算し評価することができる。そして、ρｄ／ρｃ＜＜１であれば（例えば０．２以下）、本発明の実施の形態におけるシミュレーション方法を適用しても問題がないと言える。尚、本発明者は、少なくともＶＡＤ法の火炎について比較検証を行い、流体の運動に粒子の運動が与える影響が無視できることを確認済である。
【００７６】
次に、図５及び図６を参照しながら、本発明の実施の形態における変形例を説明する。変形例では、シミュレーション装置１は、粒子計算ステップにおいて、計算対象の微小領域に含まれる粒子の密度を、計算対象の微小領域の軌跡と隣接する微小領域の軌跡との間隔、及び計算対象の微小領域の移動速度を用いて補正する。
【００７７】
まず、変形例の技術的意義について説明する。本発明の実施の形態におけるシミュレーション方法では、流れの影響を受けて、計算対象の流線と隣接する流線との間隔、及び計算対象の微小領域の移動速度が変化する。流線間の間隔が増加すると、流線に対して垂直な面の面積を大きくする作用がある。また、移動速度が増加すると、微小領域を流線の進行方向に引き延ばす作用がある。両方の作用とも、微小領域内の粒子の密度を下げる効果がある。このため、微小領域内の粒子密度を正しく補正しないと、正確な粒子数が得られない場合がある。そこで、シミュレーション装置１は、後述する粒子密度補正処理を実行し、精度良く正確な粒子数を算出する。
【００７８】
図５は、粒子密度補正処理を説明する図である。以下では、シミュレーション装置１が実行する粒子密度補正処理について、流れが１方向で２次元軸対称とみなせる層流火炎の場合を例に説明する。
【００７９】
図５（ａ）は、各計算ステップにおいて、計算対象流線に垂直な面を基準とする場合である。まず、制御部１１は、計算対象流線３１に隣接する２つの隣接流線３２、３３を特定する。隣接流線３２は、層流火炎の中心軸３０から見て外側寄りに計算対象流線３１と隣接している。隣接流線３３は、層流火炎の中心軸３０から見て中心寄りに計算対象流線３１と隣接している。
【００８０】
次に、制御部１１は、各計算ステップ（図５（ａ）では、ｔｉ、ｔｊ、ｔｋの３つを例示）において、計算対象流線３１に垂直な面（図５（ａ）では、４１、４２、４３の３つを例示）を決定する。
【００８１】
次に、制御部１１は、各計算ステップにおいて、この２つの隣接流線３２、３３と計算対象流線３１との中点をそれぞれ求めていき、それぞれの中点を結んだ２つの中点線５１、５２を決定する。
【００８２】
そして、制御部１１は、２つの中点線５１、５２で挟まれた斜線領域を計算対象の微小領域６１とする。この例の場合、火炎が２次元軸対称なので、計算対象の微小領域６１の断面は円環状となる。
【００８３】
制御部１１は、この円環の面積を各計算時間ステップで計算し、その値にその時間ステップでの微小領域の移動速度を乗じて、その時間ステップでの補正係数とする。この補正係数は、単位時間あたりに掃引する体積に相当するので、微小領域の粒子の密度に補正係数を掛けた値は微小領域内の粒子の「実数」である。粒子の「実数」は時間ステップ間で保存しなければならない。そこで、制御部１１は、この値を次の時間ステップに引き継ぐために、次の時間ステップの補正係数で割って、次の時間ステップの粒子の密度とする。反応によって新たな粒子が生成される場合、制御部１１は、この段階で、反応によって生成される新たな粒子を加える。
【００８４】
以上、シミュレーション装置１は、微小領域の軌跡に垂直な面上であって、計算対象の流線と隣接する流線との間隔から定まる所定の範囲の面積を算出し、所定の範囲の面積と計算対象の微小領域の移動速度の積を補正係数とし、計算対象の微小領域に含まれる粒子の実数を補正係数で割った値を、計算対象の微小領域に含まれる粒子の密度としてもよい。これによって、シミュレーション装置１は、精度良く正確な粒子数を算出することができる。
【００８５】
図５（ｂ）は、各計算ステップにおいて、流れの主流方向に垂直な面を基準とする場合である。流れの主流方向に垂直な面を基準とすることによって、計算を簡略化することができる。
【００８６】
前述の図５（ａ）の例では、制御部１１は、計算対象流線３１に垂直な面を決定したが、図５（ｂ）の例では、制御部１１は、流れの主流方向、すなわち層流火炎の中心軸３０の軸方向に垂直な面（図５（ｂ）では、４４、４５、４６の３つを例示）を決定する。それ以外は、前述の図５（ａ）の例と同様である。
【００８７】
以上、シミュレーション装置１は、流れの主流方向に垂直な面上であって、計算対象の流線と隣接する流線との間隔から定まる所定の範囲の面積を算出し、所定の範囲の面積と計算対象の微小領域の移動速度の積を補正係数とし、計算対象の微小領域に含まれる粒子の実数を補正係数で割った値を、計算対象の微小領域に含まれる粒子の密度としてもよい。これによって、シミュレーション装置１は、精度良く正確な粒子数を算出することができる。
【００８８】
図６は、シリカ粒子の質量流量の比較図である。計算対象は、四塩化珪素を加水分解してシリカ粒子を成長させる火炎である。（Ａ）は、流体計算から直接算出した質量流量である。（Ｂ）は、粒子密度補正処理を行わず、単純に密度の値を次の時間ステップに引き継いで算出した質量流量である。（Ｃ）は、粒子密度補正処理を行って算出した質量流量である。図６を見ると分かるように、粒子密度補正処理を行うことによって、流体計算から直接算出した値に近づけることができている。
【００８９】
以上、シミュレーション装置１が実行するシミュレーション処理について、火炎中での粒子成長計算を例に挙げて説明した。火炎中での粒子成長計算としては、ＶＡＤ法やリチウムイオン電池の正極活物質の生成などがあるが、これらの例に限定されるわけではない。また、粒子成長計算処理の前提となる処理は、燃焼計算処理に限定されない。本発明は、粒子を構成する分子を生成させるための十分なエネルギーが与えられる環境であれば、どのような系に対しても適用できる。例えば、火炎だけでなく、電気炉などが考えられる。また、例えば、室温下でも反応が進行する物質も存在する。また、本発明は、反応により粒子を構成する分子を生成する場合だけではなく、粒子が成長する環境に微小な粒子を流入させ、それらの成長を模擬することも可能である。
【００９０】
以上、添付図面を参照しながら、本発明に係るシミュレーション方法等の好適な実施の形態について説明したが、本発明は前述の例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【符号の説明】
【００９１】
１………シミュレーション装置
１１………制御部
１２………記憶部
１３………メディア入出力部
１４………通信制御部
１５………入力部
１６………表示部
１７………周辺機器Ｉ／Ｆ部
１８………バス
２１………プリプロセッサ
２２………ソルバー
２３………ポストプロセッサ
２４………条件データ
３０………中心軸
３１………計算対象流線
３２、３３………隣接流線
４１、４２、４３………計算対象流線に垂直な面
４４、４５、４６………計算対象流線に垂直な面
５１、５２………中点線
６１………計算対象の微小領域
【図１】
2014035614.tif 000013
【図２】
2014035614.tif 000014
【図３】
2014035614.tif 000015
【図４】
2014035614.tif 000016
【図５】
2014035614.tif 000017
【図６】
2014035614.tif 000018