(57)【要約】【課題】現実的な時間でガラス溶融炉内の全ての白金族粒子の挙動を模擬することができるシミュレーション方法及びシミュレーション装置を提供する。【解決手段】液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより上記複数の白金族粒子の挙動を解析するシミュレーション方法であって、上記白金族粒子のいくつかで構成される群を粒子径が上記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を上記解析対象データに施す粗視化モデル生成工程を有し、該粗視化モデル生成工程によって得られた粗視化モデルデータに上記離散要素法に基づく演算処理を施すという構成を採用する。【選択図】図１２
2010107314.tif 000002
【特許請求の範囲】
【請求項１】
液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより前記複数の白金族粒子の挙動を解析するシミュレーション方法であって、
前記白金族粒子のいくつかで構成される群を粒子径が前記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を前記解析対象データに施す粗視化モデル生成工程を有し、
該粗視化モデル生成工程によって得られた粗視化モデルデータに前記離散要素法に基づく演算処理を施すことを特徴とするシミュレーション方法。
【請求項２】
前記粗視化モデル生成工程では、前記固液混相流において粒子径の大きさに依存して作用する所定の力を、粒子径の大きな前記粗視化粒子と前記粗視化粒子を構成する粒子径の小さな複数の前記白金族粒子との間で一致させるモデル化を施すことを特徴とする請求項１に記載のシミュレーション方法。
【請求項３】
前記所定の力は、ファンデルワールス力を含むことを特徴とする請求項２に記載のシミュレーション方法。
【請求項４】
前記粗視化モデル生成工程では、前記白金族粒子及び前記粗視化粒子を球形でモデル化し、前記粗視化粒子の粒子径が前記白金族粒子の粒子径よりもＬ倍大きく、前記粗視化粒子をＬ３個の前記白金族粒子で構成する場合、前記粗視化粒子に作用するファンデルワールス力をＦＶＷ ｉとし、前記粗視化粒子を構成する前記白金族粒子のそれぞれに作用するファンデルワールス力をＦＶＷ Ｏとしたときに、
【数１】
2010107314.tif 000047
の関係を満足させることを特徴とする請求項３に記載のシミュレーション方法。
【請求項５】
前記所定の力は、電気二重層力を含むことを特徴とする請求項２〜４のいずれか一項に記載のシミュレーション方法。
【請求項６】
前記粗視化モデル生成工程では、前記白金族粒子及び前記粗視化粒子を球形でモデル化し、前記粗視化粒子の粒子径が前記白金族粒子の粒子径よりもＬ倍大きく、前記粗視化粒子をＬ３個の前記白金族粒子で構成する場合、前記粗視化粒子に作用する電気二重層力をＦｅ ｉとし、前記粗視化粒子を構成する前記白金族粒子のそれぞれに作用する電気二重層力をＦｅ Ｏとしたときに、
【数２】
2010107314.tif 000048
の関係を満足させることを特徴とする請求項５に記載のシミュレーション方法。
【請求項７】
液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより前記複数の白金族粒子の挙動を解析するシミュレーション装置であって、
前記白金族粒子のいくつかで構成される群を粒子径が前記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を前記解析対象データに施す粗視化モデル生成手段を有し、
該粗視化モデル生成手段によって得られた粗視化モデルデータに前記離散要素法に基づく演算処理を施すことを特徴とするシミュレーション装置。
【請求項８】
前記粗視化モデル生成手段では、前記固液混相流において粒子径の大きさに依存して作用する所定の力を、粒子径の大きな前記粗視化粒子と前記粗視化粒子を構成する粒子径の小さな複数の前記白金族粒子との間で一致させるモデル化を施すことを特徴とする請求項７に記載のシミュレーション装置。
【請求項９】
前記所定の力は、ファンデルワールス力を含むことを特徴とする請求項８に記載のシミュレーション装置。
【請求項１０】
前記粗視化モデル生成手段では、前記白金族粒子及び前記粗視化粒子を球形でモデル化し、前記粗視化粒子の粒子径が前記白金族粒子の粒子径よりもＬ倍大きく、前記粗視化粒子をＬ３個の前記白金族粒子で構成する場合、前記粗視化粒子に作用するファンデルワールス力をＦＶＷ ｉとし、前記粗視化粒子を構成する前記白金族粒子のそれぞれに作用するファンデルワールス力をＦＶＷ Ｏとしたときに、
【数３】
2010107314.tif 000049
の関係を満足させることを特徴とする請求項９に記載のシミュレーション装置。
【請求項１１】
前記所定の力は、電気二重層力を含むことを特徴とする請求項８〜１０のいずれか一項に記載のシミュレーション装置。
【請求項１２】
前記粗視化モデル生成手段では、前記白金族粒子及び前記粗視化粒子を球形でモデル化し、前記粗視化粒子の粒子径が前記白金族粒子の粒子径よりもＬ倍大きく、前記粗視化粒子をＬ３個の前記白金族粒子で構成する場合、前記粗視化粒子に作用する電気二重層力をＦｅ ｉとし、前記粗視化粒子を構成する前記白金族粒子のそれぞれに作用する電気二重層力をＦｅ Ｏとしたときに、
【数４】
2010107314.tif 000050
の関係を満足させることを特徴とする請求項１１に記載のシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びシミュレーション装置に関するものである。
【背景技術】
【０００２】
高レベル放射性廃液の処理方法として、下記特許文献１に開示されているような、ガラス溶融炉において高レベル放射性廃液をガラス原料と一緒に溶融して溶融ガラスを生成し、
その後、固化する技術が知られている。
【０００３】
ガラス溶融炉の運転条件によっては、ガラス溶融炉の内部において高レベル放射性廃液に含まれる金属成分、特に白金族粒子（ルテニウム（Ｒｕ）、ロジウム（Ｒｈ）、パラジウム（Ｐｄ）等）が析出し、そのガラス溶融炉の内部に堆積する可能性がある。ガラス溶融炉は、そのガラス溶融炉の内部に設けられた電極間に電流を流し、発生した熱（ジュール熱）によって、ガラス原料を溶融する。ガラス溶融炉の内部に白金族粒子が堆積すると、その白金族粒子が短絡等を引き起こし、円滑な溶融処理を妨げる可能性がある。そのような不具合の発生を抑制するため、溶融したガラスとの固液混相流における白金族粒子の挙動を把握し、ガラス溶融炉の運転条件の最適化を図る必要がある。
【０００４】
固液混相流の内部機構の解明にあたり、ガラス溶融炉のような高レベル放射性廃液を扱う大規模体系では、濃度や構成粒子の移動距離の計測が容易ではなく、実験によるアプローチが困難になることがあるため、数値シミュレーションによるアプローチが期待されている（例えば非特許文献１および非特許文献２参照）。
【０００５】
これまでに開発されてきた固液混相流の解析手法として、液相と固相とをオイラー的手法で解析するオイラー・オイラー的手法、液相をオイラー的手法で固相をラグランジュ的手法で解析するオイラー・ラグランジュ的手法および液相と固相とをラグランジュ的手法で解析するラグランジュ・ラグランジュ的手法がある。これまでは、主として液相バルク領域における固相の挙動が研究されてきた。
固相のモデル化においては、ラグランジュ的手法のひとつである離散要素法（Discrete Element Method; ＤＥＭ、以下ＤＥＭと記すことがある）が固体粒子同士の接触力を精度よく評価する場合に用いられる（例えば非特許文献３参照）。
【特許文献１】特開２００２−７１８９１号公報
【非特許文献１】中島正義、三浦昭彦、堺公明、他、「ガラス固化溶融炉の解析シミュレーションシステムの開発（１）シミュレーションシステムにおける解析手法」、日本原子力学会、日本原子力学会年会・大会予稿集、ｐｐ６４６（２００６）
【非特許文献２】三浦昭彦、中島正義、堺公明、他、「ガラス固化溶融炉の解析シミュレーションシステムの開発（２）ＴＶＦガラス溶融炉運転時の物理挙動」、日本原子力学会、日本原子力学会年会・大会予稿集、ｐｐ６４６（２００６）
【非特許文献３】酒井幹夫、越塚誠一、「固気混相流における離散要素法の粗視化手法の開発」、粉体工学会、粉体工学会会誌・第４５巻、ｐｐ１２（２００８）
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、離散要素法は、まだ実用的に確立されたものではなく、特にガラス溶融炉のような大規模体系への適用などの課題が残る。大規模体系への適用の最も大きな課題のひとつが計算粒子数である。最新のＰＣ（Personal Computer）を用いて現実的な時間で計算する場合、計算粒子数はせいぜい数十万個程度であり、この数十万個の粒子数は、砂糖スプーン一杯分に相当すると言われている。このため、ガラス溶融炉内の白金族粒子のような十億個超の粒子数が必要と言われる大規模体系（一例としてガラス溶融炉内の容積は１立方メートル程度）のシミュレーションは、現実的な時間で計算することができなかった。したがって、従来の手法では、ガラス溶融炉内の全ての白金族粒子の挙動が模擬されないことが課題となっていた。
【０００７】
本発明は、上記課題点に鑑みてなされたものであり、現実的な時間でガラス溶融炉内の全ての白金族粒子の挙動を模擬することができるシミュレーション方法及びシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の課題を解決するために、本発明は、液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより上記複数の白金族粒子の挙動を解析するシミュレーション方法であって、上記白金族粒子のいくつかで構成される群を粒子径が上記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を上記解析対象データに施す粗視化モデル生成工程を有し、該粗視化モデル生成工程によって得られた粗視化モデルデータに上記離散要素法に基づく演算処理を施すという構成を採用する。
このような構成を採用することによって、本発明では、白金族粒子の群をモデル化し、計算上大きな粗視化粒子を用いて白金族粒子体系を模擬するため、実際よりも少ない粒子数でガラス溶融炉内の白金族粒子の挙動を解析することが可能となる。
【０００９】
また、本発明は、上記粗視化モデル生成工程では、上記固液混相流において粒子径の大きさに依存して作用する所定の力を、粒子径の大きな上記粗視化粒子と上記粗視化粒子を構成する粒子径の小さな複数の上記白金族粒子との間で一致させるモデル化を施すという構成を採用する。
このような構成を採用することによって、本発明では、粒子径の大きさに依存して作用する流体力や接触力等の所定の力を白金族粒子の群と粗視化粒子との間で一致させるため、実際の白金族粒子の挙動を粗視化粒子で精度良く模擬することが可能となる。
【００１０】
また、本発明は、上記所定の力は、ファンデルワールス力を含むという構成を採用する。
このような構成を採用することによって、本発明では、ファンデルワールス力を考慮し、白金族粒子同士の付着力をモデル化することで、溶融ガラス中の白金族粒子の凝集・堆積の挙動を解析することが可能となる。
【００１１】
また、本発明は、上記粗視化モデル生成工程では、上記白金族粒子及び上記粗視化粒子を球形でモデル化し、上記粗視化粒子の粒子径が上記白金族粒子の粒子径よりもＬ倍大きく、上記粗視化粒子をＬ３個の上記白金族粒子で構成する場合、上記粗視化粒子に作用するファンデルワールス力をＦＶＷ ｉとし、上記粗視化粒子を構成する上記白金族粒子のそれぞれに作用するファンデルワールス力をＦＶＷ Ｏとしたときに、下式（ａ）の関係を満足させるという構成を採用する。
【数１】
2010107314.tif 000003
このような構成を採用することによって、本発明では、粗視化粒子に作用するファンデルワールス力と粗視化粒子を構成する複数の白金族粒子に作用するファンデルワールス力とのバランスをとることができる。
【００１２】
また、本発明は、上記所定の力は、電気二重層力を含むという構成を採用する。
このような構成を採用することによって、本発明では、電気二重層力を考慮し、白金族粒子同士の斥力をモデル化することで、溶融ガラス中の白金族粒子の凝集・堆積の挙動を解析することが可能となる。
【００１３】
また、本発明は、上記粗視化モデル生成工程では、上記白金族粒子及び上記粗視化粒子を球形でモデル化し、上記粗視化粒子の粒子径が上記白金族粒子の粒子径よりもＬ倍大きく、上記粗視化粒子をＬ３個の上記白金族粒子で構成する場合、上記粗視化粒子に作用する電気二重層力をＦｅ ｉとし、上記粗視化粒子を構成する上記白金族粒子のそれぞれに作用する電気二重層力をＦｅ Ｏとしたときに、の関係を満足させるという構成を採用する。
【数２】
2010107314.tif 000004
このような構成を採用することによって、本発明では、粗視化粒子に作用する電気二重層力と粗視化粒子を構成する複数の白金族粒子に作用する電気二重層力とのバランスをとることができる。
【００１４】
また、本発明は、液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより上記複数の白金族粒子の挙動を解析するシミュレーション装置であって、上記白金族粒子のいくつかで構成される群を粒子径が上記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を上記解析対象データに施す粗視化モデル生成手段を有し、該粗視化モデル生成手段によって得られた粗視化モデルデータに上記離散要素法に基づく演算処理を施すという構成を採用する。
このような構成を採用することによって、本発明では、白金族粒子の群をモデル化し、計算上大きな粗視化粒子を用いて白金族粒子体系を模擬するため、実際よりも少ない粒子数でガラス溶融炉内の白金族粒子の挙動を解析することが可能となる。
【発明の効果】
【００１５】
本発明によれば、液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより上記複数の白金族粒子の挙動を解析するシミュレーション方法であって、上記白金族粒子のいくつかで構成される群を粒子径が上記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を上記解析対象データに施す粗視化モデル生成工程を有し、該粗視化モデル生成工程によって得られた粗視化モデルデータに上記離散要素法に基づく演算処理を施すという構成を採用することによって、白金族粒子の群をモデル化し、計算上大きな粗視化粒子を用いて白金族粒子体系を模擬するため、実際よりも少ない粒子数でガラス溶融炉内の白金族粒子の挙動を解析することが可能となる。
したがって、本発明は、現実的な時間でガラス溶融炉内の全ての白金族粒子の挙動を模擬することができるシミュレーション方法を提供できる効果がある。
【発明を実施するための最良の形態】
【００１６】
以下、本発明の一実施形態を図面および数式に基づいて説明する。先ず、本発明におけるシミュレーション装置の構成およびシミュレーション方法のフローの概略について簡単に説明し、次いで本発明のシミュレーション方法に係る白金族粒子のモデル化や計算手法について詳細に説明する。
図１は、本発明の実施形態におけるシミュレーション装置１の概略構成を示すブロック図である。
図２は、本発明の実施形態におけるシミュレーション方法を示すフローチャートである。
シミュレーション装置１は、図１に示すように、ガラス溶融炉内部の液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の挙動を解析するコンピュータシステムであり、入力装置２、出力装置３、記憶装置４およびＣＰＵ５を有する。
【００１７】
入力装置２は、外部からデータをＣＰＵ５に入力するものであり、マウス、キーボードおよび通信装置等を有する。出力装置３は、ＣＰＵ５から出力されるシミュレーション結果を画像として出力するものであり、モニタ等の表示装置を有する。記憶装置４は、固液混相流のシミュレーションに必要なデータおよびプログラムを記憶するものであり、ハードディスク、ＲＯＭおよびＲＡＭ等を有する。この記憶装置４は、固液混相流の状態を示す解析対象データと、解析対象データに粗視化処理を施す粗視化プログラム（粗視化モデル生成手段）と、粗視化処理によって得られる粗視化モデルデータと、粗視化モデルデータに離散要素法に基づく演算処理を施す解析プログラムとを記憶する構成となっている。ＣＰＵ５は、入力装置２、出力装置３および記憶装置４との間でデータ授受を行う各種入出力インターフェイス回路等を有し、固液混相流の挙動を解析する演算処理を行う構成となっている。
【００１８】
図２に示すように、本実施形態のシミュレーション方法では、先ず、白金族粒子のいくつかで構成される群を粒子径が白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を解析対象データに施す（ステップＳＴ１：粗視化モデル生成工程）。シミュレーション装置１においては、入力装置２から固液混相流の解析指令を受けたＣＰＵ５が、記憶装置４の解析対象データに含まれる白金族粒子の固相データを引き出して粗視化プログラムを実行することにより、粗視化モデルデータを得る。
そして、ＣＰＵ５が、粗視化モデル生成工程で得られた粗視化モデルデータにＤＥＭに基づく演算処理を、さらに、液相データにＭＰＳ法（後述）に基づく演算処理を施す解析プログラムを実行し、シミュレーション結果を出力装置３に出力させることにより、固相である白金族粒子と液相である溶融ガラスとからなる固液混相流の挙動を解析する（ステップＳＴ２：固液混相流解析工程）。
【００１９】
以下、上記シミュレーション方法の詳細について説明する。
本実施形態のシミュレーション方法では、液相をラグランジュ的手法で固相をラグランジュ的手法で解析するラグランジュ・ラグランジュ的手法を採用する。具体的に本実施形態では、液相には、複数の流体粒子（以下、液相粒子と記すことがある）によりその挙動を模擬するMoving Particle Semi-implicit法（以下、ＭＰＳ法と記すことがある）を適用し、固相には、複数の固体粒子（以下、固相粒子と記すことがある）によりその挙動を模擬するＤＥＭを適用する。なお、ＭＰＳ法は、ラグランジュ型の計算手法であり、説明は後述するが、より詳細には公知技術文献（「粒子法」、越塚誠一、丸善（２００５））等に開示されている。また、以下の説明では、上記のように液相にＭＰＳ法を適用した一例について説明するが、液相には、ＳＰＨ法、有限差分法、有限体積法等を適用してもよい。
【００２０】
固相のＤＥＭシミュレーションにおいてガラス溶融炉の大規模体系の解析を実行するためには、計算負荷が大きいことや計算機メモリによる粒子数の制限があるため、本実施形態では、白金族粒子（以下、オリジナル粒子と記すことがある）の群を、オリジナル粒子よりも大きなモデル粒子（以下、粗視化粒子と記すことがある）で代表させる粗視化処理を施す。すなわち、現実的な時間でガラス溶融炉内の全ての白金族粒子の挙動を模擬するために、固体粒子のスケールをオリジナル粒子よりも大きく取り扱うモデル化が必要となる。この際、固体粒子を単に大きくしてはならない。すなわち、固相−液相間相互作用力である流体力や固相−固相間相互作用力である接触力、付着力および斥力等の力（所定の力）は、固体粒子の粒子径に依存するため、サブミクロン程度の粒子径の小さな白金族粒子の挙動を、粒子径の大きな粗視化粒子で精度よく模擬するためには、上記力を両者間で一致させる必要がある。
【００２１】
図３〜図５は、粒子径の大きな粗視化粒子と該粒子を構成する粒子径の小さな複数のオリジナル粒子との挙動の違いを示す例である。図３に示すように、矢印の方向にガスを噴きつけた場合、図３（ａ）に示す粗視化粒子に比べて、図３（ｂ）に示すオリジナル粒子は容易に移動する。また、図４に示すように、同体積の両粒子が壁に衝突した場合、図４（ａ）に示す粗視化粒子と図４（ｂ）に示すオリジナル粒子とでは、固体粒子の衝突頻度が異なり、両者のエネルギ消散が異なる。また、図５に示すように、図５（ａ）に示す粗視化粒子と図５（ｂ）に示すオリジナル粒子とでは、固体粒子の凝集のし易さ（付着力）が異なる。さらに、斥力の大きさを粗視化粒子とオリジナル粒子とで一致させて、両者間のバランスをとることが必要となる。
【００２２】
また、ラグランジュ・ラグランジュ的手法の大規模体系への応用や計算効率の観点から、液相および固相の両者の離散スケールを、計算上、同じスケールで取り扱うことが望まれる。ガラス溶融炉の大規模体系では、液相と固相との空間の離散化スケールが異なり、液相に比べて固相のスケールが著しく小さくなってしまうためである。
したがって、本実施形態では、液相と固相との離散スケールを、計算上、同じスケールとして取り扱い、粒子径の大きさに依存する上記力をオリジナル粒子と粗視化粒子との間でバランスをとるように粗視化処理を施す。以下、具体的なオリジナル粒子のモデル化や計算手法について説明する。
【００２３】
ラグランジュ・ラグランジュ的手法の固液混相流において、該固液混相流を構成する各粒子は液相もしくは固相のどちらかに属するものであり、オイラー・オイラー的手法のように各格子に両者のベクトルやスカラーをもつものとは異なる。そのため、ラグランジュ・ラグランジュ的手法では、液相粒子と固相粒子に、それぞれの運動方程式を与えて、計算する必要がある。
液相と固相の連続の式および運動方程式は、それぞれ、下式（１）〜下式（５）のように表される。
【００２４】
【数３】
2010107314.tif 000005
【００２５】
【数４】
2010107314.tif 000006
【００２６】
【数５】
2010107314.tif 000007
【００２７】
【数６】
2010107314.tif 000008
【００２８】
【数７】
2010107314.tif 000009
【００２９】
ここで、太字ｕ、ｐ、ρ、μ、太字ｇ、太字ＦＣ、Ｋ、εおよびｍは、それぞれ、流速ベクトル、圧力、密度、粘性率、重力加速度、接触力、流体抵抗係数、体積分率および質量である。添字について、ｌ、ｌｌ、ｓ、ｓｓおよびｓｌは、それぞれ、液相、液相-液相間、固相、固相-固相間および固相-液相間を意味する。
液相および固相のモデルの詳細について、以下に示す。
【００３０】
［液相］
液相については、圧力項、粘性項、流体力項および重力項を計算する。液相の運動方程式の式（２）には、微分演算子として勾配とラプラシアンが含まれる。さらに、圧力のポアッソン方程式においてもラプラシアンが現れる。これらにＭＰＳ法の勾配モデルおよびラプラシアンモデルを適用する。ここでは、圧力項および粘性項で使用する勾配モデルおよびラプラシアンモデルの概要を述べる。なお、ＭＰＳ法は、非圧縮性流れのラグランジュ型の計算手法である。ＭＰＳ法では、移流は粒子の移動を追跡することにより計算されるので、数値拡散の影響を排除できる。このため、水面形状の複雑な変化を高精度に計算できる。なお、流体力については、固相のモデリングに関する情報が必要になるため、固相の説明で述べる。
【００３１】
（重み関数）
ＭＰＳ法では、勾配やラプラシアンといった微分演算子に対して、粒子間相互作用モデルを用いて微分方程式を離散化する。粒子間相互作用モデルには、重み関数ｗを下式（６）のように導入する。
【００３２】
【数８】
2010107314.tif 000010
【００３３】
ここで、ｘおよびｒeは、粒子間距離および影響半径である。式（６）は、ｘがｒeよりも短い場合のみ相互作用することを意味する。
【００３４】
（粒子数密度）
粒子iにおける粒子数密度は、下式（７）から求められる。
【００３５】
【数９】
2010107314.tif 000011
【００３６】
ここで、太字ｘは粒子の位置ベクトルである。非圧縮性流れの場合、粒子数密度は初期状態に定めた粒子数密度ｎ０（一定値）にならなければならない。
【００３７】
（圧力項）
粒子ｉの位置におけるｐは、勾配モデルを用いて、下式（８）のように表される。
【００３８】
【数１０】
2010107314.tif 000012
【００３９】
勾配モデルは、粒子ｉとその近傍粒子ｊとの間で単純な勾配ベクトルを定義し、ｗを掛けて平均をとったものである。ここで、Ｄは空間の次元数である。
【００４０】
（粘性項）
粘性項は、ラプラシアンモデルを使用し、下式（９）のように計算される。
【００４１】
【数１１】
2010107314.tif 000013
【００４２】
ラプラシアンモデルにより、粒子ｉの物理量がその近傍の粒子ｊにｗの分布で分配される。ここで、係数λは、下式（１０）のように表される。
【００４３】
【数１２】
2010107314.tif 000014
【００４４】
［固相］
固相については、圧力項、流体力項、接触力項、付着力項、斥力項および重力項を計算する。圧力項の計算手順については液相のものと同じである。ここでは、接触力項、流体力項、付着力項および斥力項について述べる。
固相のＤＥＭシミュレーションでは、オリジナル粒子群を大きな粗視化粒子で代表させて解析を行う。ＤＥＭの粗視化モデルは、固体粒子に作用する力について、粗視化粒子がオリジナル粒子群の挙動を模擬できるようにモデル化する。ＤＥＭでは、図６に示すように固体粒子（図６において符号１０で示す）間に作用する接触力をバネ１１（バネ定数ｋ）、ダッシュポット１２（粘性減衰係数η）およびフリクションスライダ１３（摩擦係数μ）から構成されるモデルを用いて並進運動および回転運動を評価する。なお、ＤＥＭの粗視化モデルの詳細については、上記した非特許文献３を参照されたい。
【００４５】
（接触力）
本実施形態では、球形の固体粒子を扱う。オリジナル粒子よりもＬ倍大きな粗視化粒子を用いる場合、粗視化粒子にはＬ３個のオリジナル粒子が含まれる。例えば、図７に示すように、粗視化粒子の粒子径がオリジナル粒子の粒子径よりも２倍大きい場合、粗視化粒子には８個のオリジナル粒子が含まれる。
粗視化モデルでは、接触力について、粗視化粒子とオリジナル粒子群（平均）の運動エネルギが一致するとして、下式（１１）のように定式化する。
【００４６】
【数１３】
2010107314.tif 000015
【００４７】
ここで、ｍｓ、Ｉ、Ｌおよび太字ωは、それぞれ、固体粒子の質量、慣性モーメント、粗視化率および角速度ベクトルである。添字のＣＧＭおよびＯは、粗視化粒子およびオリジナル粒子を意味する。並進運動について、粗視化粒子とオリジナル粒子（平均）とが一致すると仮定する（図７参照）。回転運動のエネルギについて、粗視化粒子とオリジナル粒子群を一致させるにあたり、オリジナル粒子が各々の重心で回転することを仮定する（図８参照）。
【００４８】
式（１１）より、粗視化粒子にはオリジナル粒子のＬ３倍の接触力が働くことになる。また、粒子-流体間相互作用力、付着力、斥力項や重力について、粗視化粒子に作用するこれらの力が粗視化粒子領域に存在するオリジナル粒子群のものとバランスするようにモデル化する。その結果、粗視化粒子の運動方程式は、下式（１２）のように表される。
【００４９】
【数１４】
2010107314.tif 000016
【００５０】
ここで、太字ＦＤ、Ｖｓ、太字ＦＶＷ、太字Ｆｅおよび太字Ｆｇは、それぞれ、流体力、固体粒子の体積、ファンデルワールス力、斥力（電気二重層力）および重力である。式（１２）は、粗視化粒子に含まれるオリジナル粒子の数とその相互作用力の平均の積（粗視化粒子に含まれるオリジナル粒子群の相互作用力の総和に相当する）で表され、粗視化粒子の固体粒子iおよびｊが二体衝突すると、それに含まれる全てのオリジナル粒子（すなわち、粗視化粒子１個あたり、Ｌ３個のオリジナル粒子）が二体衝突したことに相当する（図９参照）。
回転運動について、オリジナル粒子と粗視化粒子との間には、下式（１３）の関係が成り立つ。
【００５１】
【数１５】
2010107314.tif 000017
【００５２】
ここで、太字Ｔおよびｒは、トルクおよび固体粒子の半径である。
オリジナル粒子体系と粗視化粒子体系との関係について、粗視化粒子とオリジナル粒子（平均）が一致すると仮定した式（１２）と式（１３）とにより、式（１１）は満たされる。また、粗視化モデルに作用する接触力は、上記のオリジナル粒子体系と粗視化粒子体系の関係から、法線方向成分、接触点で滑らない場合の接線方向成分および接触点で滑る場合の接線方向成分について、それぞれ、下式（１４）〜下式（１６）のように表される。
【００５３】
【数１６】
2010107314.tif 000018
【００５４】
【数１７】
2010107314.tif 000019
【００５５】
【数１８】
2010107314.tif 000020
【００５６】
ここで、ｋ、δ、ηおよび 太字ｄは、それぞれ、バネ定数、オーバーラップ幅、粘性減衰係数および並進運動による変位である。添字のｎおよびｔは、垂直方向成分および接線方向成分を意味する。なお、粗視化モデルでは、オリジナル固体粒子の物性値を使用する。
【００５７】
（流体力）
粒子-流体間の流体力について、粗視化した固体粒子ｉに作用する流体力は、下式（１７）から求められる。
【００５８】
【数１９】
2010107314.tif 000021
【００５９】
ただし、流体粒子の速度の太字ｕｌは、下式（１８）に示すように、重み関数に基づいて平均されたものを使用する。
【００６０】
【数２０】
2010107314.tif 000022
【００６１】
流体抵抗係数は、下式（１９）のように表され、粒子径にはオリジナルのものを設定する。ただし、粉体濃度が高い場合には他の係数（Ergunの式）を用いることもある。
【００６２】
【数２１】
2010107314.tif 000023
【００６３】
ここで、ＣＤおよびｄｓは、それぞれ、抗力係数および固体粒子の直径である。ＣＤは、粒子レイノルズ数Ｒｅｐの関数であり、下式（２０）のように表される。
【００６４】
【数２２】
2010107314.tif 000024
【００６５】
Ｒｅｐは、下式（２１）のように与えられる。
【００６６】
【数２３】
2010107314.tif 000025
【００６７】
式（２１）で示すＲｅｐ内の粒子径についてもオリジナルのものを使用する。
流体粒子に作用する流体力についても同様の手順により求められ、粒子および流体に作用する流体力を考慮することにより、two-way couplingが可能になる。
【００６８】
（付着力）
本実施形態では、付着力としてファンデルワールス力を考慮する。粗視化した固体粒子ｉに作用するファンデルワールス力は、下式（２２）のように与えられる。
【００６９】
【数２４】
2010107314.tif 000026
【００７０】
ここで、Ｈ、Ｘ、およびｈは、それぞれ、ハマカー定数、換算粒子径および粒子の表面間距離である。ＸＯは、下式（２３）のように与えられる。
【００７１】
【数２５】
2010107314.tif 000027
【００７２】
ＸＯｉおよびＸＯｊは、粗視化した固体粒子ｉおよび固体粒子ｊに含まれるオリジナル粒子径に基づく粒子径である。
【００７３】
（斥力）
本実施形態では、斥力（電気二重層力）を考慮する。粗視化した固体粒子ｉに作用する斥力は、下式（２４）のように与えられる。
【００７４】
【数２６】
2010107314.tif 000028
【００７５】
ただし、γおよびκは、それぞれ、下式（２５）および下式（２６）のように与えられる。
【００７６】
【数２７】
2010107314.tif 000029
【００７７】
【数２８】
2010107314.tif 000030
【００７８】
ここで、ｒｓ、ｎ、ｋｂ、Ｔ、ｈ、ｚ、ｅ、ε０、εｒおよびψ０は、それぞれ、粒子半径、電解質イオン数密度、ボルツマン定数、液相の溶媒温度、粒子表面間距離、電解質イオン価数、単位電荷、真空の誘電率、溶媒の比誘電率および電位である。
【００７９】
以上のように、流体力と固体粒子間の接触力、付着力および斥力（電気二重層力）にモデル化を施すことにより、液相と固相の空間の離散化スケールが異なる体系に対して、計算上、両者を同じスケールで扱うことができる。
ただし、高塩濃度または低温の場合では電気二重層が薄くなるため、電気二重層力が小さくなることから、該斥力を無視できることがある。
なお、ここで、上記力に加えてランダム力を考慮する構成であってもよい。
【００８０】
続いて、図１０を参照して、本実施形態における計算アルゴリズムについて説明する。なお、該計算アルゴリズムの概要は上記公知技術文献等に開示されているため、詳細は割愛することとする。
【００８１】
［液相］
ＭＰＳ法では、液相の運動方程式の時間積分過程において、ＳＭＡＣ法（A.A.Amsden and F.H.Harlow: A Simplified MAC Technique for Incompressible Fluid Flow Calculations. J.Comp.Phys. 6, 322 (1970)の公知技術文献参照）と同様の半陰的アルゴリズムである２段解法が用いられる。ここで、時間ステップｎにおける流れの計算について説明する（ステップＳＴ２１〜ステップＳＴ２９）。第１段階は、陽的なプロセスであり、重力項（ステップＳＴ２２）、粘性項（ステップＳＴ２３）および流体力（ステップＳＴ２４）に基づいて、粒子の仮の速度の太字ｕｌ＊および仮の位置のｘｌ＊が得られる（ステップＳＴ２５）。すなわち、ｕｌ＊およびｘｌ＊は、下式（２７）および下式（２８）のように表される。
【００８２】
【数２９】
2010107314.tif 000031
【００８３】
【数３０】
2010107314.tif 000032
【００８４】
陽的なプロセスが終了した段階の粒子数密度ｎ＊は、初期の粒子数密度ｎ０になっていないので、質量保存則が満たされない。そこで、第２段階では、圧力項と粒子数密度で記述された質量保存則から導かれた圧力のポアッソン方程式を陰的に計算して、圧力を更新する（ステップＳＴ２６〜ステップＳＴ２８）。ポアッソン方程式は、下式（２９）のように表される。
【００８５】
【数３１】
2010107314.tif 000033
【００８６】
速度の修正量は、下式（３０）のように表される。
【００８７】
【数３２】
2010107314.tif 000034
【００８８】
修正後の速度と位置は、それぞれ、下式（３１）と下式（３２）のようになる。
【００８９】
【数３３】
2010107314.tif 000035
【００９０】
【数３４】
2010107314.tif 000036
【００９１】
粒子数密度は、下式（３３）となる。
【数３５】
2010107314.tif 000037
【００９２】
したがって、時間ステップｎ＋１のものは、ｎ０と一致する。
液相の計算を行うにあたり、時間刻みΔtは、下式（３４）を満たすように設定する。
【００９３】
【数３６】
2010107314.tif 000038
【００９４】
ここで、ｚ０、Ｃｒｍａｘおよびｕｍａｘは、初期粒子間距離、クーラン数の上限値および流体速度の最大値である。
【００９５】
［固相］
固相の時間ステップｎにおける流れの計算については、液相の時間ステップｎにおける流れの計算と対応している（ステップＳＴ３１〜ステップＳＴ３９）。
ＤＥＭの計算プロセスについて述べる。固相の運動方程式を計算するにあたり、式（４）で示した流体力、粒子-粒子間接触力および外力に基づいた加速度および角加速度を求め、粒子の仮の速度、仮の位置、角速度および角度変位を、それぞれ、下式（３５）〜下式（３８）のように陽的に計算する。
【００９６】
【数３７】
2010107314.tif 000039
【００９７】
【数３８】
2010107314.tif 000040
【００９８】
【数３９】
2010107314.tif 000041
【００９９】
【数４０】
2010107314.tif 000042
【０１００】
ただし、太字ｕｓ＊およびｘｓ＊は、式（３１）および式（３２）と同様に、液相の第２段階の結果を使って、下式（３９）および下式（４０）のように修正される。
【０１０１】
【数４１】
2010107314.tif 000043
【０１０２】
【数４２】
2010107314.tif 000044
【０１０３】
ＤＥＭの粗視化モデルにおけるΔtは、下式（４１）となり、オリジナル粒子体系のものと同じになる。安定に解析できるように時間刻みを設定する必要がある。
【０１０４】
【数４３】
2010107314.tif 000045
【０１０５】
続いて、図１１および図１２を参照して本手法を適用したシミュレーション結果について説明する。
図１１は、本発明の実施の形態における固液混相流の解析体系の初期配置を示す図である。図１２は、本発明の実施の形態における固体粒子の凝集の状態を示すシミュレーション結果である。なお、図１１および図１２において、符号２０は周期境界線を、符号２１は粗子化した固体粒子を、符号２２はオリジナルの固体粒子を示す。
図１１に示すように初期配置された固体粒子に対して、液相の流れを矢印の方向（図１１おいて紙面左側から右側に向かう方向）与えると、固体粒子は図１２（ｂ）に示すように凝集することが確認された。したがって本手法を適用することにより、図１２（ｂ）に示す粗子化した固体粒子の凝集状態は、図１２（ａ）に示すオリジナルの固体粒子の凝集状態を定性的に模擬することが可能となる。
【０１０６】
したがって、上述の本実施形態によれば、液相であるガラスと固相である複数の白金族粒子とからなる固液混相流の状態を示す解析対象データに離散要素法に基づく演算処理をコンピュータを用いて施すことにより上記複数の白金族粒子の挙動を解析するシミュレーション方法であって、上記白金族粒子のいくつかで構成される群を粒子径が上記白金族粒子の粒子径よりも大きな粗視化粒子としてモデル化する粗視化処理を上記解析対象データに施す粗視化モデル生成工程を有し、該粗視化モデル生成工程によって得られた粗視化モデルデータに上記離散要素法に基づく演算処理を施すという構成を採用することによって、白金族粒子の群をモデル化し、計算上大きな粗視化粒子を用いて白金族粒子体系を模擬するため、実際よりも少ない粒子数でガラス溶融炉内の白金族粒子の挙動を解析することが可能となる。
したがって、本実施形態では、現実的な時間でガラス溶融炉内の全ての白金族粒子の挙動を模擬することができるシミュレーション方法を提供できる効果がある。
【０１０７】
また、本実施形態は、上記粗視化モデル生成工程で、上記固液混相流において粒子径の大きさに依存して作用する流体力、接触力、付着力（ファンデルワールス力）および斥力（電気二重層力）を、粒子径の大きな上記粗視化粒子と上記粗視化粒子を構成する粒子径の小さな複数の上記白金族粒子との間で一致させるモデル化を施すことによって、粒子径の大きさに依存して作用する所定の力を白金族粒子の群と粗視化粒子との間でバランスをとることで、実際の白金族粒子の挙動を粗視化粒子で精度良く模擬することが可能となる。
【０１０８】
以上、図面を参照しながら本発明の好適な実施形態について説明したが、本発明は上記実施形態に限定されるものではない。上述した実施形態において示した各構成や組み合わせ等は一例であって、本発明の主旨から逸脱しない範囲において設計要求等に基づき種々変更可能である。
【０１０９】
例えば、上述した本実施形態では、本粗視化モデルを、固相をラグランジュ的手法で、液相をラグランジュ的手法で計算するラグランジュ・ラグランジュ的手法に適用したが、液相をオイラー的手法で計算するオイラー・ラグランジュ的手法にも適用することができる。
また、例えば、上述した本実施形態では、液相にＭＰＳ法を適用したが、液相にＳＰＨ法、有限差分法、有限体積法等を適用してもよい。
【０１１０】
また、例えば、上述した本実施形態では、粗視化モデル生成工程で、上記固液混相流において粒子径の大きさに依存して作用する流体力、接触力、付着力（ファンデルワールス力）および斥力（電気二重層力）を、粒子径の大きな上記粗視化粒子と上記粗視化粒子を構成する粒子径の小さな複数の上記白金族粒子との間で一致させるモデル化を施すと説明したが、ここにさらにランダム力を加える構成であってもよい。ランダム力（Ｒ２）は、下式（４２）のように表される。
【０１１１】
【数４４】
2010107314.tif 000046
【０１１２】
ここで、ξは、ストークス抵抗係数である。
【図面の簡単な説明】
【０１１３】
【図１】本発明の実施の形態におけるシミュレーション装置の概略構成を示すブロック図である。
【図２】本発明の実施の形態におけるシミュレーション方法を示すフローチャートである。
【図３】本発明の実施の形態における粒子径の大きな粗視化粒子と該粒子を構成する粒子径の小さなオリジナル粒子群との挙動の違いを示す図である。
【図４】本発明の実施の形態における粒子径の大きな粗視化粒子と該粒子を構成する粒子径の小さなオリジナル粒子群との挙動の違いを示す図である。
【図５】本発明の実施の形態における粒子径の大きな粗視化粒子と該粒子を構成する粒子径の小さなオリジナル粒子群との挙動の違いを示す図である。
【図６】本発明の実施の形態における固体粒子に作用する接触力のモデルを示す図である。
【図７】本発明の実施の形態における固体粒子の並進運動において粗視化粒子のモデルとオリジナル粒子群のモデルとを示す図である。
【図８】本発明の実施の形態における固体粒子の回転運動において粗視化粒子のモデルとオリジナル粒子群のモデルとを示す図である。
【図９】本発明の実施の形態における固体粒子の二体衝突において粗視化粒子のモデルとオリジナル粒子群のモデルとを示す図である。
【図１０】本発明の実施の形態におけるシミュレーション方法の計算アルゴリズムを示す図である。
【図１１】本発明の実施の形態における固液混相流の解析体系の初期配置を示す図である。
【図１２】本発明の実施の形態における固体粒子の凝集の状態を示すシミュレーション結果である。
【符号の説明】
【０１１４】
１…シミュレーション装置、ＳＴ１…ステップ（粗視化モデル生成工程）
【図１】
2010107314.tif 000051
【図２】
2010107314.tif 000052
【図３】
2010107314.tif 000053
【図４】
2010107314.tif 000054
【図５】
2010107314.tif 000055
【図６】
2010107314.tif 000056
【図７】
2010107314.tif 000057
【図８】
2010107314.tif 000058
【図９】
2010107314.tif 000059
【図１０】
2010107314.tif 000060
【図１１】
2010107314.tif 000061
【図１２】
2010107314.tif 000062