特許第6561233号(P6561233)IP Force 特許公報掲載プロジェクト 2022.1.31 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ ダッソー システムズ シムリア コーポレイションの特許一覧

特許6561233等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子
<>
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000144
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000145
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000146
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000147
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000148
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000149
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000150
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000151
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000152
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000153
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000154
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000155
  • 特許6561233-等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 図000156
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】6561233
(24)【登録日】2019年8月2日
(45)【発行日】2019年8月21日
(54)【発明の名称】等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子
(51)【国際特許分類】
   G16Z 99/00 20190101AFI20190808BHJP
   G06F 17/18 20060101ALI20190808BHJP
【FI】
   G06F19/00 110
   G06F17/18 Z
【請求項の数】45
【全頁数】41
(21)【出願番号】特願2016-530050(P2016-530050)
(86)(22)【出願日】2014年7月24日
(65)【公表番号】特表2016-534436(P2016-534436A)
(43)【公表日】2016年11月4日
(86)【国際出願番号】US2014048004
(87)【国際公開番号】WO2015013507
(87)【国際公開日】20150129
【審査請求日】2017年7月24日
(31)【優先権主張番号】61/858,051
(32)【優先日】2013年7月24日
(33)【優先権主張国】US
(73)【特許権者】
【識別番号】519084700
【氏名又は名称】ダッソー システムズ シムリア コーポレイション
(74)【代理人】
【識別番号】100086771
【弁理士】
【氏名又は名称】西島 孝喜
(74)【代理人】
【識別番号】100088694
【弁理士】
【氏名又は名称】弟子丸 健
(74)【代理人】
【識別番号】100094569
【弁理士】
【氏名又は名称】田中 伸一郎
(74)【代理人】
【識別番号】100067013
【弁理士】
【氏名又は名称】大塚 文昭
(74)【代理人】
【識別番号】100109070
【弁理士】
【氏名又は名称】須田 洋之
(74)【代理人】
【識別番号】100109335
【弁理士】
【氏名又は名称】上杉 浩
(74)【代理人】
【識別番号】100120525
【弁理士】
【氏名又は名称】近藤 直樹
(74)【代理人】
【識別番号】100196612
【弁理士】
【氏名又は名称】鎌田 慎也
(72)【発明者】
【氏名】チェン フードン
(72)【発明者】
【氏名】チャン ラオヤン
(72)【発明者】
【氏名】ゴパラクリシュナン プラディープ
【審査官】 宮地 匡人
(56)【参考文献】
【文献】 特表2010−500654(JP,A)
【文献】 国際公開第2013/085567(WO,A1)
【文献】 米国特許出願公開第2010/0185420(US,A1)
【文献】 米国特許出願公開第2012/0296615(US,A1)
【文献】 米国特許第06915245(US,B1)
(58)【調査した分野】(Int.Cl.,DB名)
G16Z 99/00
G06F 17/18
(57)【特許請求の範囲】
【請求項1】
1又は2以上のコンピュータシステムにより、格子速度セットを用いて流体容積中の粒子の移動をシミュレートするステップであって、前記粒子の移動が前記粒子間の衝突を引き起こすものである、ステップと、
前記1又は2以上のコンピュータにより、前記粒子の前記シミュレートされた移動に基づいて、前記容積内の特定位置における粒子の粒子相対速度を決定するステップであって、前記粒子相対速度が(i)前記容積内の前記特定位置における前記容積のゼロフロー基準粒子絶対速度と(ii)前記容積内の前記特定位置における粒子平均速度との間の差である、ステップ、及び
前記1又は2以上のコンピュータにより、前記粒子相対速度に基づいて、前記衝突を表す非平衡衝突後分布を決定するステップ、により
前記1又は2以上のコンピュータシステムにより、特定の次数は格子速さのセットによりサポートされ、前記格子速度セットの前記格子速さのセットによりサポートされないより高次の項を除外する衝突演算子を含む、前記衝突を表す該特定の次数の非平衡衝突後分布関数を評価するステップと、
を含むコンピュータ実装方法。
【請求項2】
前記格子速度セットは、粒子速度の或る次数まで流体力学運動をサポート、前記シミュレートするステップは、前記1又は2以上のコンピュータシステムによってシミュレートするステップである、請求項1に記載の方法。
【請求項3】
前記格子速度セットに関する前記サポートされる次数は、前記非平衡衝突後分布関数の前記特定の次数より小さくかつこれとは異なり、
前記非平衡衝突後分布関数の前記特定の次数は、前記粒子速度の前記次数によって決定される、請求項2に記載の方法。
【請求項4】
前記容積内の前記特定位置における前記粒子平均速度は、前記特定位置における特定タイプ粒子の粒子平均速度を含む、請求項1に記載の方法。
【請求項5】
前記格子速度セットは、格子ボルツマン法に関連する状態ベクトルのセットである、請求項1に記載の方法。
【請求項6】
前記非平衡衝突後分布関数は、(i)予め定義された物理量に対する非平衡モーメントを維持し、(ii)未定義の物理量に対する非平衡モーメントを前記特定の次数まで除去する、請求項1に記載の方法。
【請求項7】
前記特定の次数は、前記流体速度の格子音速に対する比率に関連した指数値であり、前記格子速度セットは、前記指数値をサポートする、請求項1に記載の方法。
【請求項8】
前記格子速度セットは、格子に限定される空間における運動量状態のセットを含む、請求項1に記載の方法。
【請求項9】
前記非平衡衝突後分布関数は、ガリレイ不変フィルタ処理演算子である、請求項1に記載の方法。
【請求項10】
前記非平衡衝突後分布関数に基づいて、前記流体容積中での前記粒子の衝突過程をモデル化するステップをさらに含む、請求項1に記載の方法。
【請求項11】
前記非平衡衝突後分布関数は、流体力学モーメントに対して1次のサポートを提供する格子速度セットに対するマッハ数に関する1次のガリレイ不変衝突演算子
であり、前記衝突演算子は、
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
iは衝突前の前記粒子の速度ベクトルであり、
u(x,t)は時間tでの特定位置xにおける前記粒子間の平均速度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
iは一定重み係数であり、
は非平衡運動量流束である、請求項1に記載の方法。
【請求項12】
前記非平衡衝突後分布関数は、流体力学モーメントに関して無限次数のサポートを与える、格子速度セットに対する衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
は粒子相対速度であり、
ρは流体密度であり、
は平衡分布関数であり、
は非平衡運動量流束である、請求項1に記載の方法。
【請求項13】
前記非平衡衝突後分布関数は、流体力学モーメントに関して2次のサポートを提供する、格子速度セットに対するマッハ数に関して2次のガリレイ不変衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
iは衝突前の前記粒子の速度ベクトルであり、
u(x,t)は時間tでの特定位置xにおける前記粒子間の平均速度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
iは一定重み係数であり、
は非平衡運動量流束である、請求項1に記載の方法。
【請求項14】
前記予め定義された物理量は、前記特定の容積中の前記流体の質量、前記特定の容積中の前記流体の運動量、及び前記特定の容積中の前記流体のエネルギを含む、請求項に記載の方法。
【請求項15】
前記非平衡衝突後分布関数は、エネルギ流束に関連する衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
は相対的な粒子速度であり、
は平衡分布関数であり、
は非平衡エネルギ流束である、請求項1に記載の方法。
【請求項16】
オペレーションを実行する1又は2以上の処理装置で実行可能な命令を記憶する1又は2以上の機械可読ハードウェア記憶装置であって、前記オペレーションは、
格子速度セットを用いて流体容積中の粒子の移動をシミュレートするステップであって、前記粒子の移動が前記粒子間の衝突を引き起こすものである、ステップと、
前記粒子の前記シミュレートされた移動に基づいて、前記容積内の特定位置における粒子相対速度を決定するステップであって、前記粒子相対速度が(i)前記容積内の前記特定位置における前記容積のゼロフロー基準粒子絶対速度と(ii)前記容積内の前記特定位置における粒子平均速度との間の差である、ステップ、及び
前記粒子相対速度に基づいて、前記衝突を表す非平衡衝突後分布を決定するステップ、により
前記1又は2以上のコンピュータシステムにより、特定の次数は格子速さのセットによりサポートされ、前記格子速度セットの前記格子速さのセットによりサポートされないより高次の項を除外する衝突演算子を含む、前記衝突を表す該特定の次数の非平衡衝突後分布関数を評価するステップと、
を含む1又は2以上の機械可読ハードウェア記憶装置。
【請求項17】
前記格子速度セットは、粒子速度の或る次数まで流体力学運動をサポートする、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項18】
前記格子速度セットに関する前記サポートされる次数は、前記非平衡衝突後分布関数の前記特定の次数より小さくかつこれとは異なり、
前記非平衡衝突後分布関数の前記特定の次数は、前記粒子速度の前記次数によって決定される、請求項17に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項19】
前記容積内の前記特定位置における前記粒子平均速度は、前記特定位置における特定タイプ粒子の粒子平均速度を含む、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項20】
前記格子速度セットは、格子ボルツマン法に関連する状態ベクトルのセットである、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項21】
前記非平衡衝突後分布関数は、(i)予め定義された物理量に対する非平衡モーメントを維持し、(ii)未定義の物理量に対する非平衡モーメントを前記特定の次数まで除去する、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項22】
前記特定の次数は、前記流体速度の格子音速に対する比率に関連した指数値であり、前記格子速度セットは、前記指数値をサポートする、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項23】
前記格子速度セットは、格子に限定される空間における運動量状態のセットを含む、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項24】
前記非平衡衝突後分布関数は、ガリレイ不変なフィルタ処理演算子である、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項25】
前記非平衡衝突後分布関数に基づいて、前記流体容積中での前記粒子の衝突過程をモデル化するステップをさらに含む、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項26】
前記非平衡衝突後分布関数は、流体力学モーメントに対して1次のサポートを提供する格子速度セットに対するマッハ数に関する1次のガリレイ不変衝突演算子
であり、前記衝突演算子は、
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
iは衝突前の前記粒子の速度ベクトルであり、
u(x,t)は時間tでの特定位置xにおける前記粒子間の平均速度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
iは一定重み係数であり、
は非平衡運動量流束である、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項27】
前記非平衡衝突後分布関数は、流体力学モーメントに関して無限次数のサポートを提供する格子速度セットに対する衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
は粒子相対速度であり、
ρは流体密度であり、
は平衡分布関数であり、
は非平衡運動量流束である、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項28】
前記非平衡衝突後分布関数は、流体力学モーメントに関して2次のサポートを提供する、格子速度セットに対するマッハ数に関して2次のガリレイ不変衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
iは衝突前の前記粒子の速度ベクトルであり、
u(x,t)は時間tでの特定位置xにおける前記粒子間の平均速度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
iは一定重み係数であり、
は非平衡運動量流束である、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項29】
前記予め定義された物理量は、前記特定の容積中の前記流体の質量、前記特定の容積中の前記流体の運動量、及び前記特定の容積中の前記流体のエネルギを含む、請求項21に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項30】
前記非平衡衝突後分布関数は、エネルギ流束に関連する衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
は相対的な粒子速度であり、
は平衡分布関数であり、
は非平衡エネルギ流束である、請求項16に記載の1又は2以上の機械可読ハードウェア記憶装置。
【請求項31】
1又は2以上の処理装置と、
オペレーションを実行する1又は2以上の処理装置で実行可能な命令を記憶する1又は2以上の機械可読ハードウェア記憶装置と、
を備えるシステムであって、前記オペレーションは、
格子速度セットを用いて流体容積中の粒子の移動をシミュレートするステップであって、前記粒子の移動が前記粒子間の衝突を引き起こすものである、ステップと、
前記容積内の特定位置における粒子相対速度を決定するステップであって、前記粒子相対速度が(i)前記容積内の前記特定位置における前記容積のゼロフロー基準粒子絶対速度と(ii)前記容積内の前記特定位置における粒子平均速度との間の差である、ステップ、及び
前記粒子相対速度に基づいて、前記衝突を表す非平衡衝突後分布を決定するステップ、により
前記1又は2以上のコンピュータシステムにより、特定の次数は格子速さのセットによりサポートされ、前記格子速度セットの前記格子速さのセットによりサポートされないより高次の項を除外する衝突演算子を含む、前記衝突を表す該特定の次数の非平衡衝突後分布関数を評価するステップと、
を含む、システム。
【請求項32】
前記格子速度セットは、粒子速度の或る次数まで流体力学運動をサポートする、請求項31に記載のシステム。
【請求項33】
前記格子速度セットに関する前記サポートされる次数は、前記非平衡衝突後分布関数の前記特定の次数より小さくかつこれとは異なり、
前記非平衡衝突後分布関数の前記特定の次数は、前記粒子速度の前記次数によって決定される、請求項32に記載のシステム。
【請求項34】
前記容積内の前記特定位置における前記粒子平均速度は、前記特定位置における特定タイプ粒子の粒子平均速度を含む、請求項31に記載のシステム。
【請求項35】
前記格子速度セットは、格子ボルツマン法に関連する状態ベクトルのセットである、請求項31に記載のシステム。
【請求項36】
前記非平衡衝突後分布関数は、(i)予め定義された物理量に対する非平衡モーメントを維持し、(ii)未定義の物理量に対する非平衡モーメントを前記特定の次数まで除去する、請求項31に記載のシステム。
【請求項37】
前記特定の次数は、前記流体速度の格子音速に対する比率に関連した指数値であり、前記格子速度セットは、前記指数値をサポートする、請求項31に記載のシステム。
【請求項38】
前記格子速度セットは、格子に限定される空間における運動量状態のセットを含む、請求項31に記載のシステム。
【請求項39】
前記非平衡衝突後分布関数は、ガリレイ不変フィルタ処理演算子である、請求項31に記載のシステム。
【請求項40】
前記非平衡衝突後分布関数に基づいて、前記流体容積中での前記粒子の衝突過程をモデル化するステップをさらに含む、請求項31に記載のシステム。
【請求項41】
前記非平衡衝突後分布関数は、流体力学モーメントに対して1次のサポートを提供する格子速度セットに対するマッハ数に関する1次のガリレイ不変の衝突演算子
であり、前記衝突演算子は、
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
iは衝突前の前記粒子の速度ベクトルであり、
u(x,t)は時間tでの特定位置xにおける前記粒子間の平均速度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
iは一定重み係数であり、
は非平衡運動量流束である、請求項31に記載のシステム。
【請求項42】
前記非平衡衝突後分布関数は、流体力学モーメントに関して無限次数のサポートを提供する格子速度セットに対する衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
は粒子相対速度であり、
ρは流体密度であり、
は平衡分布関数であり、
は非平衡運動量流束である、請求項31に記載のシステム。
【請求項43】
前記非平衡衝突後分布関数は、流体力学モーメントに関して2次のサポートを提供する、格子速度セットに対するマッハ数に関して2次のガリレイ不変衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
iは衝突前の前記粒子の速度ベクトルであり、
u(x,t)は時間tでの特定位置xにおける前記粒子間の平均速度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
iは一定重み係数であり、
は非平衡運動量流束である、請求項31に記載のシステム。
【請求項44】
前記予め定義された物理量は、前記特定の容積中の前記流体の質量、前記特定の容積中の前記流体の運動量、及び前記特定の容積中の前記流体のエネルギを含む、請求項36に記載のシステム。
【請求項45】
前記非平衡衝突後分布関数は、エネルギ流束に関連する衝突演算子
であり、前記衝突演算子は
に従って定義され、
xは前記容積内の前記特定位置であり、
tは時間における特定時点であり、
iは前記格子速度セットでの格子速度の指数であり、
0は一定格子温度であり、
Iは2階の単位テンソルであり、
τは衝突緩和時間であり、
は相対的な粒子速度であり、
は平衡分布関数であり、
は非平衡エネルギ流束である、請求項31に記載のシステム。
【発明の詳細な説明】
【技術分野】
【0001】
(優先権の主張)
本願は、2013年7月24日出願の米国特許仮出願第61/858,051号に対して米国特許法第119条(e)の下で優先権を主張するものであり、その内容全体は、これにより引用によって組み込まれる。
【背景技術】
【0002】
衝突過程は多粒子系において2つの基本的な動力学的過程の内の1つであり、もう1つは移流過程である。衝突過程は、個々の粒子が相互作用して集団的な挙動を成すために本質的なものである。衝突過程中、質量、運動量及びエネルギは保存則に従って粒子間で交換される。これらの保存則により、関係する粒子間の全体的な質量及び運動量(及び時にはエネルギー)は衝突の前後で不変であることが保証される。
【先行技術文献】
【特許文献】
【0003】
【特許文献1】米国特許第5594671号明細書
【発明の概要】
【課題を解決するための手段】
【0004】
一般に、本明細書は、格子速度セット(lattice velocity set)において流体容積中の粒子の移動をシミュレートする技術を説明しており、移動が粒子間の衝突を引き起こし、シミュレートされた移動に基づいて容積内の特定位置における粒子の粒子相対速度(relative particle velocity)を決定するが、その粒子相対速度は(i)容積内の特定位置において容積の流れがゼロの下で測定される粒子の絶対速度と(ii)容積内の特定位置における1又は2以上(複数)の粒子の平均速度と、の間の差であり、粒子相対速度に基づいて衝突を表す特定の次数(specified order)の非平衡衝突後分布関数(non-equilibrium post-collide distribution function)を決定する。本態様の他の実施形態は、対応するコンピュータシステム、装置、機械可読ハードウェア記憶装置、及び1又は2以上のコンピュータ記憶装置に記録されたコンピュータプログラムが含まれ、各々が本方法の動作及び特徴を実行するように構成されている。1又は2以上のコンピュータのシステムは、作動時、システムにその動作を実行させるソフトウェア、ファームウェア、ハードウェア、又はそれらの組み合わせをシステムにインストールすることによって、特定の処理又は動作を実行するように構成することができる。1又は2以上のコンピュータプログラムは、データ処理装置による実行時にデータ処理装置に特定の処理又は動作を実行させる命令を備えることによってその動作を実行するように構成することができる。
【0005】
前述の及び別の実施形態は、各々選択的に以下の特徴のうちの1又は2以上を単独で又は組み合わせて含むことができる。特に、1つの実施形態は、以下の特徴の全てを組み合わせて含むことができる。特徴は、1又は2以上のコンピュータシステムによって流体力学運動を粒子速度の或る次数までサポートすることができる格子速度セットを提供するステップを含み、シミュレートすることは、1又は2以上のコンピュータシステムでシミュレートすることを含む。また、特徴は、格子速度セットのサポートされる次数(supported order)が、非平衡衝突後分布関数の特定の次数より小さくかつこれとは異なること、非平衡衝突後分布関数の特定の次数が、粒子速度の次数によって決定されることを含む。
【0006】
また、特徴は、容積内の特定位置における1又は2以上の粒子の平均速度がその特定位置における特定種類の粒子の平均速度を有することを含む。また、特徴は、格子速度セットが格子ボルツマン法と関連する状態ベクトルのセットであることを含む。また、特徴は、非平衡衝突後分布関数が、(i)予め定義された物理量に対する非平衡モーメントを維持し、(ii)未定義の物理量に対する非平衡モーメントを特定の次数まで除去することを含む。また、特徴は、特定の次数が流体速度の格子音速に対する比率に関連した指数値(exponential value)であり、格子速度セットが指数値をサポートすることを含む。また、特徴は、格子速度セットが、格子に限定される空間で運動量状態のセットを有することを含む。また、特徴は、粒子相対速度が、容積内の特定位置において容積の流れがゼロの下で測定される粒子の絶対速度から容積内の特定位置における1又は2以上の粒子の平均速度を差し引いたものであることを含む。また、特徴は、非平衡衝突後分布関数がガリレイ不変フィルタ処理演算子(Galilean invariant filtered operator)であることを含む。また、特徴は、非平衡衝突後分布関数に基づいて流体容積中での粒子の衝突過程をモデル化するステップを含む。また。特徴は、非平衡衝突後分布関数が格子速度セットに対するマッハ数に関して1次のガリレイ不変の衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
xは容積内の特定位置であり、tは時間における特定時点であり、iは格子速度セットでの格子速度の指数であり、T0は一定格子温度であり、u(x,t)は時間tでの特定位置xにおける粒子間の平均速度であり、Iは2階の単位テンソルであり、τは衝突緩和時間であり、wiは一定重み係数であり、Πneqは非平衡運動量流束である。
【0007】
また、特徴は、非平衡衝突後分布関数が、流体力学モーメントに関して無限次をサポートする、格子速度セットに対する衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
xは容積内の特定位置であり、tは時間における特定時点であり、iは格子速度セットでの格子速度の指数であり、T0は一定格子温度であり、Iは2階の単位テンソルであり、τは衝突緩和時間であり、
は相対的な粒子速度であり、ρは流体密度であり、
は平衡分布関数であり、
は非平衡運動量流束である。また、特徴は、非平衡衝突後分布関数が、流体力学モーメントに関して2次のサポートを提供する、格子速度セットに対するマッハ数に関して2次のガリレイ不変の衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
xは容積内の特定位置であり、tは時間における特定時点であり、iは格子速度セットでの格子速度の指数であり、T0は一定格子温度であり、ciは衝突前の粒子の速度ベクトルであり、u(x,t)は時間tでの特定位置xにおける粒子間の平均速度であり、Iは2階の単位テンソル(second rank unity tensor)であり、τは衝突緩和時間であり、wiは一定重み係数であり、
は非平衡運動量流束である。また、特徴は、予め定義された物理量が特定の容積中の流体質量、特定の容積中の流体運動量、及び特定の容積中の流体エネルギを有することを含む。
【0008】
また、特徴は、非平衡衝突後分布関数が、エネルギ流束に関連する衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
xは容積内の特定位置であり、tは時間における特定時点であり、iは格子速度セットでの格子速度の指数であり、T0は一定格子温度であり、Iは2階の単位テンソルであり、τは衝突緩和時間であり、
は相対的な粒子速度であり、
は平衡分布関数であり、
は非平衡エネルギ流束である。
【0009】
上述の技術の実装は、方法又は過程、システム又は装置、又はコンピュータアクセス可能媒体上のコンピュータソフトウエアを含むことができる。
【0010】
システム及び方法及び技術は、多相流れに対するShan−Chen方法と格子ボルツマン方式とのような様々なタイプの数値シミュレーション手法を使用して実施することができる。格子ボルツマン方式に関する追加の情報を本明細書で以下に説明する。しかし、本明細書に説明するシステム及び技術は、格子ボルツマン方式を使用するシミュレーションに限定されず、他の数値シミュレーション手法に適用することができる。
【0011】
システム及び技術は、格子ボルツマン方式を使用する格子ガスシミュレーションを使用して実施することができる。従来の格子ガスシミュレーションは、各格子部位での限られた数の粒子を仮定し、粒子は、ビットのショートベクトルによって表される。各ビットは、特定の方向に移動する粒子を表している。例えば、ベクトル中の1ビットは、特定の方向に沿って移動する粒子の存在(1に設定時)又は不在(0に設定時)を表すことができる。そのようなベクトルは、6ビットを有することができ、例えば、値110000は、2つの粒子がX軸に沿って反対方向に移動し、かつY軸及びZ軸に沿って移動する粒子がないことを示している。1組の衝突規則が、各部位での粒子間の衝突の挙動を支配する(例えば、110000ベクトルは、001100ベクトルになることができ、X軸に沿って移動する2つの粒子間の衝突は、Y軸に沿って離れる2つの粒子を生成したことを示す)。規則は、ビットに対して置換を実行する(例えば、110000を001100に変形する)ルックアップテーブルに状態ベクトルを供給することによって実施される。粒子は、その後に、隣接部位に移動される(例えば、Y軸に沿って移動する2つの粒子は、Y軸に沿って左右に隣接部位に移動されると考えられる)。
【0012】
強化されたシステムにおいて、各格子部位での状態ベクトルは、粒子エネルギ及び移動方向の変動を与えるためにより多くのビットを含み(例えば、亜音速流れに対して54ビット)、完全な状態ベクトルの部分集合を伴う衝突規則が使用される。更に強化されたシステムにおいて、1つよりも多い粒子が、各格子部位又はボクセル(これらの2つの用語は、本明細書を通して交換可能に使用される)で各運動量状態で存在することが許される。例えば、8ビット実装において、0〜255個の粒子は、特定のボクセルで特定の方向に移動している可能性がある。状態ベクトルは、1組のビットである代わりに、1組の整数であり(例えば、1組の8ビットバイトが、0〜255の範囲の整数を与える)、その各々は、与えられた状態にある粒子の数を表する。
【0013】
更なる強化において、格子ボルツマン方法(LBM)は、従来の計算流体力学(「CFD」)手法で可能であるよりも深いレベルで複合的な幾何学形状における3D不安定圧縮性乱流過程をシミュレートするために流体のメゾスコピック表現を使用する。LBM方法の概要を以下に与える。
【0014】
ボルツマンレベルのメゾスコピック表現(Boltzmann-Level Mesoscopic Representation)
流体システムは、いわゆる「メゾスコピック」レベルに対する運動方程式によって表すことができることが統計物理学で公知である。このレベルに対しては、個々の粒子の詳細な運動は決定する必要はない。これに代えて、流体の特性は、単一粒子位相空間、
を使用して定義される粒子分布関数によって表され、xは、空間座標であり、νは、粒子速度座標である。質量、密度、流量速度、及び温度のような典型的な流体力学的な量は、粒子分布関数の単純なモーメントである。粒子分布関数の力学は、ボルツマン方程式に従う。
方程式(1)
ここで、F(x,t)は、(x,t)で外部又は自己矛盾なく発生された物体力を表している。衝突項Cは、様々な速度及び位置の粒子の相互作用を表している。衝突項Cに対して特定の形態を指定することなく、上述のボルツマン方程式は、希薄ガスの公知の状況(ボルツマンによって本来構成されたような)だけでなく全ての流体システムに適用可能であることを強調することが重要である。
【0015】
一般的に、Cは、二点相関関数の複雑な多次元積分を含む。分布関数fのみで閉じたシステムを形成するために、並びに効率的な計算のために、最も便利かつ物理的に一貫した形態の1つは、公知のBGK演算子である。BGK演算子は、衝突の詳細がたとえ何であろうとも、分布関数は、衝突:
方程式(2)
を通じて
によって与えられる明確な局所平衡に近づくという物理的議論に従って構成され、ここで、パラメータτは、衝突を通じた平衡までの固有緩和時間を表している。粒子(例えば、原子又は分子)を扱うと、緩和時間は、典型的には定数として取られる。「混成」(水力−動力学)表現において、この緩和時間は、歪み率及び乱流運動エネルギなどのような流体力学的変数の関数である。すなわち、乱流は、局所的に決定された固有特性を有する乱流粒子(「渦」)のガスとして表すことができる。
【0016】
ボルツマン−BGK方程式の数値解は、ナビエ−ストークス方程式の解に優るいくつかの計算上の利点を有する。第1に、複雑な非線形項又は高次空間導関数が方程式にないことを直ちに認識することができ、すなわち、移流不安定性に関する問題がほとんどない。このレベルの記述では、方程式は、圧力を扱う必要がないので局所的であり、これは、アルゴリズム並列化に対してかなりの利点を提供する。2次空間導関数を有する拡散演算子がないという事実と共に線形移流演算子の別の望ましい特徴は、流体偏微分方程式(「PDE」)の数学的な条件ではなく、実際に粒子が固体面と真に相互作用する方法を模倣するように非スリップ面又はスリップ面のような物理的境界条件を実現する際の容易さである。直接の利益の1つは、固体面上でインタフェースの移動を処理する問題がないということであり、これは、格子ボルツマンベースのシミュレーションソフトウエアが失敗なく複合乱流空気力学をシミュレートすることを可能にするのを補助する。これに加えて、有限粗度面のような境界からのある一定の物理特性も、力に組み込むことができる。更に、BGK衝突演算子は、純粋に局所的であり、一方、自己矛盾のない物体力の計算は、近−隣接情報のみを通して達成することができる。その結果、ボルツマン−BGK方程式の計算は、並列処理に実質的に適応させることができる。
【0017】
格子ボルツマン方式
連続体ボルツマン方程式を解くことは、それが位置及び速度位相空間における積分微分方程式の数値評価を伴うということで有意な課題を表している。位置だけでなく速度位相空間も離散化することができることが観測された時に大規模な簡素化が行われ、それは、ボルツマン方程式の解のための効率的な数値アルゴリズムをもたらした。流体力学的な量は、高々最近隣接情報に依存する単純和の項で書くことができる。従来的には、格子ボルツマン方程式の方式は、速度の離散集合:
に対する粒子の推移を規定する格子ガスモデルに基づいていたが、この方程式は、連続体ボルツマン方程式の離散化としての第1の原理から系統的に導出することができる。その結果、LBEは、格子ガス手法に関連付けられた公知の問題を被らない。従って、位相空間における連続体分布関数f(x,ν,t)を取り扱う代わりに、離散速度指数をラベル付けする下付き文字を有する離散分布の有限集合
を追跡することが必要なだけである。巨視的記述の代わりにこの動態方程式を取り扱う重要な利点は、システムの位相空間の増大が、問題の局所性によってオフセットされるということである。
【0018】
対称性の考慮に起因して、速度値の集合は、それらが構成空間内で張られた時にある一定の格子構造を形成するように選択される。そのような離散システムの力学は、形態:
を有するLBEに従い、ここで、衝突演算子は、上述のように通常はBGK形態を取る。平衡分布形態の適正な選択により、格子ボルツマン方程式は正しい流体力学及び熱−流体力学を生じさせることを理論的に示すことができる。すなわち、
から導出された流体力学モーメントは、巨視的限界においてナビエ−ストークス方程式に従う。これらのモーメントは、
方程式(3)
して定義され、ここで、ρ、u、及びTは、それぞれ、流体密度、速度、及び温度であり、Dは、離散化された速度空間の次元(物理空間次元に全く等しくない)である。
【0019】
他の特徴及び利点は、図面及び特許請求の範囲を含めて、以下の説明から明らかとなる。
【図面の簡単な説明】
【0020】
図1】LBMモデルの速度成分を例示する。
図2】LBMモデルの速度成分を例示する。
図3】物理過程シミュレーションシステムが従う手順の流れ図である。
図4】マイクロブロックの斜視図である
図5A図3のシステムで使用する格子構造を例示する。
図5B図3のシステムで使用する格子構造を例示する。
図6】可変分解能技術を例示する。
図7】可変分解能技術を例示する。
図8】表面のファセットにより影響を受ける領域を例示する。
図9】ボクセルから表面への粒子の移動を例示する。
図10】表面から表面への粒子の移動を例示する。
図11】表面動力学を実行するための手順の流れ図である。
図12】特定次数の非平衡衝突後分布関数を決定するするプロセスの流れ図である。
図13】特定次数の非平衡衝突後分布関数を決定するするシステムの構成要素のブロック図である。
【発明を実施するための形態】
【0021】
A.特定の非平衡移動を保持し他の非平衡移動を除去する衝突演算子
格子ボルツマン・シミュレーションなどのシミュレーションシステムでは、シミュレーション空間は多数の離散点に分割され、それらは直線で連結されて離散的な数の点及び方向を与える。シミュレーションはまた、時間ステップの離散集合に制約される。このようなシステムでは、シミュレーションが現実世界の流れを近似するために、多数の異なる量が保存されなければならない。例えば、システムは質量、運動量及びエネルギを保存する。従って、シミュレーションは、適切な質量流束、運動量流束、及びエネルギ流束を有するように構成する必要がある。これらの保存量はその流束と共に、真実の物理世界と関係するシミュレーションシステムにおいて本質的なモーメントである。しかしながら、これらの量を保存する場合に、シミュレーションは離散速度空間(例えば、粒子が所定の時間ステップで進むことのできる方向及び距離の離散集合)に起因して非意図的に付加的なモーメント量を作り出すことがある。非意図的に生成された量(本明細書では、意図しない又は望ましくない不変量、保存される又は非平衡モーメントとして言及する)は、シミュレーション結果に悪影響を与えることがある。例えば、このような望ましくない量は、誤った流体の動的挙動と計算結果の数値的不安定性とをもたらす場合がある。
【0022】
非意図的に生成される不変量の影響を低減するために、保存される物理量に関してのみ非平衡モーメントを保持するが、残りの非平衡モーメントを所望の次数まで除去する、衝突演算子を本明細書に説明する。本明細書では、衝突演算子の次数は格子マッハ数(流体速度と格子音速の比)の指数を用いて定義される。望ましくない非平衡モーメントをフィルタ処理することを確実にするためのトレードオフは、計算時間と処理手続きの増大である。本明細書に説明する衝突演算子では、理論的な形式は運動量及びエネルギの非平衡流束の両方に対する任意の次数まで体系的に構築される。従って、これらの衝突演算子は、質量、運動量及びエネルギの保存を満たし、正確な質量流束、運動量流束及びエネルギ流束を選択した次数まで保証するが、同時に、望ましくない非平衡モーメントのすべてを選択した次数まで除去する。本シミュレーションシステムでは、高速流又は低粘性を取り扱う際に、フィルタ処理のスキームに関してより高次オーダーを選択することが、より低次オーダーの場合に比べて非物理的効果又は影響をより少なくする。
【0023】
B.シミュレーション空間のモデル化
LBMベースの物理過程シミュレーションシステムでは、離散速度ciの集合で評価される分布関数値fiによって、流体の流れを表現することができる。分布関数の力学は、方程式4により支配され、ここにfi(0)は平衡分布関数として知られ、
方程式(4)
として定義され、ここで、
である。
方程式(5)
この方程式は、分布関数fiの時間発展を説明する、よく知られた格子ボルツマン方程式である。左辺はいわゆる「流動過程」による分布の変化を表す。流動過程は、流体ポケットが或るグリッド位置から出発して、それから速度ベクトルの1つに沿って次のグリッド位置へ移動する過程である。その地点で、「衝突演算子」つまり近接した流体ポケットの開始流体ポケットへの影響を計算する。流体は別のグリッド位置に移動することだけが可能なので、全速度の成分全体は共通した速度の倍数であるように速度ベクトルの適切な選択が必要である。
【0024】
第1の方程式の右辺は、流体ポケット間の衝突による分布関数の変化を表す前述の「衝突演算子」である。本明細書で使用する衝突演算子の特別な形式は、Bhatnagar、Gross及びKrook(BGK)によるものである。衝突演算子は、分布関数を、「平衡」形である第2の式により与えられる規定値へ向かうように強いる。
【0025】
このシミュレーションから、密度ρ及び流速uなどの従来の流体変数は、式(3)の単純和として得ることができる。ここで、ci及びwiの集団的な値はLBMモデルを規定する。LBMモデルは、大規模に実現可能なコンピュータ・プラットフォームで効率的に実行することができ、時間的に不安定な流れ及び複雑な境界条件に対して高い堅牢性で走らせることができる。
【0026】
ボルツマン方程式から流体システムに対する巨視的な運動方程式を得るための標準的な技術は、Chapman−Enskog法であり、そこでは完全なボルツマン方程式の逐次近似がなされる。
【0027】
流体システムでは、密度の小さな擾乱が音速で伝播する。気体系では、音速は一般に温度で決まる。流れにおける圧縮性の影響の重要性は、特徴的な速度と音速との比によって評価され、それはマッハ数として知られる。
【0028】
図1を参照すると、第1モデル(2D−1)100は21個の速度を含む2次元モデルである。これら21個の速度の内、1個(105)は動いていない粒子を表し、3組の4速度は、正規化速度(r)(110−113)、正規化速度の2倍(2r)(120−123)又は正規化速度の3倍(3r)(130−133)で、格子のx軸又はy軸に沿って正方向又は負方向に動いている粒子を表し、2組の4速度は、正規化速度(r)(140−143)又は正規化速度の2倍(2r)(150−153)で、2組の4速度は、格子のx軸とy軸の両方に対して動いている粒子を表す。
【0029】
さらに図2に図示するように、第2モデル(3D−1)200は39個の速度を含む3次元モデルであり、各速度は図2の矢印の内の1つにより表されている。これら39個の速度の内、1個は動いていない粒子を表し、3組の6速度は、正規化速度(r)、正規化速度の2倍(2r)又は正規化速度の3倍(3r)で、格子のx軸、y軸又はz軸に沿って正方向又は負方向に動いている粒子を表し、8個は、正規化速度(r)でx、y、zの全3格子軸に対して動いている粒子を表し、12個は、正規化速度の2倍(2r)でx、y、zの格子軸の内の2つに対して動いている粒子を表す。
【0030】
例えば、3D−2は101個の速度を含み、2D−2モデルは37個の速度を含むなど、さらに複雑なモデルも使用することができる。
【0031】
3次元モデル3D−2に関しては、101個の速度の内、1個は動いていない粒子を表し(グループ1)、3組の6速度は、正規化速度(r)、正規化速度の2倍(2r)又は正規化速度の3倍(3r)で、格子のx軸、y軸又はz軸に沿って正方向又は負方向に動いている粒子を表し(グループ2、4及び7)、3組の8速度は、正規化速度(r)、正規化速度の2倍(2r)又は正規化速度の3倍(3r)で、x、y、zの全3格子軸に対して動いている粒子を表し(グループ3、8及び10)、12個は、正規化速度の2倍(2r)でx、y、zの格子軸の内の2つに対して動いている粒子を表し(グループ6)、24個は、正規化速度(r)及び正規化速度の2倍(2r)でx、y、zの格子軸の内の2つに対して動いており、残る1軸に対しては動いていない粒子を表し(グループ5)、24個は、正規化速度(r)でx、y、zの格子軸の内の2つに対して動いており、残る1軸に対しては正規化速度の3倍(3r)で動いている粒子を表す(グループ9)。
【0032】
2次元モデル2D−2に関しては、37個の速度の内、1個は動いていない粒子を表し(グループ1)、3組の4速度は、正規化速度(r)、正規化速度の2倍(2r)又は正規化速度の3倍(3r)で、格子のx軸又はz軸に沿って正方向又は負方向に動いている粒子を表し(グループ2、4及び7)、2組の4速度は、正規化速度(r)又は正規化速度の2倍(2r)で、格子のx軸とy軸の両方に対して動いている粒子を表し、8速度は、正規化速度(r)で格子のx軸とy軸の内の1つに対して、他方の軸に対しては正規化速度の2倍(2r)で動いている粒子を表し、8速度は、正規化速度(r)で格子のx軸とy軸の内の1つに対して、他方の軸に対しては正規化速度の3倍(3r)で動いている粒子を表す。
【0033】
上記のLBMモデルは、2次元及び3次元の流れの数値シミュレーション用の効率的で堅牢な、離散速度の動力学的モデルの特定のクラスを提供する。この種のモデルには、離散速度とこれらの速度に関係する重みとの特定の組が含まれる。速度は速度空間においてデカルト座標のグリッド点と一致し、これにより、離散速度モデル、特に格子ボルツマンモデルとして知られる種類の正確で効率的な実行が容易になる。このようなモデルを用いて、流れを高い忠実度でシミュレートすることができる。
【0034】
図3を参照すると、物理過程シミュレーションシステムは、流体の流れなどの物理過程をシミュレートするために手順300に従って作動する。シミュレーションに先立って、シミュレーション空間をボクセルの集合としてモデル化する(ステップ302)。通常、シミュレーション空間はコンピュータ支援設計(CAD)を用いて生成される。例えば、CADプログラムは、風洞に置かれた微小デバイスを描くために使用することができる。その後に、適切な分解能を備える格子構造を加えるために、及びシミュレーション空間内の物体及び表面に対処するために、CADプログラムにより作成されたデータを処理する。
【0035】
格子の分解能は、シミュレートするシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘度(ν)、流れ中の物体の特徴的な長さ(L)、及び流れの特徴的な速度に関係し、
方程式(6)
である。
【0036】
物体の特徴的な長さは、物体の大規模な特徴を表す。例えば、微小デバイスの周りの流れがシミュレートされる場合には、微小デバイスの高さは、特徴的な長さであると見なすことができる。物体の小領域(例えば、自動車のサイドミラー)周りの流れが関心領域である場合には、シミュレーションの分解能を増大させることができ、又は関心領域の周りには分解能が増大した領域を使用することができる。ボクセルの寸法は、格子の分解能が増加するにつれて減少する。
【0037】
状態空間は、fi(x,t)として表され、fiは、時間tにおいて3次元ベクトルxにより示される格子部位で状態iにある単位容積当たりの要素又は粒子の数(つまり、状態iの粒子密度)を表す。既知の時間増分に対しては、粒子数を単にfi(x)と呼ぶ。1つの格子部位の全ての状態の組み合わせをf(x)として示す。
【0038】
状態の数は、各エネルギ準位内の可能な速度ベクトルの数により決定される。速度ベクトルは、3次元x、y、及びzを備える空間において整数値の線速度から構成される。状態の数は、複数の種シミュレーションに対して増大される。
【0039】
各状態iは、特定エネルギ準位での異なる速度ベクトルを表す(つまり、エネルギ準位0、1又は2)。各状態の速度ciは、以下のように3つの次元の各々におけるその「速度」で示す。
方程式(7)
【0040】
エネルギ準位0状態は、いずれの次元でも動いていない、つまりcstopped=(0、0、0)である静止粒子を表す。エネルギ準位1状態は、3つの次元の内の1つで±1の速度、他の2つの次元では速度0を有する粒子を表す。エネルギ準位2状態は、全ての3つの次元で±1の速度、又は3つの次元のうちの1つで±2の速度及び他の2つの次元で速度0を有する粒子を表す。
【0041】
3つのエネルギ準位の可能な順列の全てを生成することにより、合計39個の可能な状態(1個のエネルギ0状態、6個のエネルギ1状態、8個のエネルギ3状態、6個のエネルギ4状態、12個のエネルギ8状態、及び6個のエネルギ9状態)が得られる。
【0042】
各ボクセル(つまり、各格子部位)は、状態ベクトルf(x)により表される。状態ベクトルはボクセルのステータスを完全に規定し、39個のエントリを含む。39個のエントリは、1個のエネルギ0状態、6個のエネルギ1状態、8個のエネルギ3状態、6個のエネルギ4状態、12個のエネルギ8状態、及び6個のエネルギ9状態に対応する。この速度セットを使用することによって、システムは、実現された平衡状態ベクトルに関してマックスウェル−ボルツマン統計を生成することができる。
【0043】
処理効率のために、ボクセルをマイクロブロックと呼ばれる2x2x2ボリュームにグループ化する。マイクロブロックは、ボクセルの並列処理を可能にし、かつデータ構造に関連づけされたオーバーヘッドを最小化するように構成される。マイクロブロック内のボクセルの略記法は、Ni(n)として規定され、nは、マイクロブロック内の格子部位の相対位置を表し、n∈{0,1,2,...,7}である。マイクロブロックを図4に示す。
【0044】
図5A及び5Bを参照すると、表面Sはシミュレーション空間においてファセットFαの集合として表され、
方程式(8)
ここで、αは、特定のファセットを列挙する指標である。ファセットは、ボクセルの境界に限定されず、通常は、ファセットが比較的小さな数のボクセルに影響を与えるように、ファセットを隣接するボクセルの大きさの程度に又はそれより僅かに小さくサイズ設定される。表面動力学を実行するために、ファセットに特性を割り当てる。特に、各ファセットFαは、単位法線(nα)、表面積(Aα)、中心位置(xα)、及びファセットの表面動力学特性を説明するファセット分布関数(fi(α))を有する。
【0045】
図6を参照すると、処理効率を向上させるために、異なるレベルの分解能をシミュレーション空間の異なる領域で使用することができる。通常は、物体655の周りの領域650が最大関心領域なので、最高分解能でシミュレートされる。粘性の影響は物体からの距離と共に減少するので、低レベルの分解能(つまり、拡張されたボクセルボリューム)を使用して、物体655から増大する距離で離間する領域660、665をシミュレートする。同様に、図7に図示するように、より低レベルの分解能を使用して、物体775のより有意でない特徴部の周りの領域770をシミュレートすることができるが、最高レベルの分解能を使用して、最も有意な特徴部(例えば、先端面及び後端面)の周りの領域780をシミュレートする。周辺領域785は、最低レベル分解能と最大のボクセルとを用いてシミュレートする。
【0046】
C.ファセットに影響を受けるボクセルの特定
図3を再び参照すると、シミュレーション空間をモデル化した状態で(ステップ302)、1又は2以上のファセットにより影響を受けるボクセルを特定する(ステップ304)。ボクセルは、多くの点でファセットの影響を受ける場合がある。最初に、1又は2以上のファセットと交差するボクセルは、交差しないボクセルと比べて低減された容積を有するという点で影響を受ける。これは、ファセット及びファセットにより表される表面下にある物質がボクセルの一部を占めるので発生する。分画因子Pf(x)は、ファセットの影響を受けないボクセルの部分(つまり、流れがシミュレートされる流体又は他の物質が占める可能性がある部分)を示す。交差しないボクセルに対しては、Pf(x)は1に等しい。
【0047】
粒子をファセットに移動するか又は粒子をファセットから受け取ることによって1又は2以上のファセットと相互作用するボクセルも、ファセットにより影響を受けるボクセルとして特定される。ファセットと交差する全てのボクセルは、ファセットから粒子を受け取る少なくとも1つの状態と、粒子をファセットへ移動する少なくとも1つの状態とを含むことになる。ほとんどの場合に、付加的なボクセルは、このような状態を含むことになる。
【0048】
図8を参照すると、0以外の速度ベクトルciを有する各状態iに対して、ファセットFαは、平行六面体Gにより画定される領域から粒子を受け取るか又はこの領域へ粒子を移動するが、この平行六面体は速度ベクトルciとファセットの単位法線とのベクトルドット積(|cii|)の大きさにより規定される高さとファセットの表面積Aαにより規定される底面とを有するので、平行六面体Gの容積Vは次のようになる。
方程式(9)
【0049】
ファセットFαは、状態の速度ベクトルがファセットに向けられた時に(|cii|<0)、容積Vから粒子を受け取り、状態の速度ベクトルがファセットから離れるように向けられた時に(|cii|>0))、その領域へ粒子を移動する。以下に説明するように、この式は、別のファセットが平行六面体Giαの一部、すなわち、内部のコーナのような非凸面特徴部の近くに発生する可能性がある状態を占有した時には修正する必要がある。
【0050】
ファセットFαの平行六面体Gは、複数のボクセルの内の一部又は全てと重なる場合がある。ボクセル又はその部分の数は、ボクセルのサイズに対するファセットのサイズ、状態のエネルギ、及び格子構造に対するファセットの向きに依存する。影響を受けるボクセルの数は、ファセットのサイズと共に増加する。従って、ファセットのサイズは、上述のように、通常はファセット近くに位置するボクセルのサイズと同程度に又はそれより小さいように選択される。
【0051】
平行六面体Gと重なるボクセルN(x)の一部分は、V(x)として定義する。この項を使用して、ボクセルN(x)とファセットFαとの間を移動する状態i粒子の流束Γ(x)は、ボクセル(Ni(x))における状態i粒子の密度にボクセルと重なる領域の容積(V(x))を乗じたものに等しい。
方程式(10)
【0052】
平行六面体Gが1又は2以上のファセットと交差する場合、以下の条件が正しく、
方程式式(11)
ここで、第1の和は、Gと重なる全ボクセルに相当し、第2の項は、Gと交差する全てのファセットに相当する。平行六面体Gが別のファセットと交差しない場合、この式は以下に帰する。
方程式(12)
【0053】
D.シミュレーションの実行
1又は2以上のファセットにより影響を受けるボクセルが特定された状態で(ステップ304)、シミュレーションを開始するためにタイマを初期化する(ステップ306)。シミュレーションの各時間増分中に、粒子のボクセルからボクセルへの移動は、粒子の表面ファセットとの相互作用に対応する移流段階(ステップ308−316)によりシミュレートされる。次いで、衝突段階(ステップ318)は、各ボクセル内の粒子の相互作用をシミュレートする。次に、タイマを増分する。増分されたタイマがシミュレーションの完了(ステップ322)を示さない場合、移流段階及び衝突段階(ステップ308−320)が繰り返される。増分されたタイマがシミュレーションの完了(ステップ322)を示す場合、シミュレーションの結果を記憶及び/又は表示する(ステップ324)。
【0054】
1.表面に対する境界条件
表面との相互作用を正確にシミュレートするために、各ファセットは、4つの境界条件を満たす必要がある。第1に、ファセットが受け取る粒子の合計質量は、ファセットが移送する粒子の合計質量に等しくなければならない(つまり、ファセットに対する正味の質量流束は0でなければならない)。第2に、ファセットが受け取る粒子の合計エネルギは、ファセットが移送する粒子の合計エネルギに等しくなければならない(つまり、ファセットに対する正味のエネルギ流束は0でなければならない)。これら2つの条件は、各エネルギ準位(つまり、エネルギ準位1及び2)での正味質量流束が0に等しいことを必要とすることにより満たすことができる。
【0055】
他の2つの境界条件は、ファセットと相互作用する粒子の正味運動量に関連する。本明細書では滑り面として言及する、表面摩擦のない表面に対しては、正味接線方向運動量流束は0でなければならず、正味法線方向運動量流束はファセットでの局所的な圧力に等しくなければならない。従って、受け取った合計運動量と移送した合計運動量の、ファセットの法線nαに対して垂直な成分(つまり、接線方向成分)は等しくなければならないが、受け取った合計運動量と移送した合計運動量の、ファセットの法線nαに対して平行な成分(つまり、法線方向成分)の差はファセットでの局所的な圧力に等しくなければならない。非滑り面に対しては、表面の摩擦は、ファセットが受け取る粒子の合計接線方向運動量に対するファセットが移送する粒子の合計接線方向運動量を摩擦量に関連する係数だけ減少させる。
【0056】
2.ボクセルからファセットへの集まり
粒子と表面間の相互作用をシミュレートする場合の第1のステップとして、粒子は、ボクセルから集めてファセットへ供給する(ステップ308)。前述のように、ボクセルN(x)とファセットFαとの間の状態i粒子の流束は、以下のようになる。
方程式(13)
【0057】
この式から、ファセットFα(ciα<0)に向けられた各状態iに対して、ボクセルからファセットFαに供給される粒子数は、以下のようになる。
方程式(14)
【0058】
(x)が0以外の値を有するボクセルのみを合計する必要がある。前述のように、ファセットのサイズは、V(x)が少数のボクセルに対してのみ0以外の値を有するように選択される。V(x)及びPf(x)が非整数値を取る場合があるので、Γα(x)は、実数として記憶及び処理される。
【0059】
3.ファセットからファセットへ移動
次に、粒子をファセット間で移動させる(ステップ310)。ファセットFαの流入状態(ciα<0)の平行六面体Gが別のファセットFβと交差する場合、ファセットFαが受け取る状態iの粒子の一部は、ファセットFβから到来することになる。特に、ファセットFαは、以前の時間増分中にファセットFβにより生成される状態i粒子の一部を受け取ることになる。この関係を図10に示すが、ファセットFβと交差する平行六面体Gの部分1000は、ファセットFαと交差する平行六面体Gの部分1005に等しい。前述のように、交差する部分は、V(β)として示す。この項を使用して、ファセットFβとファセットFαとの間の状態i粒子の流束は、以下のように説明することができ、
方程式(15)
ここで、Γi(β,t−1)は、前回の時間増分中にファセットFβによって生成された状態i粒子の尺度である。この式から、ファセットFα(ciα<0)に向けられた各状態iに対して、他のファセットによりファセットFαに供給される粒子数は、
方程式(16)
であり、ファセットへの状態i粒子の全流束は、以下のようになる。
方程式(17)
【0060】
ファセット分布関数とも呼ばれるファセットに対する状態ベクトルN(α)は、ボクセル状態ベクトルのMエントリに対応するMエントリを有する。Mは、離散格子速度の数である。ファセット分布関数N(α)の入力状態は、それらの状態への粒子の流束を容積Vで割ったものに等しく設定され、ciα<0に対して以下のようになる。
方程式(18)
【0061】
ファセット分布関数は、ファセットからの出力流束を生成するためのシミュレーションツールであり、必ずしも実際の粒子を表すものではない。正確な出力流束を生成するために、他状態の分布関数に値を割り当てる。外向きの状態は、内向きの状態をポピュレートするために前述の技術を使用してポピュレートされ、ciα≧0に対して
方程式(19)
となり、ここで、Γiother(α)は、ΓiIN(α)を生成するための前述の技術を用いるが、流入状態(ciα<0)以外の状態(ciα≧0)にこの技術を適用して決定する。代替技術では、前回の時間ステップからΓiOUT(α)の値を用いて、Γiother(α)を生成することができるので、
方程式(20)
となる。
【0062】
平行状態(ciα=0)に対しては、V及びV(x)は共に0である。Ni(α)の式において、Viα(x)は分子に現れ(Γiother(α)の式から)、Vは、分母に現れる(Ni(α)の式から)。従って、平行状態に対するNi(α)は、V及びV(x)が0に近づく時のNi(α)の極限値として決定される。
【0063】
0速度を有する状態(つまり、静止状態及び状態(0、0、0、2)及び(0、0、0、−2))の値は、温度と圧力に関する初期条件に基づいてシミュレーションの最初に初期化される。これらの値は、その後、時間と共に調整される。
【0064】
4.ファセットの表面力学の実行
次に、各ファセットが前述の4つの境界条件を満たすように、表面力学を実行する。ファセットに対して表面力学を実行するための手順を図11に示す。最初に、ファセットFαに垂直な合計運動量は、ファセットにおける粒子の合計運動量P(α)を確定することによって、以下のように決定され(ステップ1105)、全てのiに対して
方程式(21)
これから、法線方向運動量Pn(α)は、以下のように決定される。
方程式(22)
【0065】
その後、この法線方向運動量は、Nn-(α)を生成するために押込/引抜技術(ステップ1110)を用いて除去する。この技術に従って、法線方向運動量にのみ影響を及ぼすように状態間で粒子を移動させる。この押込/引抜技術は米国特許第5,594,671号に記載されており、この特許は、引用により組み込まれている。
【0066】
その後、Nn-(α)の粒子を衝突させてボルツマン分布Nn-β(α)を生成する(ステップ1115)。流体力学の実行に関して以下で説明するように、ボルツマン分布は、1組の衝突規則をNn-(α)に適用することによって実現することができる。
【0067】
その後、ファセットFαに関する流出流束分布を流入流束分布及びボルツマン分布に基づいて決定する(ステップ1120)。最初に、流入流束分布Γi(α)とボルツマン分布との差を以下のように決定する。
方程式(23)
【0068】
この差分を用いると、流出流束分布は、nαi>0に対して
式(24)
であり、ここで、i*は、状態iと反対の方向を有する状態である。例えば、状態iが(1、1、0、0)の場合、状態i*は、(−1、−1、0、0)である。表面摩擦及び他の要因を考慮するために、流出流束分布を以下のようにさらに精緻化することができる。
αi>0に対して、
方程式(25)
ここで、Cfは、表面摩擦関数であり、tは、nαに垂直な第1の接線ベクトルであり、tは、nα及びtに垂直な第2の接線ベクトルであり、ΔNj,1及びΔNj,2は、状態i及び表示される接線ベクトルのエネルギ(j)に対応する分布関数である。分布関数は以下の方程式に従って決定する。
方程式(26)
ここで、jは、エネルギ準位1に対して1、エネルギ準位2に対して2に等しい。
【0069】
ΓiOUT(α)の方程式の各項の関数は以下の通りである。第1及び第の2項は、衝突がボルツマン分布を生み出すことに有効な範囲で法線方向運動量流束の境界条件を実行するが、接線方向運動量流束偏差を含む。第4及び第5の項は、不十分な衝突による離散効果又は非ボルツマン構造に起因して生ずる場合があるこの偏差を補正する。最後に、第3の項は、表面上の接線方向運動量の所望の変化を強いるために特定の量の表面分画を追加する。摩擦係数Cfの生成を以下に説明する。ベクトル演算を伴う全ての項は、シミュレーションの開始に先立って計算することのできる幾何形状因子であることに留意されたい。
【0070】
これにより、接線方向速度を以下のように決定する。
方程式(27)
ここで、ρは、ファセット分布の密度である。
方程式(28)
【0071】
前と同様に、流入流束分布とボルツマン分布との差を以下のように決定する。
方程式(29)
【0072】
従って、流出流束分布は、以下のようになる。
方程式(30)
この方程式は、前回の技術によって求めた流出流束分布の最初の2行に対応するが、特異な接線方向流束に対する補正を必要としない。
【0073】
あらゆる手法を用いても、結果として生じる流束分布は、運動量流束の条件の全てを満たし、つまり、
方程式(31)
であり、ここで、pαは、ファセットFαにおける平衡圧力であり、ファセットに粒子を与えるボクセルの平均密度及び温度値に基づいており、uαは、ファセットの平均速度である。
【0074】
質量及びエネルギの境界条件が満たされることを保証するために、入力エネルギと出力エネルギとの差を、以下のように各エネルギ準位jに対して測定する。
方程式(32)
ここで、指標jは、状態iのエネルギを示す。その後、このエネルギ差を使用して差の項を生成する。
jiα>0に対して、
方程式(33)
この差の項は、流束が以下になるように流出流束を修正するために使用する。
jiα<0に対して、
方程式(34)
この演算は、質量及びエネルギの流束を補正するが、接線方向の運動量流束を不変のままにする。この調整は、流れがファセットの近傍でほぼ均一で平衡状態に近い場合は小さい。この調整の後に、結果として得られた法線方向運動量流束は、近傍の不均一性又は非平衡特性に起因する補正を加えて近傍の平均特性に基づいて平衡圧力である値に僅かに変更される。
【0075】
5.ボクセルからボクセルへの移動
再び図3を参照すると、粒子を3次元直線格子に沿ってボクセル間で移動させる(ステップ314)。このボクセルからボクセルへの移動は、ファセットとは相互作用しないボクセル(つまり、表面近くに位置しないボクセル)で実行される唯一の移動演算である。典型的なシミュレーションでは、表面と相互作用するほどには近くに位置していないボクセルは、ボクセルの大多数を構成する。
【0076】
異なる状態の各々は、3次元x、y及びzの各々において整数速度で格子に沿って移動する粒子を表す。整数速度は、0、±1及び±2を含む。速度の符号は、粒子が対応する軸に沿って移動している方向を示す。
【0077】
表面と相互作用しないボクセルに関して、移動演算は、計算的にはきわめて簡単である。各時間増分中に、状態の全粒子数を現在のボクセルから目的地ボクセルへ移動させる。同時に、目的地ボクセルの粒子をボクセルから固有の目的地ボクセルへ移動させる。例えば、+1x方向(1、0、0)に移動しているエネルギ準位1の粒子は、現在のボクセルからx方向に+1、他の方向に対して0のボクセルに移動させる。粒子は、結局、移動(1、0、0)の前に有していたのと同じ状態で目的地ボクセルに行く。ボクセル内の相互作用は、他の粒子及び表面との局所的な相互作用に基づいて、その状態に対して粒子総数を変化させることになる。そうでない場合は、粒子は同じ速度及び同じ方向で格子に沿って移動し続けることになる。
【0078】
移動演算は、1又は2以上の表面と相互作用するボクセルに対しては少し複雑となる。これにより、1又は2以上の分画粒子は、ファセットに移送される。このような分画粒子のファセットへの移送により分画粒子はボクセル内に残る。これらの分画粒子は、ファセットで占有されるボクセルへ移送される。例えば、図9を参照すると、ボクセル905の状態i粒子の部分900がファセット910に移動した場合(ステップ308)、残りの部分915は、ボクセル920へ移動させるが、そこにはファセット910が位置しており、状態iの粒子はそこからファセット910へ向けられる。従って、状態粒子数が25に等しく、V(x)が0.25に等しい場合(つまり、ボクセルの4分の1が平行六面体と交差する)、6.25個の粒子は、ファセットFαに移動させることになり、18.75個の粒子は、ファセットFαにより占有されたボクセルへ移動させることになる。複数のファセットが単一のボクセルと交差することができるので、1又は2以上のファセットにより占有されたボクセルN(f)に移送される状態の粒子数は、以下のようになる。
方程式(35)
ここで、N(x)は、ソースボクセルである。
【0079】
6.ファセットからボクセルへの散乱
次に、各ファセットからの流出粒子をボクセルに散乱させる(ステップ316)。本質的に、このステップは、粒子をボクセルからファセットへ移動させた収集ステップの逆である。ファセットFαからボクセルN(x)へ移動する状態iの粒子数は、以下の通りである。
方程式(36)
ここで、Pf(x)は、部分的なボクセルの容積減少に対処する。これから、各状態iに対して、ファセットからボクセルN(x)に向けられる全粒子数は以下のようになる。
方程式(37)
【0080】
粒子をファセットからボクセルまで散乱させ、周囲のボクセルから流入した粒子と結合して結果を整数化した後に、特定のボクセルの特定の方向がアンダーフローするか(負になる)、又はオーバーフロー(8ビット実装で255を超える)する可能性がある。これは、これらの量が許容範囲の値に適合するように切り捨てられた後に、質量、運動量、及びエネルギの取得又は損失をもたらす。このような事象の発生を防ぐために、範囲を越えた質量、運動量、及びエネルギが、問題のある状態の切り捨てに先立って蓄積される。状態が属するエネルギに関して、得られる(アンダーフローによる)又は失われる(オーバーフローによる)値に等しい量の質量は、同じエネルギを有しかつそれ自体オーバーフロー又はアンダーフローを受けないランダムに(又は連続的に)選択された状態に追加される。この質量及びエネルギの追加から生じる追加の運動量は蓄積され、切り捨てられた運動量に追加される。同じエネルギ状態にだけ質量を追加することにより、質量カウンタが0に到達した場合に質量及びエネルギの両方が補正される。最後に、運動量アキュムレータが0に戻るまで、運動量は、押込/引抜技術を利用して補正される。
【0081】
7.流体力学の実行
最後に、流体力学を実行する(ステップ318)。このステップは、微視的動力学又はボクセル内演算と呼ぶことができる。同様に、移流手順は、ボクセル内演算と呼ぶことができる。また、以下に説明する微視的動力学演算を使用して、ファセットで粒子を衝突させてボルツマン分布を生成することができる。
【0082】
流体力学は、BGK衝突モデルとして知られる特定の衝突演算子によって格子ボルツマン方程式において保証される。この衝突モデルは、実際の流体システムの分布の動力学を模倣している。衝突過程は、方程式1の右辺及び式2によりうまく説明することができる。移流ステップの後で、流体システムの保存量、具体的には、密度、運動量、及びエネルギは、方程式3を用いて分布関数から取得する。これらの量から、方程式(2)のfeqで説明される平衡分布関数は、方程式(4)により完全に特定される。速度ベクトルセットci及び重みの選択は表1に列挙され、方程式2と一緒に、巨視的挙動が正確な流体力学の方程式に従うことを保証する。
【0083】
E.衝突過程
関連する流体物理を再現するために、格子ボルツマンシステムでの衝突過程は、実際の流体システムと同様の基本的な役割を果たし、同様の基本的な保存要件に従う。便宜上、以下の方程式の番号は、方程式1.1)から始まる。
は、「衝突前」粒子分布関数(つまり、時間tにおいて位置xにある単位容積中の粒子数)とし、この分布関数は衝突後に
に変わる(つまり、「衝突後」分布関数)。質量の保存は以下の方程式を満たす。
ここで、ρ(x、t)は、位置x及び時間tでの全速度ベクトル値の粒子数密度に等しい流体密度である。(1.1)(及び以下の方程式)の合計は、格子ボルツマンモデルにおける可能な粒子速度ベクトル値の全てに対するものである。運動量の保存は、
によって与えられ、ここでu(x,t)は流体速度であり、単に、位置x及び時間tにおける粒子間の平均速度である。
【0084】
粒子の運動エネルギも保存される特定の流体システム(多粒子で構成される)に関しては、(1.1)及び(1.2)に加えて、次の関係も定義される。
ここで、
であり、定数Dは粒子運動の次元である。T(x,t)は、x及びtにおける流体システムの温度である。
【0085】
保存則から明らかであるが、ρ(x、t)及びu(x,t)(及びエネルギ保存システムではT(x、t))は、衝突過程中には不変である。所定のρ(x,t)及びu(x,t)(及びT(x,t))に対して、平衡分布関数とも呼ばれる、特別な形式の粒子分布関数
が存在する。平衡分布関数は、式(1.1)及び(1.2)(及び(1.3))に定義されるものと同じ質量及び運動量(及びエネルギー)を有し、ρ(x,t)及びu(x,t)(及びT(x,t))の関数として完全に決定される。
【0086】
分布関数の合計(粒子速度ベクトル値に対する)に関する量は、一般に、分布関数のモーメントと呼ばれる。質量、運動量、及びエネルギの巨視的保存量に対応する、(1.1)−(1.3)の3つの基本的なモーメントの他に、それらの流束に対応するモーメントは、同様に重要である。これらは、衝突後の分布関数の合計によって定めることができる。
ここで、Π(x,t)及びQ(x,t)はそれぞれ、位置x及び時間tにおける運動量及びエネルギの流束テンソルである。(1.2)の定義を想起されたい。運動量流束及びエネルギ流束だけが上記の巨視的保存量((1.1)−(1.3)で定義される)に無関係である。
【0087】
平衡状態での分布に関して(つまり、(1.5)及び(1.6)で代わりに
を用いる)、平衡状態の運動量流束及びエネルギ流束は、それぞれ基本的な(連続)気体分子運動理論から公知の形態
を有し、ここで、p(x,t)は、圧力であり、理想気体の法則に基づいてp(x,t)=ρ(x,t)T(x,t)であり、(1.7)の「I」は2階の単位テンソルを示す。Tr[ ]はトレース演算である。格子ボルツマン法での中心テーマは、(1.7)及び(1.8)の形態を取り戻すことである。
【0088】
(1.5)及び(1.6)のモーメントは、平衡部及び非平衡部の項で以下のように表すこともできる。
【0089】
明らかに、
である。非平衡運動量流束は、流体システムでの輸送挙動を決定する際に重要な役割を果たす。例えば、非平衡運動量流束(1.9)は、ニュートン流体システムにおいて粘性を直接決定するが、(1.10)の非平衡エネルギ流束は、温度伝導率を決定する。粒子分布関数から無数のモーメントを構成することができるが、(1.1)−(1,10)によるモーメントだけが物理的な流体システムにおいて巨視的に関連する。
【0090】
本システムは、平衡状態からの偏差の程度を
として示し、従って、物理的に安定な衝突過程は、常にこの偏差を低減する方向に働く。すなわち、平衡状態からの衝突後の偏差は、衝突前の偏差よりも小さい。
【0091】
実際の多粒子システムでの衝突は非常に複雑な場合がある。1つの最も簡単な衝突モデル(演算子)は、保存則要件と(1.11)の平衡状態への収束要件との両方を満たす。「BGK」衝突過程は、数学的には以下のように表され、
衝突演算子は、以下の通りである。
【0092】
3つの分布関数の全てが同じρ(x,t)及びu(x,t)(及びエネルギ保存システムに対してはT(x,t))を与えるので、質量及び運動量(及びエネルギー)の保存(つまり、(1.1)−(1.3))は、自動的に満たされる。さらに、(1.12)から、衝突後分布関数に関する平衡状態からの偏差は、衝突前分布関数に関する平衡状態からの偏差に係数(1−1/τ)で比例する。従って、収束条件(1.11)は、パラメータ値τ(衝突緩和時間と呼ぶ)が1/2よりも大きい限り満たされ、衝突後偏差は、τ=1の時にゼロになる。
【0093】
BGK関係式(1.12)及び(1.13)から、方程式(1.9)及び(1.10)は、以下のように書き換えることができ、
ここで、Πneq(x,t)及びQneq(x,t)はそれぞれ、以下に定義される衝突前の非平衡状態の運動量流束及びエネルギ流束である。
【0094】
τの値が選択された状態で、BGK衝突演算子(1.12)を用いて格子ボルツマン流体における動粘性率値は、以下のように決定される。
【0095】
ここで、定数T0は標準格子温度であり、格子単位系での時間増分が1となるように、
いわゆる格子単位規則を用いる。BGK衝突演算子は、格子ボルツマンモデルにおいて最も一般的に用いられる演算子である。BGK衝突演算子は、過去20数年で種々レベルの成功を収めた。一方で、BGK演算子にはいくつかの本質的な制限がある。単一プラントル数(=粘性率と熱拡散率の比)の他に、(1.9)−(1.10)の基本的な問題に加えて、1つの問題は、τが1でない場合は、全ての衝突後非平衡モーメントが生成される点にある。実際にBGKに関して、衝突前後の非平衡モーメントは、以下の関係を示し、
ここで、Mn′(x,t)及びMn(x,t)はそれぞれ、衝突前及び衝突後の任意のn次モーメントを表す。
【0096】
F.フィルタ衝突演算子(Filtered Collision Operator)
任意の格子ボルツマンモデルの最も一般的な特徴は、有限かつ一定の粒子速度ベクトル値セットを有することである。所定の格子ボルツマンモデルに対して、定数ベクトル値セットが特定される。その結果、離散速度セットから構成された粒子分布関数のモーメントの有限セットだけが、実際の流体におけるそれらの対照物を取り戻すことができる。実際の流体のモーメントを任意の所定の次数まで取り戻すための一般的な構成は、厳密に規定される。格子ボルツマンベクトル値の異なるセットは、物理的モーメントを異なる次数までサポートすることができる。例えば、いわゆるD3Q19及びD3Q15は、平衡運動量流束と低マッハ数の制限下での線形偏差の(ニュートンの)非平衡運動量流束までのモーメントをサポートするだけである。一方では、D3Q39等のいわゆる高次格子ボルツマンモデルは、有限マッハ数及びクヌーセン数での線形偏差を超える平衡エネルギ流束及び非平衡運動量流束までモーメントをサポートすることができる。
【0097】
物理的に関連するモーメントは、平衡及び非平衡運動量、並びにエネルギのモーメント及びそれらの流束だけなので、正確な物理的流体の挙動を実現するためには、非平衡運動量流束(場合によってはエネルギー流束)以外の全ての衝突後非平衡モーメントがゼロになるように、衝突演算子を設計することが望ましい。具体的には、非平衡運動量流束及びエネルギ流束は、以下の方程式で表されるが、
他の衝突後非平衡モーメントはゼロである。
【0098】
(1.19)のBGK衝突は、全ての非平衡モーメントを生成することを想起されたい。また、(2.1)−(2.3)を実現する衝突演算子は、非本質的な非平衡モーメントをフィルタ除去するので、フィルタ衝突演算子と呼ばれる。フィルタ衝突演算子は、他の非平衡モーメントを全て除去しながら、非平衡運動量流束を(2.1)のように維持する。このフィルタ衝突演算子の明示的表現は以下の通りである。
ここで、定数wiは、特定の格子速度セットが選択された状態で決まる重み係数値である。格子速度セットは、格子に制限される空間での微視的速度又は運動量の離散セットである。重み係数値のセットは、異なる格子ボルツマン粒子速度セットに対して異なり、その目的は、期待される次数までモーメントの等方性を実現することである。新しい(フィルタ処理された)衝突演算子(2.4)は、BGKと同じ衝突後非平衡運動量流束(1.9)をもたらし、自動的にBGKと同じ(1.18)の粘性作用を引き起こす。
【0099】
特定の格子システムに関して、1ではないτeを有する方程式(2.2)に示すように、非ゼロの非平衡エネルギ流束を得ることが望ましい。このような目的を実現するために示された特定の衝突演算子の形式は以下の通りである。
【0100】
値τeは、τと等しい必要はないので、プラントル数は、BGKの場合とは対照的に単一の値に限定されないことに留意されたい。運動量流束とエネルギ流束との間のモーメント直交性により、より一般化されたフィルタ衝突演算子形式は(2.4)及び(2.5)の2つの形式の直和により構成されるので、(2.1)及び(2.2)の非平衡流束は、両方とも自動的に満たされるが、残りの非平衡モーメントはゼロになる。
【0101】
所定の格子速度セットではサポートされない不必要な非平衡モーメントをフィルタ除去することによって、一般化された衝突演算子((2.4)、(2.5)又はこれらを組み合わせたもの)は、BGK衝突演算子よって著しく改善された流体流れ等方性を示し、BGK衝突演算子と同様に所望の運動量流束及びエネルギ流束を維持する。一般化された衝突演算子の形式(2.4)及び(2.5)(及びこれらを組み合わせたもの)の態様は、ナビエ−ストークス流体の粘性率及び熱拡散率に適用可能であるだけでなく、有限なクヌーセン数を含むより広範な流体レジームにおいて正確な流体力学を保証する。
【0102】
G.N次のガリレイ不変フィルタ演算子を導出する手順
本開示に適合するシステムは、どれだけの量の速度がシミュレーションによってサポート可能かに基づいて、格子マッハ数の特定の指数を有するフィルタ衝突演算子を生成する。このフィルタ衝突演算子では、絶対値(例えば、速度又はエネルギの)ではなく相対値(例えば速度又はエネルギーの)を使用する。本明細書では、不必要なモーメントを含まない衝突演算子に対する正しい理論形式を適切に構成する手順を説明する。非平衡モーメントの総数は、シミュレーションによってサポート可能な速度の総数に対応する。従って、フィルタ衝突演算子は、所定の格子ボルツマンモデル粒子速度セットによってサポートされる非平衡モーメントを表すので、非平衡モーメントは物理的世界で実際に起こるものに対応する。相対速度を用いることにより、シミュレーションの速度によってサポートされる所望の非平衡状態の寄与部を備えたフィルタ衝突演算子の構成が可能となる。
【0103】
すなわち、非平衡状態の運動量流束及びエネルギ流束を相対的な粒子速度及びエネルギを用いて決定することにより、シミュレーションにおいて格子ボルツマンモデルの格子速さのセットによりサポートされない、より高次の項をフィルタ除去する(例えば、除外する)フィルタ衝突演算子を構成することができる。1つの実施例では、フィルタ衝突演算子は拡張形式で特定されるので、多数の格子速度に関して、ユーザ又はシステムは、どの次数までの項を保存する必要があり、どの次数を超える項を除去する必要があるか(フィルタ衝突演算子から)を知る方法を有する。そうするために、拡張部(例えば、フィルタ衝突演算子)を相対速度で表す必要がある。
【0104】
1つの実施例では、速度モデルは、粒子速度の有限セットである。従って、粒子移動の本当の物理モーメントだけが、特定の次数までこのモデルで正確に表すことができる。ガリレイ不変衝突演算子を有しかつより速い流速で速度をシミュレートするために、衝突演算子の形式は相対速度に基づいており、ここでは粒子速度は自身の流速に対して測定される。少なくともこの理由から、パラメータciを以下の方程式3.5で使用するのではなく、パラメータci’(相対速度又は相対エネルギーを表す)を方程式3.5で使用する。方程式3.5でのパラメータciの使用は、流速(例えば、u(x,t))の累乗でこの衝突演算子形式を展開することを許さないので、局所的な流速に対する相対速度に関してガリレイ不変性に従う適切な衝突演算の構成をもたらさない。しかしながら、方程式3.5でのパラメータci’の使用は、流速(例えば、u(x,t))の累乗(例えば、N次)でこの衝突演算子形式を展開することを可能にするので、N次のガリレイ不変フィルタ演算子をもたらす。
【0105】
1つの実施例では、BGKモデルの下での衝突演算子は、無限次数の非平衡モーメントを含む。特定の速度モデルは無限次数まで正確にモーメントをサポートできないので、特定の次数までの非平衡モーメントは、非物理的なものとなる。従って、このシステムは、高次の非平衡モーメントを除去して物理的人工物を含まないようにする必要がある。例えば、19速度モデルに関連するモーメントは、一次までである。従って、より高次のモーメントは不適切である。しかしながら、BGK演算子を用いて研究シミュレーションを実行する場合、モデルによりサポートされるものを超える非平衡モーメントが含まれる。従って、本明細書に説明するフィルタ衝突演算子は、物理的世界で起こっているものにより密接に対応する。
【0106】
前述のフィルタ衝突演算子形式は、非常に低い流速の流れ状況に対応するのみである。ガリレイ不変性の基本原理から、物理的流体の統計的特性の全ては、平均流束に対する粒子相対速度だけの関数である必要がある。具体的には、物理的流体システムにおける粒子非平衡分布関数及び関連の非平衡モーメントは、絶対速度ci(格子ボルツマンシステムのゼロフロー基準座標系で測定される)ではなく、相対速度(ci−u(x,t))にのみ依存するべきである。従って、(1.14)及び(1.15)の代わりに、より物理的に意味のある非平衡の運動量流束及びエネルギ流束は、それぞれ以下の通りである。
ここで、粒子相対速度及びエネルギは、
で与えられる。上記方程式(3.3)に示すように、粒子絶対速度は、相対速度で置換される。興味深いことに、質量及び運動量の保存により、非平衡運動量流束(3.1)は(1.16)と同じであることが分かる。一方で、非平衡エネルギ流束(3.2)は、(1.17)に変換できない。
【0107】
本開示に適合するシステムは、非平衡運動量流束に関するガリレイ不変衝突演算子を生成する。気体分子運動論の基本法則に従って、流速の不均一性に起因する、1次(leading order)の非平衡粒子分布関数は、以下の通りである。
ここで、Uは、連続体気体分子運動論での局所的平均流れに対する粒子相対速度である。この概念に端を発して、本明細書に説明したシステムは、ボルツマン法において、ガリレイ不変性に一致する非平衡分布関数に類似した式を特定する。本システムは、完全に対応する非平衡衝突後分布関数として、以下の明示形式を特定する。
ここで、Πneq(x,t)は、(1.14)又は(3.1)で与えられる。上記の方程式3.2に示すように、分布関数は、平衡成分及び非平衡成分を含む。さらに、上記の方程式は、粒子速度セットを無限次数まで必要とする。(3.5)の非平衡分布関数は、任意のマッハ数に対する完全形式である。
【0108】
それぞれ非平衡及び平衡の分布関数に対する(3.5)及び(3.6)の完全形式は、格子速度セットが全次数まで正確な流体力学モーメントをサポートする場合にのみ実現可能である。換言すると、粒子速度値の有限セットを備える任意の所与の速度セットに対して、完全形式は実現できない。しかしながら、以下に詳細に説明するように、これらを所定の有限格子速度セットを用いて対応する次数Nまで実現することが可能である。
【0109】
式(3.5)及び(3.6)は、u(x,t)の累乗で多項式形式に展開することができる。平衡分布関数(3.6)は、エルミート多項式級数を用いて表現できることはとはよく知られている。
ここで、V(x,t)[n]は、ベクトルV(x,t)n次直積に関する簡略表記である。
n次エルミート多項式H(n)(ξi)は、標準(スカラー)n次エルミート関数のn階テンソル一般化である。裏付けはないが、非平衡の寄与分(3.5)は、エルミート多項式級数で表現することもできる。
【0110】
エルミート多項式の直交特性を用いると、m≦Nに対してum(x、t)に比例する項を維持しながら、2つの無限級数(3.7)及び(3.8)を単純に切り捨てることによって、システムは、これらの形式に対するN次近似を得る。1つの速度モデルは、各格子点をその第1及び第2近傍と連結する19速度のD3Q19立方格子である。具体的には、D3Q19(又はD3Q15)型の格子速度セットに対して、システムは(3.6)のfieq(x,t)をu3(x,t)まで、(3.5)のCi(x,t)をu(x,t)の第1(線形)累乗まで展開する。
【0111】
上記の方程式3.9に示すように、この衝突演算子は、19速度をサポート可能なシミュレーションシステムに対して正確である。さらに、この衝突演算子は、19速度をサポート可能なシミュレーションに対して物理的世界で起こるものに対応するように、(3.5)のCi(x,t)を第1(線形)累乗まで展開した形式である。この衝突演算子は、方程式1.12で使用して分布関数を修正することができる。上記の方程式3.9において、xは容積内の特定位置、tは時間における特定時点、iは格子速度セットでの格子速度の指数、T0は一定格子温度、ciは衝突前の粒子速度ベクトル、u(x,t)は時間t及び特定位置xにおける各粒子間の平均速度、Iは2階の単位テンソル、τは衝突緩和時間、wiは一定重み係数、及びΠneqは非平衡運動量流束である。
【0112】
別の速度モデルはD3Q39格子であり、これは39速度までサポートする。D3Q39に関して、展開は、u(x,t)の2乗又は3乗まで行うことができる。u2(x,t)まででの切り捨ては、以下に明示する通りである。
【0113】
上記の方程式3.10に示すように、衝突演算子の展開形式は、2乗を超える、従ってシミュレーションでサポートされない高次のモーメントをその展開形式中に含まないことにより、不必要なモーメントを除外する。本実施例では、xは容積内の特定位置、tは時間における特定時点、iは格子速度セットでの格子速度の指数、T0は一定格子温度、ciは衝突前の粒子速度ベクトル、u(x,t)は時間t及び特定位置xにおける各粒子間の平均速度、Iは2階の単位テンソル、τは衝突緩和時間、wiは一定重み係数、及びΠneqは非平衡運動量流束である。
【0114】
従って、システムは、(3.5)からu(x,t)と無関係の項だけを含むので、(2.4)の衝突演算子形式を0次近似(つまり、(3.5)のCi(0)(x,t))として再解釈する。D3Q39等のより高次の格子速度セットに関しては、u5(x,t)までの平衡分布関数(3.6)の項を維持できるが、u3(x,t)までの非平衡分布関数(3.5)の項は維持できない。一般に、格子速度セットが流体力学モーメントに十分に対応する次数のサポートを提供する場合、この体系的な手順は、任意次数(u(x,t)の累乗)まで行うことができる。所与の有限次数においてガリレイ不変性は厳密には満たされないが、より高い格子速度セットを使用する場合、及び(3.5)及び(3.6)においてより高次の展開形式を用いる場合、誤差はより高次に向かう。
【0115】
上記では、流体速度の累乗での任意の所定の次数に対する一般化されたフィルタ衝突演算子を構成するために、システムは、体系的な手順を生成して実行する。特に、衝突後非平衡流束に関しては、既存のフィルタ演算子(2.4)に対する1次及び2次の補正は、(3.9)及び(3.10)に明示的に示される。フィルタ衝突演算子は、所望のモーメント以外の非平衡モーメントをフィルタ除去する。フィルタ演算子(2.4)及び(3.9)(又は(3.10)は、非平衡運動量流束を維持する目的を果たす。一方で、(3.9)(及び(3.10))は、BGK及び(2.4)と同じ非平衡運動量流束を提供するが、この衝突演算子は、改善された数値的安定性及びガリレイ不変性を実現する。
【0116】
同じ手順をエネルギ流束に関係するフィルタ衝突演算子に対しても策定することができる。一般的な完全ガリレイ不変の形態は、(3.5)に類似しており、相対速度で表される。
ここで、Wneq(x,t)は、(3.1)及び(3.2)のQneq(x,t)及びu(x,t)Πneq(x,t)の適切な線形結合である。同じ手順に従って、所定の格子速度セットによって十分にサポートされるu(x,t)の累乗での任意の有限次数の形式は、体系的に得ることができる。運動量流束とエネルギ流束との間のモーメント直交性により、所望の粘性率及び熱拡散率を独立して実現できる一般的な衝突後演算子形式は、(3.5)及び(3.11)の追加として(適切な展開形式で)簡単に生成される。上記の方程式(3.11)では、xは容積内の特定位置、tは時間における特定時点、iは格子速度セットでの格子速度の指数、T0は一定格子温度、Iは2階の単位テンソル、τは衝突緩和時間、ci’(x,t)は粒子相対速度、fieqは平衡分布関数、及びWneqは非平衡エネルギ流束である。
【0117】
図12を参照すると、本開示に適合するシステムは、非平衡衝突後分布関数を決定するために処理1200を実行する。作動時、本システムは、粒子速度の特定の次数まで流体力学運動をサポートする格子速度セットを提供する(そうでない場合には取得する)(1202)。例えば、本システムは、1次までの流体力学運動の19速度をサポートするD3Q19を取得する。本実施例では、格子速度セットに関するサポートされる次数(例えば19速度)は、非平衡衝突後分布関数の特定の次数より小さくこれとは異なり(例えば、1次又は線形)、非平衡衝突後分布関数の特定の次数は、粒子速度の次数によって決定される。すなわち、19速度をサポートする速度モデルに対して、非平衡衝突後分布関数は線形であるか又は1次までサポートされるだけである(19速度は、1次の流速(つまり、u(x,t))と関係するので)。本実施例において、本システムは、特定の速度モデルに対して流速のサポートされる次数にアクセスするように構成される。本実施例において、本システムは、この情報を有する一覧表又は別の関数にアクセスする。
【0118】
1つの実施例では、格子速度セットは、格子ボルツマン法と関係する状態ベクトルセットである。本実施例において、状態ベクトルは、格子部位での粒子間の衝突挙動を表す一連の2進ビットである。別の実施例において、格子速度セットは、格子に限定される空間で運動量状態のセットから成る。
【0119】
本システムは、格子速度セットにおいて、移動が粒子間の衝突を引き起こす、流体容積中の粒子の移動をシミュレートする(1204)。シミュレーション過程は前述の通りである。シミュレートされた移動に基づいて、容積内の特定位置における粒子の相対速度を決定し(1206)、粒子相対速度は、(i)容積内の特定位置において容積の流れが0の下で測定される粒子の絶対速度と(ii)容積内の特定位置における1又は2以上の粒子の平均速度との間の差である。1つの実施例では、粒子相対速度は、容積内の特定位置において容積の流れが0の下で測定される粒子の絶対速度から容積内の特定位置における1又は2以上の粒子の平均速度を差し引いたものを含む。1つの実施例では、容積内の特定位置における1又は2以上の粒子の平均速度は、その特定位置における特定タイプの粒子の平均速度から成ることが含まれる。例えば、流体容積には多様な異なるタイプの粒子が含まれる場合がある。本実施例において、本システムは、少なくとも特定タイプの粒子から成る部分セットの平均速度を決定するように構成される。前述のように、本システムは、式(3.3)に基づいて粒子の相対速度を決定する。
【0120】
また、本システムは、粒子の相対速度に基づいて、衝突を表す特定の次数の非平衡衝突後分布関数を決定する(1208)。1つの実施例では、非平衡衝突後分布関数は、(i)予め定義された物理量に対する非平衡モーメントを維持し、(ii)未定義の物理量に対する非平衡モーメントを特定の次数まで除去する。本実施例では、非平衡衝突後分布関数は、非平衡衝突後分布関数の展開形式においてこれらの予め定義された物理量を表す項を含むことにより、予め定義された物理量に対する非平衡モーメントを維持する。非平衡衝突後分布関数は、その展開形式において未定義の物理量を表す項を含まないことにより、例えば、方程式(3.7)、(3.8)に示される無限級数を切り捨てて特定の次数に比例する項を維持することにより、未定義の物理量に対する非平衡モーメントを特定の次数まで除去する。本実施例では、特定の次数は、流体速度の格子音速に対する比に関係した指数値であり、格子速度セットは、指数値をサポートする。
【0121】
図13は、ネットワーク環境1300の構成要素のブロック図である。ネットワーク環境1300にはさらにシステム1302が含まれ、そのシステムにはメモリ1304、バスシステム1306、及びプロセッサ1308が含まれる。メモリ1304は、ハードドライブと、ダイナミック・ランダム・アクセス・メモリ、機械可読ハードウェア記憶装置、機械可読媒体又は他タイプの固定機械可読記憶装置などのランダムアクセスメモリとを備えることができる。例えば、データバス及びマザーボードを含むバスシステム1306は、システム1302の構成要素間のデータ通信を確立し制御するために使用することができる。プロセッサ1308は、1又は2以上のマイクロプロセッサ及び/又は処理装置を含むことができる。一般に、プロセッサ1308は、データの受信と記憶、及びネットワーク上の通信(図示せず)が可能な任意の適切なプロセッサ及び/又は論理回路を含むことができる。
【0122】
システム1302は、サーバ、分散型計算システム、デスクトップ型コンピュータ、ラップトップ、携帯電話機、ラックマウント・サーバなどのデータ受信が可能な種々の計算装置のいずれでもよい。システム1302は、単一のサーバとすること、或いは同じ場所に又は別の場所にある一群のサーバとすることができる。図示するシステム1302は、入力/出力(「I/O」)インタフェース1310を経由してデータを受信することができる。I/Oインタフェース1310は、イーサネット(登録商標)、無線ネットワーク・インタフェース、光ファイバ・ネットワーク・インタフェース、モデムなどの、ネットワーク上のデータを受信可能な任意タイプのインタフェースとすることができる。システム1302は、速度モデル、シミュレーションデータなどを記憶するように構成することができるデータ格納部1312との通信のために設定される。
【0123】
本明細書に説明する技術を用いて、システムは、非平衡衝突後分布関数、例えばガリレイ不変フィルタ演算子を生成するように説明される。これらの技術を用いて、種々のタイプの非平衡衝突後分布関数は、例えば、方程式(3.5)、(3.9)、(3.10)及び(3.11)に示すように生成される。生成された非平衡衝突後分布関数は、流体容積中の粒子衝突過程をモデル化する際に、例えば方程式(1.12)に示すように使用する
【0124】
実施形態は、デジタル電子回路、又はコンピュータハードウエア、ファームウェア、ソフトウェア、又はその組合せに実施することができる。本明細書に説明する技術の装置は、プログラマブルプロセッサによる実行のために機械可読媒体(例えば、ストレージデバイス)に有形に具現化又は記憶されたコンピュータプログラム製品に実施することができ、方法アクションは、入力データに対して演算して出力を生成することにより、本明細書に説明する技術の演算を実行する命令のプログラムを実行するプログラマブルプロセッサによって実行することができる。本明細書に説明する技術は、データ及び命令をデータストレージシステム、少なくとも1つの入力デバイス、及び少なくとも1つの出力デバイスから受信して、データ及び命令をそれらに送信するように結合された少なくとも1つのプログラマブルプロセッサを含むプログラマブルシステム上で実行可能である1つ又はそれよりも多くのコンピュータプログラムで実行することができる。各コンピュータプログラムは、高レベル手順又はオブジェクト指向プログラミング言語、又は必要に応じてアセンブリ又はマシン語で実行することができ、いずれにしても、言語は、コンパイル又は解釈された言語とすることができる。
【0125】
適切なプロセッサには、一例として、汎用及び専用マイクロプロセッサの両方が挙げられる。一般的に、プロセッサは、命令及びデータを読取専用メモリ及び/又はランダムアクセスメモリから受信することになる。一般的に、コンピュータは、データファイルを記憶する1つ又はそれよりも多くの大容量ストレージデバイスを含むことになり、そのようなデバイスには、内蔵ハードディスク及び取外し可能ディスクのような磁気ディスク、光磁気ディスク、光ディスクが挙げられる。コンピュータプログラム命令及びデータを有形に具現化するのに適するストレージデバイスには、一例として、EPROM、EEPROM、及びフラッシュメモリデバイスのような半導体メモリデバイス、内部ハードディスク及び着脱式ディスクのような磁気ディスク、光磁気ディスク、及びCD−ROMディスクを含む全ての形態の不揮発性メモリが挙げられる。上述のいずれも、ASIC(特定用途向け集積回路)によって補足されるか、又はASICに組み込むことができる。
【0126】
いくつかの実施を説明した。それにも関わらず、本発明の精神及び範囲から逸脱することなく様々な修正を行うことができることは理解されるであろう。従って、他の実施も以下の特許請求の範囲内である。
【符号の説明】
【0127】
1300 ネットワーク環境
1302 システム
1304 メモリ
1306 バスシステム
1308 処理装置
1310 インタフェース
1312 データ格納部
図1
図2
図3
図4
図5A-5B】
図6
図7
図8
図9
図10
図11
図12
図13