JP2011191848A - 拡散現象解析方法、拡散現象解析装置およびプログラム - Google Patents

Abstract

【課題】従来等方的な拡散現象を解析する用途のみに用いていた格子ボルツマン法を、非等方的な拡散現象に対しても適用できるようにする。
【解決手段】格子ボルツマン法における衝突式に含まれる衝突行列に、非等方運動を反映させる成分を与える。これにより、等方的な拡散現象はもとより、非等方的な拡散現象を解析することが可能となる。
【選択図】図９

Description

本発明は、格子ボルツマン法を用いた拡散現象解析方法に関する。

拡散方程式やナビエ−ストークス（Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ）方程式などの拡散現象に関する運動方程式を解くことにより拡散状態を求める拡散現象解析が従来から行われている。これらの運動方程式の解法の一つとして、空間、時間を離散化し、離散化された空間、時間における仮想粒子と呼ばれる粒子の挙動を求めることにより、拡散状態（例えば流体及び固体内の物質の拡散状態や熱拡散状態など）を求める格子ボルツマン法（Ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ Ｍｅｔｈｏｄ、ＬＢＭ）と呼ばれる手法が知られている。

空間の離散化は、所定の解析空間を複数の格子点で表すことにより行われる。また、時間の離散化は、所定の解析期間を複数の時間ステップに分けることにより行われる。

格子ボルツマン法による拡散現象の解析は、仮想粒子の分布を求めることにより行われる。格子ボルツマン法では、仮想粒子は並進過程と衝突過程を繰り返し、衝突のたびに非平衡状態の粒子が減少するものと定めている。ここで、非平衡状態とは所定の速度方向成分を持った粒子が格子点を通過するときに、格子点に入るときと出たときとで当該速度方向成分を持った粒子の数が異なる（収支が一致しない）状態を指している。したがって、衝突のたびに非平衡状態の粒子が減少するとは、衝突のたびに粒子の収支の差が小さくなることを指している。

格子ボルツマン法には並進過程における仮想粒子の分布を求める並進式と、衝突過程における仮想粒子の分布を求める衝突式とが設けられており、さらに衝突式には粒子同士の衝突状態を規定する衝突行列が含まれている。拡散現象を表わす運動方程式を解く際には、当該運動方程式に対応する値を予め並進式、衝突式の各パラメータに入力する。

各パラメータを入力した並進式と衝突式とを演算することにより、並進過程及び衝突過程における仮想粒子の分布を求めることができる。並進式と衝突式を演算する処理を、全時間ステップについて行うことにより、解析期間における仮想粒子分布の時間変化を求めることができる。さらに、この粒子分布の時間変化を疎視化することで拡散状態を把握することができる。

ここで、拡散現象（例えば流体の流れや熱の拡散）は着目する方向によって物理量（移動量やフラックス等）が変わらない等方的な運動と、着目する方向によって物理量が変化する非等方的な運動とに大別することができる。例えば、所定の場（ｆｉｅｌｄ）の勾配や湧き出しが中心点から同心円状に広がっていれば等方的な場であり、楕円など非同心円状に広がっていれば非等方的な場である。

一般的に等方的な運動はその対象性から非等方的な運動に比べて計算が容易であることが知られている。格子ボルツマン法は格子点の数×時間刻み分の演算を必要とし、過大な計算負荷がかかることから、従来は格子ボルツマン法による解析対象を計算負荷の小さい等方的な運動に限定し、衝突式および並進式も等方的な運動に対応するように各パラメータが定められている。

ドミニク デュミエール、イリーナ ギンズバーグ、マンフレッド クラフィチェク、ピエール ラレマン、リシ ルオー「三次元空間における多緩和時間格子ボルツマンモデル（Ｍｕｌｔｉｐｌｅ−ｒｅｌａｘａｔｉｏｎ−ｔｉｍｅ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｍｏｄｅｌｓ ｉｎ ｔｈｒｅｅ ｄｉｍｅｎｓｉｏｎｓ）」、ロイヤル ソサエティー（Ｔｈｅ Ｒｏｙａｌ Ｓｏｃｉｅｔｙ）、（英国）、平成１４年、第３６０巻、ｐ．４３７−４５１ ルイ デュー、バオチェン シ、シンワン チェン「非圧縮性流体流れのための多緩和時間格子ボルツマンモデル（Ｍｕｌｔｉｐｌｅ−ｒｅｌａｘａｔｉｏｎ−ｔｉｍｅ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｍｏｄｅｌ ｆｏｒ ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ ｆｌｏｗ）」、フィジクス レターズ（Ｐｈｙｓｉｃｓ Ｌｅｔｔｅｒｓ）、（オランダ）平成１８年、セクションＡ、第３５９巻、ｐ．５６４−５７２

ところで、時間ステップの間隔を従来よりも長くすることのできる手法が近年提案されている。例えば非特許文献１、２には、仮想粒子の分布を表す分布関数ベクトルをモーメントベクトルと呼ばれるベクトルに変換することで、時間ステップの間隔を従来よりも長くすることができることが記載されている。時間ステップの間隔が長くできるということは、解析期間あたりの時間ステップ数を減らすことができるということになる。時間ステップ数を低減できる分、計算負荷は従来よりも軽減される。この計算負荷の軽減に伴い、等方的な拡散現象よりも計算負荷の高い非等方的な拡散現象に対しても格子ボルツマン法による拡散現象の解析を行いたいという要望が高まっている。

そこで、本発明は非等方的な拡散現象の解析を可能とするように格子ボルツマン法を改良するとともに、当該格子ボルツマン法を用いた拡散現象解析方法を提供することを目的とする。

本願発明は、並進過程及び衝突過程で表わされた粒子の挙動を算出することにより拡散現象を解析する格子ボルツマン法を用いた、拡散現象解析方法に関するものである。衝突過程を表す衝突式には、粒子の衝突状態を定める衝突行列が含まれている。本願発明に係る拡散現象解析方法においては、この衝突行列に非等方成分が含まれた衝突式を用いて拡散現象の解析を行う。

また、上記発明において、解析対象となる拡散現象は拡散方程式によって表わされ、非等方成分は、拡散方程式における拡散係数テンソルの成分に応じて定められることが望ましい。

また、上記発明において、衝突式は、粒子の分布を表わす分布関数ベクトルと、分布関数ベクトルをモーメントベクトルに変換する変換行列とを含み、モーメントベクトルはフラックス成分を含んでいる。衝突行列の行列成分のうち、モーメントベクトルのフラックス成分に対応する行列成分に非等方成分が含まれることが望ましい。

また、本願発明は、並進過程及び衝突過程で表わされた粒子の挙動を算出することにより拡散現象を解析する格子ボルツマン法を用いて拡散現象解析を行う、拡散現象解析装置に関するものである。衝突過程を表す衝突式は、粒子の衝突状態を定める衝突行列を有している。拡散現象解析装置の構成のうち、格子ボルツマン法の演算を行う演算部は、衝突行列に非等方成分が含まれた衝突式を用いて拡散現象の解析を行う。

また、本願発明は、並進過程及び衝突過程で表わされた粒子の挙動を算出することにより拡散現象を解析する格子ボルツマン法を用いて拡散現象解析を行う拡散現象解析装置に、拡散現象解析を実行させるプログラムに関するものである。衝突過程を表す衝突式は、粒子の衝突状態を定める衝突行列を有している。プログラムは、衝突行列に非等方成分が含まれた衝突式を用いて、拡散現象解析装置に拡散現象解析を実行させる。

本願発明者らは、格子ボルツマン法における衝突式に含まれる衝突行列に非等方運動を反映させる成分を与えることにより、非等方的な拡散現象を演算できるという知見を得た。非等方成分を含む衝突行列を備えた衝突式を用いることで、従来は解析できなかった非等方的な拡散現象の解析が可能となる。

本実施形態に係る拡散現象解析装置を例示する図である。 本実施形態に係る拡散現象解析のフローチャートを示す図である。 解析対象の境界条件について説明する図である。 本実施形態に係る拡散現象解析装置の機能ブロック図である。 ３次元の格子空間を例示する図である。 本実施形態に係る拡散現象解析を行った結果を示す図である。 本実施形態に係る拡散現象解析を行った結果を示す図である。 本実施形態に係る拡散現象解析を行った結果を示す図である。 本実施形態に係る拡散現象解析を行った結果を示す図である。 本実施形態に係る拡散現象解析の結果と、厳密解とを比較した図である。

本発明の実施の形態における拡散現象解析装置の構成を図１に示す。拡散現象解析装置１０は、図１に示すように、演算部１２、入力部１４、表示部１６及び記憶部１８を含んで構成される。

演算部１２は、いわゆるＣＰＵであり、記憶部１８に予め記憶させてある拡散現象解析プログラムを実行することによって後述する拡散現象解析処理を行う。入力部１４は、キーボード等の入力手段を含んで構成され、計算空間である格子空間における格子点の設定、格子ボルツマン法における並進式や衝突式のパラメータの設定、全格子点の仮想粒子分布の初期値、境界条件、時間刻みΔｔ、終了時刻Ｔ等の入力をユーザから受け付ける。入力された情報は記憶部１８に格納及び保持される。表示部１６は、ディスプレイやプリンタ等の情報出力手段を含んで構成され、拡散現象解析処理の対象となる格子空間、設定されたパラメータ、処理経過及び処理結果等の情報を表示出力する。記憶部１８は、半導体メモリやハードディスク等の記憶手段を含んで構成され、拡散現象解析プログラム、拡散現象解析処理のために設定された情報、解析処理での演算結果等を格納及び保持する。記憶部１８は、演算部１２から適宜アクセス可能である。

次に、演算部１２が行う拡散現象解析処理について説明する。演算部１２は格子ボルツマン法による数値計算を行うことにより、拡散現象解析を行う。まず格子ボルツマン法について概説する。格子ボルツマン法においては、仮想粒子の挙動は格子点から他の格子点に移動する並進過程と、仮想粒子同士が衝突する衝突過程の２つの過程で表される。仮想粒子は並進過程と衝突過程を繰り返すと仮定する。また、並進過程において、時間刻みΔｔごとに仮想粒子は必ず格子点上に存在し、格子線の途中でとどまることはできないものと仮定する。また、衝突は瞬時に行われ、かつ格子点上でのみ起こるものと仮定する。

さらに、並進過程と衝突過程における仮想粒子の分布を示す分布関数をそれぞれ定める。この分布関数を用いて、所定の格子点について、並進過程における仮想粒子の分布と衝突過程における仮想粒子の分布を時間刻みΔｔごとに算出する。

なお、格子ボルツマン法においては、衝突のたびに非平衡状態の粒子が減少するものと定めている。このような衝突のルールを反映させるために、衝突過程を表す衝突式には粒子の衝突状態を規定する衝突行列が含まれている。

並進過程および衝突過程における仮想粒子の分布の算出を、定められた範囲の格子点（例えば全格子点）について行う。さらに定められた範囲の格子点についての演算を初期時刻ｔ＝０から終了時刻ｔ＝Ｔまで行う。最終的に得られた仮想粒子の分布の空間的、時間的な変化は疎視化すると拡散状態を表わすものとして捉えることができ、これにより拡散現象（流れの変化や、固体内における物質拡散、熱拡散等）を把握することが可能となる。

演算部１２は上述した格子ボルツマン法を用いて拡散現象解析を行う。演算部１２が行う演算処理のフローチャートを図２を用いて説明する。まず解析対象となる拡散現象の初期状態を入力部１４から入力する（Ｓ１）。初期状態とは初期時刻ｔ＝０における計算空間（例えば全格子点）の場の状態を表す関数を示す。

次に、演算部１２は衝突過程の演算を行い（Ｓ２）、続いて並進過程の演算を行う（Ｓ３）。並進過程の演算後、演算部１２は境界条件の演算を行う（Ｓ４）境界条件とは、図３に示すような計算領域Ωにおける境界部分∂Ω_d、∂Ω_nにおける仮想粒子の挙動を定めるものである。この境界部分においては、例えば仮想粒子が境界部分に当たって跳ね返るなど、通常の（つまり、境界から十分に離れた領域における）並進過程とは異なる振る舞いをするように予め記憶部１８に境界条件の数式が記憶されている。

さらに演算部１２は演算した時刻における各格子点の仮想粒子の分布を算出する（Ｓ５）。その後演算部１２は演算を行った時刻が終了時刻Ｔであるか否かを判定する（Ｓ６）。もし終了時刻Ｔであれば演算を終了する（Ｓ７）。この演算結果は表示部１６に出力され、入力者やその他のユーザに対して拡散現象解析の結果を表示する。もしまだ終了時刻Ｔに達していない場合には時間刻みΔｔを一つ増やして（時刻を進めて）再び衝突過程の演算を行う（Ｓ２）。

以上説明した拡散現象解析処理を実行する拡散現象解析装置１０の機能ブロック図を図４に示す。拡散現象解析装置１０は、Ｓ１を実行する入力データ取得手段２０、Ｓ２を実行する衝突過程演算手段２２、Ｓ３を実行する並進過程演算手段２４、Ｓ４を実行する境界条件演算手段２６、Ｓ５を実行する物理量計算手段、Ｓ６を実行する終了時刻確認手段２８、Ｓ７を実行する解析結果出力手段３０を含んで構成される。

以上、演算部１２が行う拡散現象解析処理をフローチャートを用いて概説した。次に、演算部１２が行う演算の内容を詳細に説明する。この演算内容と上述したフローチャートの各工程との対応関係については演算内容の説明後に示す。まず、拡散現象解析の対象となる拡散方程式を下記数式１に例示する。

ここで、φは求めようとする未知数、例えば温度、濃度など、拡散方程式に支配される物理量を示す。また、ｔは時間、ｘは位置座標、ｖは既知の流速ベクトルを示し、例えば空気の流れや、溶媒の流れを表している。さらにＤ_ijは拡散係数テンソルを表わしている。なお、数式１では表記を省略したが、例えば化学反応による発熱や物質生成を表わすために右辺に生成項ベクトルＲを加えても良い。

上記の数式１を格子ボルツマン法によって解く。格子ボルツマン方程式を数式２に示す。

ここで、ｔは時刻、Δｔは時間刻み、ｘは格子点位置、ｅαはα方向の単位ベクトル、Δｘは格子点間隔、ｆα（ｔ，ｘ）は分布関数を表わしている。分布関数とは仮想粒子の数密度を表わし、任意の実数値を取り得る。

また、分布関数ｆαのサフィックスαは仮想粒子の速度方向を表わし、速度方向の単位ベクトルｅαに対応している。このサフィックスαについて説明する。計算空間である格子空間が図５に示すような三次元の場合、仮想粒子の速度方向ベクトルｅαは下記数式３で表わされる。この数式３から、αは０〜６の値を取りうる。

さらに、数式２の右辺第２項のχαは衝突演算子を表わしている。これは、仮想粒子同士が衝突することにより、衝突前と衝突後とで速度方向が変わる（静止も含む）ことを考慮したものである。すなわち衝突演算子χαは仮想粒子の衝突状態を規定している。所定の速度方向成分αについて（例えばα＝１）、右辺第２項は、（１）衝突によって非α方向の仮想粒子がα方向の粒子に変わった場合と、（２）衝突によってα方向の仮想粒子が非α方向の粒子に変わった場合とを考慮している。

すなわち、数式２は、時刻ｔ、位置ｘにおける速度方向αを持った粒子が、Δｔ後には位置ｘ＋ｅαΔｘに移動することを表わしている。さらに、速度方向αを持った粒子について、右辺第２項では、（１）粒子同士の衝突により速度方向がαになったもの、あるいは、（２）αでなくなったものの粒子数を反映している。この（１）（２）をより明確にするために、数式２は下記数式４のように書き換えることができる。ただし、αの取り得る値は０〜ｋ（ｋは整数）の整数値としている。

数式４は下記数式５のように表すことができる。

数式５において、Ａαβは衝突演算子χをまとめた衝突行列を表わしている。また、右辺第２項は仮想粒子の平衡、非平衡の概念を導入することによりさらに書き換えることができる。平衡状態とは、ある速度成分の粒子がある格子点に入ったときと出たときとで数が変わらない（収支が一致する）状態を指している。平衡状態の分布関数をｆ^eq（ｔ，ｘ）で表わす。また、ある速度成分の粒子が格子点から出たときに粒子数が入ったときよりも減っている又は増えている（収支が一致しない）状態を非平衡状態と呼び、この収支の差分を非平衡分布関数ｆ^neq（ｔ，ｘ）で表わす。数式５は、分布関数ｆ^neq（ｔ，ｘ）を用いて数式６のように書き直すことができる。

なお、記号ＢはＡでない新たな衝突行列を示している。さらに非平衡状態の分布関数ｆ^neq（ｔ，ｘ）は平衡状態の分布関数ｆ^eq（ｔ，ｘ）を用いて下記数式で表わすと、下記数式７のように表すことができる。

また、数式６の右辺第２項について、仮想粒子の衝突により非αの速度成分を持った粒子がαに変わる粒子数よりも、αの速度成分を持った粒子が非αに変わる粒子数の方が多い場合、つまり衝突によりαの粒子数が減る場合、数式６の右辺第２項に負符号を与えると理解が容易になる。このことから、数式６は下記数式８のように変形することができる

次に、数式８を衝突過程を表わす衝突式と並進過程を表わす並進式とに分ける。衝突後の分布関数を、分布関数ｆの頭部にハット記号（＾）が付された関数で表わすとき、衝突式は下記数式９で表わされ、並進式は下記数式１０で表わされる。

なお、数式２〜６、８、９においては表記を省略したが、それぞれの右辺第３項に定数項や変数項を設けても良い。例えば数式１を説明する際に付言した生成項ベクトルＲを数式に反映させるために、それぞれの右辺第３項としてΔｔＲωαを加えても良い。係数ωαについては後述する。

さらに衝突式（数式９）について説明する。数式９を展開すると下記数式１１のようになる。ただし、格子空間を三次元とする。

ここで、本実施形態においては、分布関数ベクトルＦβ＝（ｆ₀ ^eq（ｔ，ｘ）−ｆ₀（ｔ，ｘ），…，ｆ₆ ^eq（ｔ，ｘ）−ｆ₆（ｔ，ｘ））^Tをモーメント空間を呼ばれる空間に写像する。分布関数ベクトルＦβをモーメント空間に写像し、モーメントベクトルｍβと呼ばれるベクトルに変換することで、非線形な拡散現象を表現することが可能となる。モーメント空間に写像する行列を変換行列Ｍで表わすと、数式９は下記数式１２のように表わすことができる。

ここで、Ｓは衝突行列であり、数式９の行列Ｂに対応し、粒子の衝突状態を規定するものである。さらに、モーメント空間に写像されたモーメントベクトルｍβに衝突行列Ｓを掛けた後に、当該モーメントベクトルｍβをもとの分布関数ベクトルＦβに変換するために変換行列Ｍの逆行列Ｍ^-1が設けられている。

格子空間が三次元である場合、変換行列Ｍは下記数式１３のように表わすことができる。

ここで、ａ、ｂ、ｃ、ｄ、ｇ_n、ｈ_n、ｌ_n（ｎ＝０〜６）は予め定めた定数である。なお、ｇ_n、ｈ_n、ｌ_nは、変換行列Ｍ中の行ベクトルがすべて直交するように定める。

変換行列Ｍの２〜４行目の成分について説明する。２行目の成分は２列目にｂ、３列目に成分−ｂが定められ、それ以外の成分は０となっている。一方、変換行列Ｍにより変換される分布関数ベクトルＦβの各成分に着目すると、ｆ₁の成分が変換行列Ｍの２行目２列目の成分ｂに対応し、またｆ₂の成分が変換行列Ｍの２行目３列目の成分−ｂに対応する。したがって、変換行列Ｍにより変換されたモーメントベクトルｍβの第２成分は（ｂｆ₁−ｂｆ₂）となる。さらに、数式３から、ｆ₁およびｆ₂はｘ軸方向の単位ベクトルｅ₁およびｅ₂の速度方向を有する粒子に対応する。つまり、（ｂｆ₁−ｂｆ₂）はｘ軸方向のフラックス（流束）を表わしていることが理解される。このことから、変換行列Ｍの２行目成分によって変換されたモーメントベクトルｍβの第２成分ｍ₂は、ｘ軸方向のフラックス成分を意味していることが理解される。

また、変換行列Ｍの３行目の成分は４列目にｃ、５列目に−ｃが定められ、分布関数ベクトルＦβのｙ軸方向の移動を表わすｆ₃のｆ₄に対応している。このことから、変換行列Ｍの３行目成分によって変換されたモーメントベクトルｍβの第３成分ｍ₃は、ｙ軸方向のフラックス成分を意味していることが理解される。さらに、変換行列Ｍの４行目の成分は、６列目にｄ、７列目に−ｄが定められ、分布関数ベクトルＦβのｚ軸方向の移動を表わすｆ₅、ｆ₆に対応している。このことから、変換行列Ｍの４行目成分によって変換されたモーメントベクトルｍβの第４成分ｍ₄は、ｚ軸方向のフラックス成分を意味していることが理解される。変換行列Ｍの各成分に数値が与えられた例を下記数式１４に示す。

さらに数式１２に則り、モーメントベクトルｍβに衝突行列Ｓを掛ける。本発明者らは、この衝突行列Ｓに非等方成分を含めることにより、非等方な拡散現象を解析できるという知見を得た。

本実施形態に掛かる三次元格子空間の衝突行列Ｓを下記数式１５に示す。

ここで、各γの値は衝突行列Ｓの逆行列Ｓ^-1に基づいて定められる。逆行列Ｓ^-1を下記数式１６に示す。

ここで、係数τ₀、τ₄、τ₅、τ₆の値は予め定められる。例えば０．５以上の値を取り、１であると好適である。また、係数τ_ijの値は下記数式１７に示されるように数式１の拡散係数テンソルＤ_ijの値に応じて定められる。

ここで、δ_ijはクロネッカーのデルタ（ｉ＝ｊのとき１、ｉ≠ｊのとき０、）を表わしている。また、係数εは下記数式１８によって定められる。

ここで、α_maxは分布関数ｆαおよび単位ベクトルｅαのサフィックスαが取りうる最大値を示す。例えば格子空間が三次元であればα_max＝６である。衝突行列Ｓとモーメントベクトルｍβの積を下記数式１９に示す。なお、モーメントベクトルｍβの第２から第４成分がそれぞれｘ方向、ｙ方向、ｚ方向のフラックスを表わしていることから、これらの成分をｍ_x、ｍ_y、ｍ_zと記載する。

数式１９の右辺を参照すると、ｘ方向のフラックスを表わす２番目の成分にｙ方向およびｚ方向のモーメント成分が含まれていることが理解される。以下同様に、ｙ方向のフラックスを表わす３番目の成分、ｚ方向のフラックスを表わす４番目の成分に異方向のモーメント成分が含まれていることが理解される。すなわち、衝突行列Ｓの非対角成分に非ゼロ値の係数を入れることにより、所定方向のフラックスに対して異方向のフラックスの影響が反映される。この数式１９から、衝突行列Ｓに非等方成分を入れることによって、非等方的な拡散現象を表現できることが理解される。

さらに、衝突行列Ｓの対角成分に着目すると、数１９よりγ_xx、γ_yy、γ_zzはそれぞれ異なる値を取りうる。このことから、γ_xx、γ_yy、γ_zzにそれぞれ異なる値を入れることにより、非対角成分を０にしても、各方向でフラックスの異なる対角非等方な運動を表現することができる。

数式１２に戻り、Ｓｍβに対して変換行列Ｍの逆行列Ｍ^-1が掛けられる。これによりモーメント空間の各点はもとの空間（速度空間と呼ぶ）に写像され、モーメントベクトルｍαは分布関数ベクトルＦβに変換される。さらに数式１２の右辺第１項における分布関数ｆαと変換後の分布関数ベクトルＦβとの差を取ることにより、衝突過程の分布関数を求めることができる。

次に、上述した演算過程と、図２のフローチャートとの対応関係について説明する。演算部１２は、入力部１４から入力され、または記憶部１８に予め記憶された初期条件に基づいて初期分布を作成する。ここで、解析対象が数式１の拡散方程式である場合、初期条件は下記数式２０で表わされる。

ここで、ψは初期時刻における場（ｆｉｅｌｄ）を表わす関数である。この初期分布を格子ボルツマン法の式に当てはめる。当てはめた数式の例を下記数式２１に示す。なお、初期分布の式は数式２１に限られず、種々の形態にて表現することが可能である。

ここで、ωαは格子点空間の次元により異なる値を取る定数であり、格子点空間が三次元の場合、ωαは下記数式２２により定められる。

σは数式１８におけるものと同一のものである。次に、演算部１２は数式１２で示された衝突式に基づき、衝突過程の演算を行う。ここで、数式１２の平衡状態における分布関数ｆ^eqを下記数式２３に例示する。

次に、演算部２１は数式１０で示された並進式に基づき、並進過程の演算を行う。続いて、演算部２１は境界条件（Ｓ４）の演算を行う。

境界条件の演算について、本実施形態においては図３における境界部分にある仮想粒子の挙動を定めている。まず、図３の境界∂Ω_dの境界条件を下記数式２４に示す。

これを並進式（数式１０）に当てはめると下記数式２５のようになる。

ここで、分布関数ｆのサフィックスβはαに対する反対方向を表わしている。また、図３における境界∂Ω_nの境界条件を下記数式２６に示す。

これを並進式に当てはめると下記数式２７のようになる。

ここで、分布関数ｆのサフィックスβはαに対する反対方向を表わしている。境界条件が演算されると、演算部１２は物理量の計算を行う（Ｓ５）。本実施形態においては時間ｔ＋Δｔにおける各格子点のφの値を算出する。φの算出式を下記数式２７に示す。

次に、現在演算処理を行っている時刻が終了時刻Ｔに達しているか否かを判定し（Ｓ６）、達していない場合には時間刻みΔｔを加えて再び衝突過程の演算を行う（Ｓ２）。終了時刻Ｔに達している場合は演算処理を終了する（Ｓ７）。

以上説明した拡散現象解析処理の結果を図６から９に示す。まず、解析対象について説明する。数式１において、流速ベクトルｖはｘ軸方向を向いているものとする。すなわち、ｖ＝（１０，０，０）。また、生成項ベクトルＲは０とする。また、解析対象の初期分布ψを下記数式２９に示す。

ここで、φ₀＝０．０１、σ₀ ²＝０．０２とする。初期時刻におけるφの分布を図６に示す。なお、図６は格子空間のうち、ｙ＝０におけるｘ−ｚ平面を示している。以下、図７〜９も同様にｙ＝０におけるｘ−ｚ平面を示している。

また、格子ボルツマン法の格子空間について、ｘ？［−１／２，３／２］、ｙ，ｚ？［−１，１］とした。また、境界面では周期境界条件を課した。また、時間刻みΔｔについて、Δｔ＝０．０００５×（３２／Ｎ）²とした。さらに拡散係数テンソルＤ_ijは下記の３種類を用意した。

ここで、数式３０は拡散現象が等法的（Ｉｓｏｔｒｏｐｉｃ）な場合を示している。また、数式３１は対角非等方（Ｄｉａｇｏｎａｌ ａｎｉｓｏｔｒｏｐｉｃ）な場合、すなわちｘ、ｙ、ｚ方向のそれぞれについてフラックスが異なる場合を示している。数式３１においてはｚ方向のフラックスが最も大きく、ｘ方向のフラックスが最も小さい。また、数式３２は非対角非等方（Ｏｆｆ−ｄｉａｇｏｎａｌ ａｎｉｓｏｔｒｏｐｉｃ）な場合を示している。つまり、所定方向（例えばｘ方向）のフラックスが他の方向（例えばｙ方向、ｚ方向）の影響を受ける場合を示している。上述したように、数式１９のτ_ij（ｉ＝１、２、３，ｊ＝１、２、３）は数式３０〜３２の拡散係数テンソルＤ_ijの各成分に基づいて定められる。また、数式１９のτ₀、τ₄、τ₅、τ₆はいずれも１とした。

数式３０の拡散係数テンソルＤ_ijを適用した際の、本実施形態に掛かる拡散現象解析の結果を図７に示す。ここで、図７はｔ＝０．０２５におけるφの分布を示している。以下、図８、図９ともｔ＝０．０２５におけるφの分布を示す。図７の等方的な分布から、本実施形態に掛かる拡散現象解析方法は、少なくとも定性的には与えられた拡散方程式を適正に解いていることが理解される。

数式３１の拡散係数テンソルＤ_ijを適用した際の、本実施形態に掛かる拡散現象解析の結果を図８に示す。上述したように、数式３１の拡散係数テンソルＤ_ijはｚ方向のフラックスが大きく、ｘ方向のフラックスが小さい。一方、図８におけるφの分布はｘ方向には延びず、ｚ方向に延びる分布を示している。この分布は数式３１の拡散係数テンソルＤ_ijのｘ方向のフラックスおよびｚ方向のフラックスの傾向と一致する。このことから、本実施形態に掛かる拡散現象解析は、数式３１の拡散係数テンソルＤ_ijを適用したときの拡散方程式を適正に解いており、対角非等方な拡散現象であっても解析が可能であることが理解される。

数式３２の拡散係数テンソルＤ_ijを適用した際の、本実施形態に掛かる拡散現象解析の結果を図９に示す。上述したように、数式３２の拡散係数テンソルＤ_ijは所定方向のフラックスは他方向のフラックスの影響を受ける。一方、図９におけるφの分布は斜め方向に延びる分布を示している。この分布は数式３２の拡散係数テンソルＤ_ijのフラックスの傾向と一致する。このことから、本実施形態に掛かる拡散現象解析は、数式３２の拡散係数テンソルＤ_ijを適用したときの拡散方程式を適正に解いており、非対角非等方な拡散現象であっても解析が可能であることが理解される。

なお、本実施例に掛かる拡散現象解析の結果が、解析対象である拡散方程式の解をよく近似していることを図１０を用いて説明する。数式２９から数式３２まで示した拡散方程式には厳密解が存在する。厳密解φ_exactを下記数式３３に示す。

ここで、σ_ij＝σ₀ ²δ_ij―２ｔＤ_ijであり、｜｜σ_ij｜｜はσ_ijの行列式の絶対値を表わしている。また、（σ^-1）_ijはσ_ijの逆行列の（ｉ，ｊ）成分である。この解φ_exactに数式３０、３１、３２のそれぞれの拡散係数テンソルＤ_ijを代入することにより、それぞれの場合におけるφを求めることができる。

厳密解により求めたφと、本実施形態における拡散現象解析により求めたφとの差異を示すグラフを図１０に示す。図１０の縦軸は厳密解により求めたφと本実施形態により求めたφとの絶対誤差の自乗平均値を示している。また、横軸は格子解像度１／Ｎ（Ｎは格子数）を示している。数式３０、３１、３２の拡散係数テンソルＤ_ijを適用したいずれの場合においても誤差の自乗平均値は最大でも１０^-3以下であり、本実施形態による拡散方程式の解法は極めて精度が高く、したがって精度の高い解析を行うことが可能であることが理解される。

なお、上述した実施形態においては、格子空間を３次元としたが、２次元の格子空間に対しても本発明は実施可能である。２次元の格子空間における各種パラメータについて下記に説明する。

２次元格子空間では、仮想粒子の速度方向ベクトルｅαは下記数式３４に示すように５つの単位ベクトルから構成される。

また、変換行列Ｍは下記数式３５のように表すことができる。

ここで、ｄ_n、ｇ_n（ｎ＝０〜４）は変換行列Ｍ中の行ベクトルがすべて直交するように定める。変換行列Ｍの一例を下記数式３６に示す。

さらに、衝突行列Ｓの逆行列Ｓ^-1を下記数式３７に示す。

また、数式２２におけるωαは下記数式３８により定められる。

１０ 拡散現象解析装置、１２ 演算部、１４ 入力部、１６ 表示部、１８ 記憶部、２０ 入力データ取得手段、２２ 衝突過程演算手段、２４ 並進過程演算手段、２６ 境界条件演算手段、２８ 終了時刻確認手段、３０ 解析結果出力手段。

Claims (5)

1. 並進過程及び衝突過程で表わされた粒子の挙動を算出することにより拡散現象を解析する格子ボルツマン法を用いた、拡散現象解析方法であって、
   前記衝突過程を表す衝突式は、前記粒子の衝突状態を定める衝突行列を有し、
   前記衝突行列に非等方成分が含まれた前記衝突式を用いて拡散現象の解析を行うことを特徴とする、拡散現象解析方法。
2. 請求項１に記載の拡散現象解析方法であって、
   解析対象となる前記拡散現象は拡散方程式によって表わされ、
   前記非等方成分は、前記拡散方程式における拡散係数テンソルの成分に応じて定められることを特徴とする、拡散現象解析方法。
3. 請求項１又は２に記載の拡散現象解析方法であって、
   前記衝突式は、前記粒子の分布を表わす分布関数ベクトルと、前記分布関数ベクトルをモーメントベクトルに変換する変換行列とを含み、
   前記モーメントベクトルはフラックス成分を含み、
   前記衝突行列の行列成分のうち、前記モーメントベクトルの前記フラックス成分に対応する行列成分に前記非等方成分が含まれることを特徴とする、拡散現象解析方法。
4. 並進過程及び衝突過程で表わされた粒子の挙動を算出することにより拡散現象を解析する格子ボルツマン法を用いて拡散現象を解析する、拡散現象解析装置であって、
   前記衝突過程を表す衝突式は、前記粒子の衝突状態を定める衝突行列を有し、
   前記格子ボルツマン法の演算を行う演算部は、前記衝突行列に非等方成分が含まれた前記衝突式を用いて拡散現象の解析を行うことを特徴とする、拡散現象解析装置。
5. 並進過程及び衝突過程で表わされた粒子の挙動を算出することにより拡散現象を解析する格子ボルツマン法を用いて拡散現象の解析を行う拡散現象解析装置に、前記拡散現象の解析を実行させるプログラムであって、
   前記衝突過程を表す衝突式は、前記粒子の衝突状態を定める衝突行列を有し、
   前記衝突行列に非等方成分が含まれた前記衝突式を用いて、前記拡散現象解析装置に拡散現象の解析を実行させることを特徴とする、プログラム。