Vol. 55, No. 233
2016. 10

◆特集：LBMの最近の進展http://www.htsj.or.jp/dennetsu/dennetsu55-j.html

図１

D3Q27 MRT LBM による多孔体界面乱流の DNS（ Reb  3000 ，格子総数 3.3 億個）

千鳥状に配置された立方体群からなる空隙率 0.71 の多孔体底面を持つ，チャネル乱流における速度勾配テ
ンソルの第２不変量の II   0.05 等値面（左図）から，多孔体界面での渦生成が圧倒的に滑面（計算領域上面）
より活発であることが分かる．また，断面内の圧力変動と変動速度ベクトルの表示（右図）からは，多孔体内部では乱れが急激に減衰するものの，圧力変動は多孔体層下部まで及ぶことが分かる 1,2)．
1) Y. Kuwata, K. Suga, Int. J. Heat Fluid Flow, (2016), DOI:10.1016/j.ijheatfluidflow.2016.03.006
2) Y. Kuwata, K. Suga, Flow Turb. Combust., (2016), DOI:10.1007/s10494-016-9759-9

図２ D3Q27 MRT LBM によるカーボンペーパー界面乱流のＬＥＳ（ Reb  3000 ，格子総数 0.8 億個）
底面が空隙率 0.8 のカーボンペーパーからなるチャネルの界面（左図： II   0.075 の等値面を瞬時速度 u で表示）でも図１で示された乱流特性が確認できる．また，カーボンペーパー内部（右図）では，不規則な構造と流速分布との関係が分かる．
（特集記事「格子ボルツマン法による乱流の数値解析（須賀一彦，桑田祐丞）」より）

伝熱2016 年 10 月

J. HTSJ, Vol. 55, No. 233

Vol.55

2016

No 233October伝目熱次〈巻頭グラビア〉
D3Q27 MRT LBM による多孔体界面乱流の DNS とカーボンペーパー界面乱流の LES········································ 須賀 一彦（大阪府立大学），桑田 祐丞（東京理科大学） ··········· 表紙裏〈特集：LBM の最近の進展〉
二相系格子ボルツマン法による混相流シミュレーション ··········· 吉野 正人（信州大学） ···········

1

二相系格子ボルツマン法の最近の研究動向 ····························· 稲室 隆二（京都大学） ···········

8格子ボルツマン法による乱流の数値解析 ·························· 須賀 一彦（大阪府立大学）
桑田 祐丞（東京理科大学） ··········· 14
GPU スパコンによる大規模 LBM 計算 ···························· 青木 尊之（東京工業大学） ··········· 22〈行事カレンダー〉 ········································································································· 29〈お知らせ〉
平成 28 年度 日本伝熱学会 学術賞・技術賞・貢献賞・奨励賞候補者推薦募集のお知らせ ··········· 31
日本伝熱学会主催

第4回

国際伝熱フォーラム ································································ 32

日本伝熱学会主催講習会「計測技術」 ··············································································· 33
第 54 回日本伝熱シンポジウム研究発表論文募集·································································· 34
優秀プレゼンテーション賞（第 54 回日本伝熱シンポジウム）について ··································· 36
事務局からの連絡・編集出版部会からのお知らせ······················································································ 37
・新入会員一覧 ··········································································································· 38〈編集出版部会ノート〉 ·································································································· 41

Vol.55，No. 233，October 2016
CONTENTS＜Opening-page Gravure:heat-page＞
Kazuhiko SUGA (Osaka Prefecture University)
Yusuke KUWATA (Tokyo University of Science) ······································

Opening Page

＜Special Issue: Recent Development of LBM＞
Numerical Simulation of Multiphase Flow by the Two-Phase Lattice Boltzmann Method
Masato YOSHINO (Shinshu University) ································································· 1
Recent Trend of Research on Two-Phase Lattice Boltzmann Method
Takaji INAMURO (Kyoto University) ··································································· 8
Turbulent Flow Simulation by the Lattice Boltzmann Method
Kazuhiko SUGA (Osaka Prefecture University)
Yusuke KUWATA (Tokyo University of Science) ····················································· 14
Large-scale Computations of Lattice Boltzmann Method on a GPU Supercomputer
Takayuki AOKI (Tokyo Institute of Technology) ····················································· 22＜Calendar＞ ···················································································································· 29＜Announcements＞ ·········································································································· 31＜Note from the JHTSJ Editorial Board＞ ········································································· 41

特集：LBM の最近の進展

二相系格子ボルツマン法による混相流
シミュレーション
Numerical Simulation of Multiphase Flow by the
Two-Phase Lattice Boltzmann Method
吉野 正人（信州大学）
Masato YOSHINO (Shinshu University)e-mail: masato@shinshu-u.ac.jp

格子ボルツマン法[1-3]（Lattice Boltzmann Method,
以下 LBM と記す）は，1990 年代以降開発が進んできた非圧縮性粘性流体の比較的新しい数値計算法
である．LBM では，気体分子運動論のアナロジーを利用して，流体を有限個の速度をもつ多数の仮想粒子の集合体（格子気体モデル）で近似し，各粒子の衝突と並進とを粒子の速度分布関数を用いて逐次計算し，そのモーメントから巨視的な流れ場や温度（濃度）場を計算する．LBM の特長は，複雑流れに対してもアルゴリズムが簡単であることや，並列計算に向いていることなどがあげられる．そのため，これまでに多孔質内流れなどの複雑な境界をもつ流れや，気液・液液二相流などの界面が複雑に変形する流れに適用され成功を収めている．
本稿では，特に微小スケールの混相流現象（液滴の挙動や変形する物体を含む流れ）に焦点を当て，筆者らの最近の研究成果を紹介する．
２．二相系流れの LBM モデル
これまでに提案されている二相系 LBM のモデル
は，大別すると (i) Color-Gradient モデル [4]， (ii)Pseudopotential モデル[5]，(iii) Free-Energy モデル（Phase-Field モデル）[6]に分類される[7]．このうち(iii)のモデルは，非平衡熱力学に基づく系の自由エネルギーを導入した方法であり，二相の界面形状
は系の自由エネルギーが最小となるように自律的に変形して決まる．Swift ら[6]が提案した最初のモ
デルは，数値安定性の問題から二相流体の密度比に限界（せいぜい 10 程度くらいまで）があり，高密度比の二相流を扱える LBM の開発が望まれていた．
これに対し Inamuro ら[8]は，Free-Energy モデルをベースに Projection 法[9]を導入した高密度比の二相系 LBM を開発し，本手法を用いることにより，例
えば水－空気系に代表される密度比 1000 程度の計算を安定に実行することができると報告している．

伝熱2016 年 10 月

-1-D

１．はじめに

図1

固体壁面上における液滴衝突

以下では，3.4 節を除き，この高密度比の二相系
LBM を用いて計算を行った．LBM の速度モデルとしては，3.1 節では 2 次元 9 速度モデル[10]，それ以外では 3 次元 15 速度モデル[10]を用いた．なお，
特に断らない限り，変数は適切な代表量（基準量）
を用いて無次元化されたものである．計算手法および定式化の詳細は，文献[8]を参照されたい．
３．数値計算例
3.1 固体壁面上での液滴の衝突・混合現象の解析
固体壁面上における液滴の衝突・混合現象は，例
えばスプレーコーティングなどで見られるように，工業的にも重要な現象の一つである．ここでは，固体壁面上で静止する液滴に別の液滴が衝突・混合する際の動的な挙動について，二次元計算を行った結果[11]を述べる．
図 1 に示すようなの計算領域に対して，0の固体壁面上に静止する直径 の液滴（液滴 A）に，同じ大きさの別の液滴（液滴 B）が上方から衝突する場合を考える．重力は 方向に作用する（重力加速度を とする）．
0の固体壁には濡れ性を考慮したすべりなし境界条件，
では各物理量の勾配が 0 のノイマン条件，左右の境界には周期境界条件を用いた．本問題の無次元パラメータは，液相密度 と気相密度 の比（ ⁄ ），液相粘度 と気相粘度 の比（ ⁄ ），液相基準のレイノルズ数⁄ ，液相基準のウェーバー数⁄ （ :界面張力），およびフルード数

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

図2

液滴形状の時間変化（

10,

⁄2）

図3

⁄

Mixing rate [%]55.6（
2.27 10 ∆ ，
4.08 10 ∆ ），∆1.07 10 ,3.32 10 （∆ ：空間の格子間隔）であり，その他の値は文献[11]を参照していただきたい．
結果の一例として，固体壁面の静的接触角を⁄2とし，10および100に対する液滴形状の時間変化をそれぞれ図 2, 図 3 に示す（各図は，計算領域の下半分 0 y/2 のみを表示している）．
これらの図は，液滴 A と B を区別するために，二色のトレーサ粒子をそれぞれの液滴に埋め込み，その挙動を四次のルンゲ・クッタ法によって計算した結果を表示したものである．
10のケース（図 2）では，液滴 B が液滴 A の上に衝突した際に，飛び散ることなく壁面上に広がる結果が得られている．一方，100のケース（図 3）では，衝
突後二つの液滴が合一した両端部分にリムと呼ばれる特徴的な突起部が形成されている．また，合一した液滴は壁面近傍で分裂し，ほぼ液滴 A で構成

2016 年 10 月

⁄2）

55.0

であるが，一般に，液滴の直径は 2～3mm 程度であることを考慮すると，ウェーバー数が支配パラメータであり，その違いによって液滴の挙動は大きく異なると考えられる．本計算で用いた主なパラメータの値は，
1500∆
450∆ ,
50，1)， /100∆ ， ⁄
50 (

伝熱

100,

液滴形状の時間変化（

-2-

52.5
50.0
47.5
45.0
0

20

40

60

80

100

We
図4

ウェーバー数と合一液滴の混合率の関係された微小液滴が生成されている．その後，この微小液滴が飛散していく間に，合一した液滴は A とB の両方の液が混合し再び分裂が発生する．
次に，合一液滴における混合率（mixing rate）を各ウェーバー数に対して計算した結果を図 4 に示
す．ここで，混合率は合一液滴における落下液滴（液滴 B）の割合を表している．本図から，
40では混合率はウェーバー数とともに増加するが，40
100ではほぼ一定値になっている．言い換えれば，飛沫した液滴中に含まれる A と B の割合は，40100ではほぼ一定となる結果が得られた．現在，本問題の三次元解析に取り組んでいるところである．

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

図5

液滴衝突の計算領域および変数の定義

図 7 液滴衝突後における挙動の時間変化 ()
2000,
40, ∆ 0.5, ∗ 2 ⁄

図6

同径液滴の衝突におけるレイノルズ数と臨界ウェーバー数の関係

3.2 液滴同士の二体衝突解析
液滴同士の衝突は，機械工学ではエンジン内の燃
料噴霧，自然界では雨滴の形成など様々な分野で見られる現象である．本研究では，種々の条件下における液滴同士の二体衝突解析を行った[12]．以下ではその結果について述べる．
図 5 に示すようなの計算領域内に，

図8）の二つの液滴を配置し，相直径が ， （
対速さ で衝突させた．境界条件は，全面に周期境界条件を適用した．この系の支配パラメータは，液相密度 と気相密度 の比（ ⁄ ），液相粘度 と気相粘度 の比（ ⁄ ），液相基準のレイノルズ数⁄ ，液相基準のウェーバー数⁄ ，および 2 つの液滴の直径比∆⁄ である．以下では0.4 ∆ 1.0とし，全領域を0.4 ∆
0.7 では0.7 ∆ 1.0 で は

246∆

3.2.1 同径液滴の衝突挙動解析
同じ大きさの液滴衝突（）において，各レイノルズ数に対してウェーバー数を変化させ計算を行った．その結果，衝突後の液滴の挙動は，同じレイノルズ数でもウェーバー数の値によって合体するケースと分離するケースが得られた．本研究では，ウェーバー数を 1 刻みで変化させ，衝突後の挙動を合体と分離の 2 つのパターンに分類し，そのとき挙動が変化する前後 2 点のウェーバー数の中間値を臨界ウェーバー数と定義した．各レイノルズ数に対する衝突挙動の形態と臨界ウェーバー数を図6 に示す．図中の実線は挙動形態を分類する境界線であり，Volume of Fluid（VOF）法[13]に基づく数値計算から得られた Saroka らによる結果[14]である．
この図より，
2000ではレイノルズ数が小さい

128∆
128∆ に，
192∆
96∆

96∆ にそれぞれ分割した．無次元パラメータを⁄⁄
50,
50, 600
4000, 2080とし，計算パラメータを

32∆ ,

50，
1, 32∆
80∆ , 4.00
0.1,
10 ∆
2.67 10 ∆ , 8.00 10 ∆5.33 10 ∆ ，0.2∆
σ 0.8∆ とした．

伝熱2016 年 10 月

液滴の直径比と臨界ウェーバー数の関係-3-

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

t* = 0.00

t* = 1.56zt* = 2.03

液滴挙動の時間変化（xt* = 3.75

4.97,

0.7）

t* = 2.03

t* = 3.75

図 10 液滴挙動の時間変化（

ほど合体しやすく，臨界ウェーバー数は増加していることがわかる．これはレイノルズ数の減少に伴い液相の粘度が相対的に増加するため，液滴同士が合体しやすくなったことが要因と考えられる．一方，
2000では，臨界ウェーバー数はレイノルズ数に依存せずほぼ一定の値を示し，Saroka らによる結果[14]と良く一致している．本結果より，
2000
では，レイノルズ数は衝突挙動に影響を与え，臨界ウェーバー数を支配するパラメータであると言える．
3.2.2 異径液滴の衝突挙動解析
異なる大きさの液滴衝突において，直径比およびレイノルズ数が衝突挙動に与える影響について調べた．結果の一例として，
2000，
40，∆
0.5のケースに対する液滴挙動の時間変化を図 7 に示す．液滴は衝突直後に半径方向に広がって傘状になり，その後界面張力により縮み，初速方向に伸び始める．瓢箪状になった液滴はさらに伸び続け，やがては二つの液滴に分離する結果となった．
次に，
600, 1000, 2000, 4000に対して同様の計算を行い衝突の際の挙動を調べた．各レイノルズ数に対して，液滴の直径比と臨界ウェーバー数の関係を図 8 に示す．この図より，いずれのレイノルズ数においても，直径比が∆ 0.7前後で臨界ウェーバ
ー数は増加しており，∆ 0.7のときに最小となることがわかる．よって，600
4000では，直径比に対する臨界ウェーバー数の依存性は，いずれのレイノルズ数に対しても同じであることがわかった．
また，直径比が0.55 ∆ 0.9では，レイノルズ数が小さいほど臨界ウェーバー数は増加しており，この伝熱

t* = 1.56zx

図9t* = 0.00

2016 年 10 月

0.994，

0.7）

傾向は同径液滴の場合と一致している．一方，∆0.55では，同径液滴の傾向とは異なることから，この関係は異径特有の結果であると言える．
3.3 雲の成長過程における微小水滴の挙動解析
雲は，雲粒子と呼ばれる直径数 mから数 mm サ
イズの非常に小さな水滴の集まりで構成され，それらが重力や上昇気流等の影響で衝突し，合体や分裂を繰り返すことで粒径が変化し，雲の成長や消滅を促すと言われている[15]．それゆえ，雲の成長や消滅現象を解明するためには，微小水滴の衝突挙動を調べることが重要である．
以下では，雲中における微小水滴（雲粒子）の挙動を計算するために，直径が5 m
50 m，相対衝突速さが0.1m/sの雲粒子を想定した衝突挙動を調べた[16]．3.2 節で述べた図 5 のような計算領域内に，直径 の同一液滴を二個配置し，相対
速さ で衝突させた．主な計算パラメータは，
128∆
64∆
96∆ ,
32∆ ,
5，
1）， ⁄
5（
1.0 10 , ⁄3.21 10 ∆ , 5.83
55 , 3.21 10 ∆10 ∆5.83 10 ∆ ， σ 1.16 10 ∆とした．ここで，液相と気相の密度比について，実際には 1000 程度であるが，液滴のサイズが非常に小さく重力（浮力）の影響があまり大きくないと予想されるため，本研究では ⁄5とした．ウェーバー数については，雲の場合10 ～10 のオー
ダであるが，数値安定性の問題から1.37とし，
液滴直径および相対衝突速さ基準のレイノルズ数⁄ を0.497
4.97で変化させた．
0.5∆ とし
また，文献[8]のモビリティの値は，-4-

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

図 12 狭窄部を過ぎる物体の変形挙動解析
図 11

レイノルズ数と衝突パラメータに対する液滴の挙動形態

た（ただし，その値が最適かどうかについては，今後検討の必要がある）．
結果の一例として，衝突パラメータ（液滴の初期直径に対する無次元オフセット量）を0.7とし，
4.97, 0.994（それぞれ，
50 m，10 mに相当する）に対する液滴の挙動をそれぞれ図 9，図 10に示す．ここで，二つの液滴には，時刻 ∗ 0で 方向の初期速度を与えた．図 9 のケースでは二つの液滴が衝突して合体しているが，図 10 のケースでは液滴同士が接近した後，それる挙動を示し衝突は起
こらない結果となった．
次に，種々のパラメータに対して同様の計算を行い，液滴挙動の違いを表したダイアグラムを図 11に示す．液滴の挙動は 2 つのパターン（合体：coalescence，それる：deviation）に分類されることがわかる．特に2では，レイノルズ数が小さくなるにつれて，衝突パラメータが小さいケースでも液滴はそれる（衝突しない）ことがわかる．これは，液滴サイズがさらに小さくなり，その軌道がまわりの流れ場の影響を受けやすくなるためだと考えられる．現在，種々の条件下での計算を実行し，液滴が受ける揚力を調べるなど詳細な検討を行っているところである．
3.4 粘弾性変形する物体を含む固液二相流解析
変形を伴いながら移動する物体を含む流れの問
題は，医学や工学をはじめとする多くの分野で見ることができる．例えば血液の流れは，赤血球などの固体成分が血漿とともに流れる固液二相流である．
赤血球は，粘弾性をもつ膜が内部流体を覆った直径約8 mの粘弾性皮膜物体（以下では単に物体と記す）であり，柔軟に変形するため，それより小さい血管径の流路においてもスムーズに流動することがで伝熱2016 年 10 月

-5-

図 13 狭窄部で変形する物体の形状変化
きる[17]．筆者らは，これまでに，質点－ばねモデルに基づく弾性力を等密度の液液二相系 LBM[18]に組み込んだ計算手法[19]を構築した．ここでは，赤血球を模擬した双凹面形状の物体が，狭窄部を通過する際の挙動について調べた結果[20]を述べる．
図 12 に示すような中央部に狭窄部をもつ円管を考える．内部に配置された物体は，赤血球を模擬した双凹面形状をしており，外部から力を受けると粘弾性変形する．領域の入口と出口には，圧力差を伴う周期境界条件を適用した．また，円管の表面には，埋め込み境界法[21] に基づくすべりなし境界条件を用いた．主な計算パラメータは，
105∆ ，
40∆ ，
10∆ ，
49.6∆ ，
22∆ ，
1，
3.93
19.2∆ ，
10 ∆ ,
9.42 10 ∆ である．ここで，および はそれぞれ，周囲流体および内部流体の密度， ,,
はそれぞれ，周囲流体，内部
流体，および膜の粘度である．その他の値は文献[20]を参照されたい．なお，重力は考慮していない．
結果の一例として，狭窄部で変形する物体の形状変化を図 13 に示す．物体は，狭窄部を通過する際に流れ方向（ 方向）に大きく変形し，上流側が凹面，下流側が凸面のパラシュート形状になっている．
次に，物体の変形度および 方向速度 の時間変化をそれぞれ図 14, 図 15 に示す．ここで，物体の, , 方向の長さをそれぞれ , , とすると，変形度は2 /で定義した．また
は，物体が存在しない場合のポアズイユ流における最J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

0.25

1.6

0.20

1.2
0.15

0.8
0.10

0.4
0.0
0

0.05

5

10

15

20

0.00

25

0

5

10

15

20

25

図 14 物体の変形度の時間変化

図 15 物体の 方向速度 の時間変化

大流速である．これらの図より，物体が狭窄部を通過する際の流れ方向（ 方向）の長さは，両スパン
方向（ および 方向）の平均長さの約 1.5～1.6 倍であることがわかる．また，物体は狭窄部に入るときに若干加速し，狭窄部内では減速し，狭窄部から出たときに再び大きく加速することもわかる．
今後の課題としては，複数個の物体の挙動解析を行い，狭窄部における物体の閉塞現象を調べることなどがあげられる．

最後に，最近 Inamuro ら[23]は，高密度比の二相流に対しても圧力のポアソン方程式を解く必要のない，画期的な改良二相系 LBM を提案した．この改良二相系 LBM を用いることにより，高密度比の二相流の場合でも数値的に安定で，かつ，計算時間もかなり短縮されることが期待できる．今後は，新しい手法を用いて，今まで計算が難しかった条件や新たな複雑流れの問題にもチャレンジしていきた
いと考えている．

４．おわりに本稿では，二相系 LBM を用いた微小スケールの混相流問題の解析結果をいくつか紹介してきた．二相系 LBM は，VOF 法[13]などの従来の混相流計算法に比べて，各相の質量および運動量の保存性に優れており，界面をシャープにとらえることができる
効率の良い数値計算法である．特に最近では，計算機環境の目ざましい発達のおかげで，LBM を用いた大規模並列計算や GPU による高速計算（例えば文献[22]など）が行われるようになってきた．したがって，今後ますます複雑な流れ問題への LBM の適用が期待できる．
その一方で，解決しなければいけない課題もまだ残されている．例えば，3.3 節で述べた雲の成長過程における微小水滴の挙動解析では，レイノルズ数を一定にしたまま非常に小さなウェーバー数（10～10 ）を実現するためには，界面張力をかなり大きく（現状の計算の 100～1000 倍程度大きく）する
必要があり，数値安定性の問題から計算は実行できていないのが現状である．また，現在の二相系 LBM
では，界面において非圧縮性流体の連続の式を満たすために圧力のポアソン方程式を解く必要があり，その高速解法が課題である．

謝 辞
本稿で紹介した数値計算は，信州大学大学院修了
生の田中義人 君，村山寿郎 君，勝見真悟 君，深谷昇弘 君，ならびに現在在籍中の大学院生 澤田純
平 君のご協力により行われました．また，本研究の一部は，JSPS 科研費 JP23560192, JP26420105,
JP15H02218 の助成, ならびに学際大規模情報基盤
共 同 利 用 ・ 共 同 研 究 拠 点 （ jh140025, jh150012,
jh160012）の支援を受けたものです．ここに記して謝意を表します．

伝熱2016 年 10 月

-6-

参考文献
[1] Succi, S., The Lattice Boltzmann Equation for Fluid
Dynamics and Beyond, Oxford University Press(2001).
[2] Inamuro, T., Lattice Boltzmann Methods for Viscous
Fluid Flows and for Two-Phase Fluid Flows, Fluid
Dyn. Res., 38 (2006) 641.
[3] Aidun, C. K. and Clausen, J. R., Lattice-Boltzmann
Method for Complex Flows, Annu. Rev. Fluid Mech.,
42 (2010) 439.
[4] Gunstensen, A. K. et al., Lattice Boltzmann Model of
Immiscible Fluids, Phys. Rev. A, 43 (1991) 4320.

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展[5] Shan, X. and Chen, H., Lattice Boltzmann Model for
Simulating Flows with Multiple Phases and[15] 大西領，高橋桂子，雲に見られる乱流現象－気相乱流中での微小水滴の衝突，ながれ，30 (2011)385.
[16]Yoshino, M. et al., Two-Phase Lattice Boltzmann
Simulation of Binary Collisions of Small Water
Droplets in Clouds, Abstracts of International
Workshop on Cloud Turbulence, Nagoya (2015).

Components, Phys. Rev. E, 47 (1993) 1815.
[6] Swift, M. R. et al., Lattice Boltzmann Simulation of
Nonideal Fluids, Phys. Rev. Lett., 75 (1995) 830.
[7] Chen, L. et al., A Critical Review of thePseudopotential Multiphase Lattice Boltzmann
Model: Methods and Applications, Int. J. Heat Mass

[17]菅原基晃，前田信治，血液のレオロジーと血流，
コロナ社 (2003).
[18]Inamuro, T. et al., Lattice Boltzmann Simulations ofDrop Deformation and Breakup in Shear Flows, Int.

Transfer, 76 (2014) 210.
[8] Inamuro, T. et al., A Lattice Boltzmann Method forIncompressible Two-Phase Flows with Large
Density Differences, J. Comput. Phys., 198 (2004)
628.
[9] Chorin, A. J., Numerical Solution of the Navier-J. Mod. Phys. B, 17 (2003) 21.
[19]Murayama, T. et al., Three-Dimensional Lattice
Boltzmann Simulation of Two-Phase Flow
Containing a Deformable Body with a Viscoelastic

Stokes Equations, Math. Comput., 22 (1968) 745.
[10] Qian, Y. H. et al., Lattice BGK Models for Navier-Membrane, Commun. Comput. Phys., 9 (2011) 1397.
[20]Yoshino, M. and Katsumi, S., Lattice Boltzmann
Simulation of Motion of Red Blood Cell in

Stokes Equation, Europhys. Lett., 17 (1992) 479.
[11] Tanaka, Y. et al., Numerical Simulation of Dynamic
Behavior of Droplet on Solid Surface by the TwoPhase Lattice Boltzmann Method, Comput. Fluids,

Constricted Circular Pipe Flow, J. Fluid Sci. Tech., 9(2014) 14-00133.
[21]Peskin, C. S., Flow Patterns around Heart Valves: A

40 (2011) 68.
[12] 澤田純平ら，二相系格子ボルツマン法を用いた
異径液滴衝突の数値解析，混相流シンポジウム
2016 講演要旨集, D323 (2016).
[13] Hirt, C. W. and Nichols, B. D., Volume of Fluid(VOF) Methods for the Dynamics of Free

Numerical Method, J. Comput. Phys., 10 (1972) 252.
[22]Feichtinger, C. et al., Performance Modeling andAnalysis of Heterogeneous Lattice Boltzmann
Simulations on CPU-GPU Clusters, Parallel
Comput., 46 (2015) 1.
[23]Inamuro, T. et al., An Improved Lattice Boltzmann
Method for Incompressible Two-Phase Flows with

Boundaries, J. Comput. Phys., 39 (1981) 201.
[14] Saroka, M. D. et al., Numerical Investigation of
Heat-on Binary Drop Collisions in a Dynamically

Large Density Differences, Comput. Fluids, 137(2016) 55.

Inert Environment, J. Appl. Fluid Mech., 5 (2012) 23.

伝熱2016 年 10 月

-7-

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

二相系格子ボルツマン法の最近の研究動向
Recent Trend of Research on Two-Phase Lattice Boltzmann Method

稲室 隆二（京都大学）
Takaji INAMURO (Kyoto University)e-mail: inamuro@kuaero.kyoto-u.ac.jp
に関する国際会議の一つである今年 7 月に深セン（中国）で開催された DSFD2016 [21]でも，多数の二相系 LBM の発表があった．
以下では，高密度比の気液二相流への適用を目指した二相系 LBM の最近の研究動向を紹介する．特に，高密度比の二相系 LBM の Key となる圧力場の計算方法に注目して，最近の提案手法を解説する．

１．はじめに格子ボルツマン法（Lattice Boltzmann Method, 以
下 LBM）が非圧縮性粘性流体の数値計算法として提案されたのは 1990 年代初めである．世界的に見ても，まだまだ利用者は多いとは言えないが，非圧
縮性粘性流体の数値計算法の一つとして認知されており，その特徴を生かして適用分野がますます広がっている．
さて，LBM とは，流体を有限個の速度（例えば，3 次元では 15 速度，19 速度，27 速度など）をもつ
多数の仮想粒子の集合体（格子気体モデル）で近似
し，各粒子の衝突と並進を表す粒子の速度分布関数の発展方程式を逐次計算し，その速度分布関数のモーメントから巨視的流れ場（流速，圧力など）を求める完全陽的な数値計算法である．LBM の最大の特徴は，圧力のポアソン方程式を解く必要がないた
めに，アルゴリズムが簡単であること，ならびに並列計算に適していることである．格子上の値は隣接する格子点の値のみから求めることができることから，究極の局所的な計算スキームである．そのため，この約 20 年の間に，複雑な流路内の単相流だ
けでなく，種々の移動境界流れ（気液二相流，流体‐構造連成流れなど）へ適用されている．それらの発展の様子については，例えば，拙著の解説記事[14]を参照されたい．
本稿では，本会誌の読者の関心の高い気液二相流に適用できる二相系 LBM を取り上げる．二相系
LBM の研究は，LBM が提案された直後の 1990 年
代半ばからすでに始まっている．当初は密度比<10程度の計算が限界であったが，2000 年代になって
密度比≒1000 の計算が可能な方法が立て続けに提案され[5-7]，広い分野から二相系 LBM に注目が集まった（ただし，その後，Zheng らの方法[7]は等密度比の方法であることが指摘されている）．それ以降 10 年以上経過しているが，高密度比の二相系
LBM の研究がさらに進展している．ちなみに，LBM

伝熱2016 年 10 月

２．二相系 LBM
一般に，二相系 LBM は，二相（気相および液相）
をどう表現するかによって以下の 3 つの方法に大別できる．なお，以下の説明は，すべて無次元変数（定義は文献[1-4]を参照）で行う．また，等温場の二相系 LBM（非熱流体モデル，athermal model と呼ばれる）を取り扱う．
2.1 Color-Gradient Model
このモデルは，もともと格子気体モデル（粒子の速度分布関数を逐次計算するのではなく，格子点上に粒子を置いて衝突と並進を追跡する本来の粒子モデル）から発展してきたもので，Gunstensen ら[8]によって LBM に導入された．このモデルでは，2つの色分けされた粒子を考え，それらの間に適切な相互作用を導入して界面を表現する．二相流の数値計算でよく用いられる VOF 法に近い考え方である．
当初は等密度比の二相流にしか適用できなかったが，その後，高密度比へ発展してきている[9]．しかし，静止液滴や静止気泡の計算では密度比が 10000の計算例もあるが，運動する液滴や気泡の計算では安定して計算できるのはせいぜい密度比が 10 程度
である．
2.2 Pseudo-Potential Model
Shan and Chen [10]によって提案されたもので，その後，多くの論文が発表されている．上記の ColorGradient Model を現象論的な物理モデルで置き換え-8-

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展たもので，粒子間に擬似ポテンシャルを導入することにより界面を表現している．なお，擬似ポテンシ
ャルは，二相を表す状態方程式（例えば，ファン・デル・ワールス方程式）から求めることができる．
この擬似ポテンシャルを粒子の速度分布関数の発展方程式の外力項として追加する．当初の方法では密度比は 10 程度が限界であるが，その後，状態方
程式の形を工夫したり[11]，粒子速度を増やす[12]などにより密度比が 1000 程度の計算もできるようになったが，二相ポアズイユ流の計算では密度比>20 で速度の不連続が発生することも報告されている[12]．
2.3 Free-Energy Model（Phase-Field Model）
Swift ら [13]によって初めて二相系 LBM に導入されたもので，非平衡熱力学に基づいて二相界面を形成するモデルである．粒子の速度分布関数の発展
方程式の平衡分布に圧力テンソルの形で導入する（あるいは，外力項として導入することもできる）．
このモデルは，Phase-Field Model とも呼ばれ，材料科学分野の界面現象の計算に用いられてきたモデルと同じである．やはり，状態方程式（例えば，ファン・デル・ワールス方程式）から気相と液相の密度 ρG および ρL が決まる．界面の表現が理論に基づいているため，他の 2 つのモデルより優れていると言われているが，高密度比の計算では圧力の計算が発散し易いため，安定した圧力の計算に工夫が必要である．前述の文献[5-7]の方法はすべてこのモデルを用いている．最近，密度比≒100 に適用できる方法[14]が提案されたが，密度比>100 では圧力計算が不安定になると言われている．
以上，界面を表現する 3 つのモデルを紹介したが，それぞれのモデルには一長一短がある．ここでは，それらの詳細には立ち入らず，以下では高密度比の二相系 LBM において Key となる圧力計
算について焦点を絞って説明する．
３．圧力計算法
よく知られているように，LBM は弱圧縮性流体を仮定（流れのマッハ数<<1）している．そのため，
非圧縮性流体の連続の式∇・u=0 を満たす圧力 p は，
速度分布関数のモーメントから計算した密度 ρ（連続の式を満たしている）を状態方程式 p=cs2ρ（ここ伝熱2016 年 10 月

-9-で，cs= 1/3は音速）に代入して求めることができる[1-3]．このように，圧力に対するポアソン方程式
を解くことなく圧力 p を求めることができるのが
LBM の最大の利点である．二相系 LBM においても圧力 p はポアソン方程式を解くことなしに求めることができるが，密度比≒1000 の場合には安定した計算が困難になる．
筆者らも 2000 年前後に，この高密度比での数値
不安定性を解決するのに苦労した挙句，取りあえずの解決策としてポアソン方程式を解いて圧力を求めることにした[5]．その後，ポアソン方程式を用いない方法を模索してきて，最近，4 章で紹介する新
しい方法を提案した．この模索の中でわかったことは，高密度比の場合には圧力計算がしばしば不安定になり，圧力 p を計算するには工夫が必要であるということである．すなわち，高密度比の二相系 LBM
では，圧力計算が Key となる．
以下では，二相系 LBM で用いられている，ポア
ソン方程式を用いずに圧力を求める２つの方法を紹介する．
3.1 連続の式および状態方程式
単相流の場合の LBM と同様に，次の連続の式(1)を満たす密度 ρ を速度分布関数のモーメントから計算し，状態方程式（例えば，ファン・デル・ワー
ルス方程式）(2)に代入して圧力 p を計算する．なお，気相と液相の密度 ρG および ρL は，マクスウェルの等面積則から決まることに注意する（図 1 参照）．すなわち，(1)および(2)式から求まる圧力 p は，熱力学的圧力とも呼ばれる．界面内の密度変化ならびに領域内の流速および圧力を１つの速度分布関数を用いて計算することができる．この方法は，単相流の LBM での圧力の求め方と全く同じ方法である．ただし，/ なので，図 1 からわかるように，気相とcs=
液相で音速 cs が異なること，ならびに状態方程式の形によっては音速 cs が小さくなって，LBM の制
限条件である，流れのマッハ数 u/cs<<1 の条件を満たすことができなることに注意しなければならない．

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

3.3 高密度比と数値不安定性
上記の 2 つの方法は原理的には等価であるが，一般に，(3)式を用いる方が安定であることが知られている[15]．二相系 LBM においても，密度比<100では(1)および(2)式で（熱力学的）圧力を計算することができるが，密度比≒1000 では(3)式を用いて（流体力学的）圧力を求める必要がある．最近発表された高密度比に適用できる方法[16]も(3)式を用いて圧力を求めている．そういう観点から，高密度比の圧力計算については，文献[6]を超える方法はまだ提案されていない．(1)および(2)式を用いる圧力計
算の改良が望まれるところである．今後の研究の進展に期待したい．また，ρL=O(1)および ρG<<1 の場
合（これまでに提案されているほとんどの二相系
LBM は，この場合を取り扱っている）には，気相の速度場が大きく変動して u/cs<<1 の条件を満たせなくなり発散することが多い．特に，密度比≒1000の気泡流の計算が困難になることが指摘されている[17]．
したがって，密度比≒1000 では，非圧縮性流体の連続の式∇・u=0 を満たす（流体力学的）圧力 p を求める工夫が必要である．

図 1 ファン・デル・ワールスの状態方程式（温度 T を決めたとき，ρG および ρL はマクスウェルの等面積則から決まる．
）理想気体の状態方程式
p=cs2ρ も表示している．
3.2 圧力発展方程式
一方，(1)および(2)式を変形すると(3)が得られる，この式は圧力発展方程式と呼ばれる．
なお，(3)式の第 3 項は他の 2 項に比べて小さいので無視されることが多い．界面内の密度を計算する
速度分布関数とは別の新しい速度分布関数を考え，(3)式を満足するように新しい速度分布関数の発展
方程式を定式化すると，その速度分布関数のモーメントから直接圧力 p を求めることができる[6]．なお，(3)式中の音速 cs は単相流と同じ 1/3が用いることが多いが，原理的には自由に選ぶことができる．
その時， (3)式から求まる圧力 p は，流体力学的圧力とも呼ばれる．すなわち，圧力 p は二相を表現す
る状態方程式(2)と独立に決まり，界面の表現と圧力との間に熱力学的一貫性がないことになる．このことは，界面内の密度変化ならびに領域内の流速および圧力を計算するのに，2 つの速度分布関数を用いることからも明らかである．これが圧力発展方程式を用いるときの弱点であるが，従来法（VOF 法，
Front Tracking 法，Level Set 法など）でも界面の表
現と圧力計算は切り離されているので，従来法と同じ取り扱いであると考えると，二相流計算法としては成立する．

伝熱2016 年 10 月

- 10 -

４．新しい方法および数値計算例
以下では，筆者らの最新の二相系 LBM とそれを用いた数値計算例[18]を紹介する．
4.1 新しい方法の紹介
新しい方法の詳細については文献[18]を参照されたい．ここでは簡単に方法を紹介する．
界 面 の 表 現 に は 文 献 [5] と 同 様 に Free-Energy
Model を用いている．圧力計算には(3)式の圧力発展式を満たす圧力 p を LBM の発展形である LKS（Lattice Kinetic Scheme）を用いて計算する．すな
わち，2 つの速度分布関数を用いて，界面内の密度変化ならびに領域内の流速および圧力を計算する．
なお，本方法では，気相の密度を ρG=1 ，液相の密度を ρL>>1 としている．したがって，(3)式から分かるように，気相に比べて液相では圧力の時間発展が大きくなる．このことを考慮して，圧力計算における LKS に過緩和係数 ω を導入することのより，両相の音速の値が O(1)になるように調整している．
以下の数値計算例で示すように，新手法では密度比≒1000 の計算が可能である．

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

4.2 数値計算例 1（ミルククラウン）
新しい方法を用いたミルククラウンの計算結果
を図 2 に示す．液滴直径 D を 240Δx，液膜を 30Δ
x とし，密度比=800，レイノルズ数 Re=2000，ウェ
ーバー数 We=1000 とした計算結果である．3 次元計算であるが，液滴中心を通る 2 次元平面上の密度分布を示している．液滴が液膜に衝突後に形成されるクラウンのシート状膜が広がる様子は良くとらえられているが，最後にシート状膜が小さな液滴に分裂する様子は実験結果と少し異なっている．図 3
は，クラウンのシート状膜の半径 rw の時間変化を簡易な理論予測値と比較したものである．計算結果の rw が理論予測値と良く一致していることがわかる．ミルククラウンの計算は，計算科学の長年の夢と言われている難しい問題である．二相系 LBM により，その夢が叶えられることを期待している．
4.3 数値計算例 2（水面に落下する球）
次に，新しい二相系 LBM と埋め込み境界法[19]を組み合わせた方法を用いて，水面に落下する球（表面は親水性）の計算結果[20]を図 4 に示す．球の直径 D=250Δx とし，気液の密度比=800，レイノ
ルズ数 Re=5000，キャピラリー数 Ca=0.033，ウェー
バー数 We=166 としたときの計算結果である．図 4より，落下する球に水面が引きずり込まれ(b-e)，その反動で水柱が形成される様子(f)が確認できる．この問題は，気・液・固の三相が絡む難しい流体問題
であるが，二相系 LBM が適した応用分野である．

図 2 ミルククラウンの計算結果（t*=tU/D，U：液滴の初期速度，D：液滴直径）

図 3 クラウンのシート状膜の半径 rw の時間変化．
実線と破線は簡易な理論予測値を示す．
（tc*：液滴が液面に接触してからの無次元時間）

伝熱2016 年 10 月

- 11 -

図 4 水面に落下する球の計算結果．(a)t*=0, (b)t*=0.28, (c) t*=0.88, (d) t*=1.48, (e) t*=1.96, (f)t*=2.80（t*=tU/D， U：初期落下速さ，D：球の直径）

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

５．まとめ
二相系 LBM の最近の研究動向について紹介した．
界面を表現する 3 つの方法と圧力計算の 2 つの方法を説明した．3 章で紹介したように，密度比<100では，どの組み合わせを用いても計算可能である．
しかし，密度比≒1000 では，現在までのところ，(3)式を用いて連続の式∇・u=0 を満たす（流体力学的）圧力 p を求める手法のみが成功している．
理想的には，Free-Energy Model と(1)および(2)式
による圧力計算法を組み合わせた，熱力学的一貫性のある二相系 LBM の開発が望ましい．現在，粒子速度を増やしたり，MRT（Multiple-Relaxation-Time，
多重緩和時間）法や Entropic LBM による圧力計算の安定化が追求されているが[12]，密度比<100 が限界である．何らかの技術的ブレークスルーによる新
しいスキームの開発が望まれる．
一方，Free-Energy Model と(3)式による圧力計算
法を組み合わせた二相系 LBM は，密度比≒1000 までの計算に成功しているが，本稿では述べなかった
圧力計算以外の課題（高レイノルズ数での安定性，界面厚さの影響，界面近傍の擬似流動など）がまだ残されている．これらの課題についても，研究が進んでいる[21]．
何れにしても，二相系 LBM の最大の特徴は，圧力のポアソン方程式を解く必要がないために，アル
ゴリズムが簡単であること，ならびに並列計算に適していることである．まだまだ解決すべき課題は多いが，今後，手法の発展とともにますます適用分野が広がることを期待している．
本稿が読者の二相系 LBM の理解への一助になれ
ば幸いである．

参考文献
二相系 LBM の論文は多数発表されているので，すべてを列挙することはできない．以下では，オリジナル論文と最近の重要論文のみを挙げている．

謝辞
本稿は，第 25 回 LBM 研究会（日本機械学会流体工学部門「格子ボルツマン法の基礎と応用に関する研究会」，2016.6.3，富山大学）で発表した内容をまとめたものである．本研究会の参加者の皆様からいただいた貴重なコメントに感謝いたします．
また，本稿で紹介した筆者らの研究の一部は，科学研究費助成金基盤研究(C)(26420108)および「京」を中核とする HPCI システム利用課題(hp120112,
hp140025, hp150087, hp160082)の支援を受けて実施した．

components, Phys. Rev. E 47 (1993) 1815.
[11] Chen, L. et al., A critical review of thepseudopotential multiphase lattice Boltzmannmodel: methods and applications, Int. J. Heat Mass

[1] 稲室隆二，格子ボルツマン法による複雑流れの数値解析，ながれ 18 (1999) 31．
[2] 稲室隆二，格子ボルツマン法‐新しい流体シミュレーション法‐，物性研究 77 (2001) 197．
[3] Inamuro, T., Lattice Boltzmann methods for viscous
fluid flows and for two-phase fluid flows, Fluid Dyn.
Res. 38 (2006) 641.
[4] Inamuro, T., Lattice Boltzmann methods for moving
boundary flows, Fluid Dyn. Res. 44 (2012) 024001.
[5] Inamuro, T. et al., A lattice Boltzmann method forincompressible two-phase flows with large density
differences, J. Comput. Phys. 198 (2004) 628.
[6] Lee, T. and Lin, C.-L., A stable discretization of thelattice Boltzmann equation for simulation ofincompressible two-phase flows at high density ratio,
J. Comput. Phys. 206 (2005) 16.
[7] Zheng, H. W. et al., A lattice Boltzmann model for
multiphase flows with large density ratio, J. Comput.
Phys. 218 (2006) 353.
[8] Gunstensen, A. K. et al., Lattice Boltzmann modelof immiscible fluids, Phys. Rev. A 43 (1991) 4320.
[9] Huang, H. et al., On simulations of high-density ratio
flows using color-gradient multiphase lattice
Boltzmann models, Int. J. Mod. Phys. C 24 (2013)
1350021.
[10] Shan, X. W. and Chen, H. D., Lattice Boltzmannmodel for simulating flows with multiple phases andTransfer 76 (2014) 210.
[12] Montessori, A. et al., Three-dimensional latticepseudo-potentials for multiphase flow simulations athigh density ratios, J. Stat. Phys. 161 (2015) 1404.
[13] Swift, M. R. et al., Lattice Boltzmann simulation of
nonideal fluids, Phys. Rev. Lett. 75 (1995) 830.
[14] Mazloomi, A. et al., Entropic lattice Boltzmann
method for multiphase flows, Phys. Rev. Lett. 114

伝熱2016 年 10 月

- 12 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展(2015) 174502.
[15] Karni, S., Hybrid multifluid algorithms, SIAM J.

[19] Suzuki, K., An immersed boundary-lattice
Boltzmann method for moving boundary flows andits application to flapping flight, PhD thesis, Kyoto
University (2014).
[20] Yano, M. and Inamuro, T., Gas-liquid-solid threephase flow simulations using two-phase lattice
Boltzmann method combined with the immersedboundary method, The 24th Int. Conf. on Discrete
Simulation of Fluid Dynamics, 13-17 July 2015, The
Royal Society of Edinburgh, UK.
[21] The 25th Int. Conf. on Discrete Simulation of Fluid
Dynamics, 4-8 July 2016, SUSTech, Shenzhen,
China.

Sci. Comput. 17 (1996) 1019.
[16] Wang, Y. et al., An improved multiphase lattice
Boltzmann flux solver for three-dimensional flows
with large density ratio and high Reynolds number,J. Comput. Phys. 302 (2015) 41.
[17] Liu M. et al., A modified pseudopotential for alattice Boltzmann simulation of bubbly flow, Chem.
Eng. Sci. 65 (2010) 5615.
[18] Inamuro, T. et al., An improved lattice Boltzmann
method for incompressible two-phase flows withlarge density differences, Comput. Fluids 137 (2016)
55.

伝熱2016 年 10 月

- 13 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

格子ボルツマン法による乱流の数値解析
Turbulent Flow Simulation by the Lattice Boltzmann Method

須賀 一彦（大阪府立大学），桑田 祐丞（東京理科大学）
Kazuhiko SUGA (Osaka Prefecture University), Yusuke KUWATA (Tokyo University of Science)e-mail: suga@me.osakafu-u.ac.jp
１．はじめに乱流の直接数値解析（DNS）やラージ・エディ・
シミュレーション（LES）に用いられる数値解析
手法には，スペクトル法や有限差分法などが従来用いられてきたことは周知の事実である．スペク
トル法は，支配方程式の空間微分を波数空間で行うため，数値誤差が格子数に対して指数関数的に減少し，高い計算精度を得られるが，流路形状と境界条件に制約がある．いっぽうで，有限差分法
や有限体積法は，空間微分を物理空間で離散化近似するため，形状や境界条件の複雑化に対して柔軟性が高く，工学的な乱流計算手法の主流となっている．しかしながら，大規模な実用解析を実行するにあたっては，格子生成や並列計算，非圧縮
性流れにおける圧力のポワソン方程式の計算コストなどに関して，いくつかの容易でない課題がある．
そこで近年，格子生成の必要が無く，複雑流路表現能力が高く，プログラムが簡単で並列計算に適し，かつ圧力のポワソン方程式を解かなくてもよい格子ボルツマン法（LBM）が脚光を浴びている．いっぽうで，LBM による乱流解析は歴史が浅
いこともあり，課題も少なくない．したがって，
筆者らはそれらの課題に取り組み，ここ数年でかなり手法の信頼性を向上させてきたと考えている．
そこで本稿では，筆者らが取り組んできた LBM
による非圧縮性流体の乱流シミュレーション法に関して紹介する．

2.1

２．格子ボルツマン法
格子 BGK(LBGK)モデル

2016 年 10 月と表したものを格子 BGK(LBGK)モデルと呼ぶ．
ここで τ は単一の値をとることから単一緩和時間（SRT）LBM とも呼ばれる．このように LBM では，速度場を解くにあたって式(1)で表される速度分布関数の時間発展方程式を陽的に解くだけなので，プログラムは極めて単純であり，かつ並列化も容易である．また，精度は空間・時間ともに二次精度[2]であり，曲面物体形状も補間法により二次精度で表現可能である[3]．局所平衡分布関数
f eq は気体定数 R ，密度  , 流体速度 u および温度
T を用いて， ξ  u 1  (ξ  u) 2 u 2   f eq  x, t   w   0  α   α 2   ,
 RT
RT   
2  ( RT )(2)と表される．ここで 0 は基準値で単位圧力に設定
し，  0   とする．非熱流体 LBM では音速と温度は cs  RT の関係にある．外力項 F は加速
度を a として，  ξ  a  ξ  u  a  u  F  x, t   w    α  1  α   ,RT  RT    RT 

格子ボルツマン方程式は，ボルツマン方程式の速度空間を Q 個の離散速度 ξ α {  0,1, , Q  1} に
よって離散化する. 離散速度モデルには多くの提案があるが，非圧縮（非熱流体）モデルではいずれのモデルでも無次元音速 cs / c は固定値として扱う．ここで, c   / t は格子間隔と時間刻みの伝熱比で，通常は単位速度に設定される．離散速度 ξ αと重み係数 w は用いる離散速度モデルによって異なる．
（表 1，2 参照）
格子ボルツマン方程式では，BGK 近似[1]を用いて，デカルト座標 x 上での速度分布関数
f (x, t ) の発展方程式を1
f   x  ξ α t , t  t   f   x, t    f   x, t   f eq  x, t  

1 

(1)
  1   t F  x, t  , 2 

(3)とする．流体の密度，速度および圧力は，
Q 1

Q 1

 0

 0

   f  , 0 u   f  ξ    0

ta, p  cs2 ,2(4)と与えられる．また，緩和時間 τ と粘性係数の間- 14 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展には   (   1 / 2) 0 cs2 t の関係がある．

年になってこの現象が指摘されるようになったが[4-6]，そのメカニズムについては，筆者ら[7]が明らかにするまでは不明であった．格子ボルツマン方程式の誤差解析[2]を離散速度モデルごとに行い，運動量の式

2.2 離散速度モデル
三次元直交立方格子上の離散速度モデルには，
D3Q7, D3Q15, D3Q19, D3Q27, D3Q39 など多く提案されている．しかし，一般には隣接するキュー
ブ上の格子点のみを参照する D3Q7～D3Q27 モデ
ルが用いられる．D3Q7 モデルは差分法と同じように，着目する格子点と直接格子線で繋がった前後上下左右にある 6 参照点の，合計 7 点に向かう
離散速度で構成される．D3Q15, D3Q19 モデルでは，図 1 に示すように斜めの格子点も参照し，
D3Q27 モデルでは全ての隣接点に向かう速度を
扱う．

ui
u
p 1  2 ui uj i  E i ,
t
x j
xi Re x 2jに現れる誤差項Ei を円筒座標系で評価すると，例
えば円周方向の誤差成分は，E  D 3Q15   2



1  1 2
5 
4

 
 2
 3  u u   O( 4 ),2
4  r rr r r 


(a)

(6)

D3Q15

εr+ S

εθ+ S

-0.02
-0.02
-0.06
0.02

0.2

(b) D3Q19model3( 62  6  1) 
3 2
1  
sin 4   3 6r r 2 4  r

3   2
1  3 4 28 
4
 2  u z  uz uz    3 
 2
 3  urur2
4  rr r r r
r r r 

３．LBM による乱流解析
3.1 D3Q27 離散速度モデル
三 次 元 流 れ の 解 析 で は 通 常 D3Q15, D3Q19,
D3Q27 のいずれかのモデルが用いられるが，どれでも結果に差がなければ，計算負荷の小さい方が良いはずである．たとえば，図 2 に示すように層流の計算では，離散速度モデルによって結果に差異はない．しかし，流れが乱流に変化すると管内速度場が離散速度モデルによって変化し，D3Q27
モデルのみが妥当な結果を示すことがわかる．近(a) D3Q15 model(5)1(c) D3Q27 model

(b)

10

y+

100

100

10

1

y+

D3Q19

εr+ S
0.02
0.06

0.2

εθ+ S
0.02
-0.02

図1

離散速度モデル
0.2

D3Q15

D3Q19

1

D3Q27

(c)

10

y+

100

10

1

y+

D3Q27

εr+ S

Laminar

100

0.2

εθ+ S
-4

0.001

1.0×10
1.0×10-4-1.0×10-1.0×10

Turbulent

0.0001
0.2

2016 年 10 月

10

y+

100

100

10

y+

1

0.2

図 3 離散速度モデルによる誤差分布; S=Δ+2(-6τ2+6τ-1)

図 2 離散速度モデルの違いによる円管内軸方
向平均流速分布の変化

伝熱1- 15 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展のようにマトリクス表記する．ここで記号 f はディラックのケット・ベクトルであり， f
 ( f 0 , f1 ,… , fQ 1 )T である．また， M は Q  Q のマトリクスであり， m = M  f により分布関数 f をモーメント m に変換する．衝突マトリクス S は対角行列 S  diag( s0 , s1 ,… , sQ 1 ) であり，緩和時間の逆数に対応する．ここで e  ξ α / c と表記すると，
D3Q27 モデルのモーメント m には，密度  ，流束
j ，運動エネルギー q :となる．層流で放物形の平均速度分布を持つ場合，
式(6)-(8)における O(Δ2)の誤差項は全てゼロになり，D3Q15, D3Q19, D3Q27 のいずれのモデルでもO(Δ4)以上の精度で速度場が計算できる．しかし，乱流になると一変し，レイノルズ応力と平均速度の勾配に応じて，円周方向に一定の周期で O(Δ2)の誤差応力が生じるようになる．O(Δ2)といえども，
本来は応力が全く生じない領域に周期的に応力が現れるので，非物理的な現象が発生すると考えられる．この誤差項の分布を Re=5400 の円管内乱流
DNS データを用いて可視化したものが図 3 である．
全てのモデルで周期的な誤差は存在するが，
D3Q15 と D3Q19 では D3Q27 に比して二桁程度絶対値が大きくなることがわかる．そのため，図 2に示されるように顕著な予測差になって現れる．
また軸対象流れに限らず，矩形ダクト流れであっても，同様の理由で D3Q15 と D3Q19 では，管軸周りに 45 度回転させた格子では大きく結果が変わり[5]，乱流解析では少なくとも D3Q27 モデル
を用いなければならないことがわかる．
3.2 D3Q27 多緩和時間モデル
式(1)の LBGK モデルで高 Re 数流れを解析する
際，解が振動したりして計算が不安定になることが知られている[8]．これは，LBGK モデルでは速度場の様々な物理量の変化に対し，全て単一の緩
和時間 τ で記述しているためである考えられており，これを物理モーメントそれぞれに独立に与えようというのが多緩和時間(MRT)モデル[9]である．
MRT モデルでは分布関数の時間発展を

m0    f   , m13    ex13 f   jx13 ,等 27 成分あり，7 次のモーメント項まで現れる．
前述のように，乱流解析のためには少なくとも離
(a)

(b)

(c)

(d)
u2
w2
v2

2016 年 10 月

Pk

+

d kp

+

d k

+



20

Sv
Sw

Spectrum (Kim et al.)Spectrum (Horiuti)(Kasagi et al.)
LBM

Fv

15
10

Fw

5

Su
0.4

0.6

0.8

Fu

0
0.0

1.0

1

0.2

0.4

10

E11

10

-1

+

-1

10

E11

-2

+

10

E33

-2

10

-3

10

0.8

1.0

0

E33
E22

+
+

-3

10

-4

+

E22

-4

10

10

-5

10

-5

10

-6

10

-7

10

-6

10

-7

10

-8

-8

10

10

Kim et al.
D3Q27 MRT

-9

10

-3

10

-2

10

+

kx

図4

0.6

10

+

0

10

y/

1

10

-10

- 16 -

25

Spectrum (Kim et al.)Spectrum (Kasagi
et al.)(Horiuti)
LBM

0.2

d kt

y/

f ( x  ξ  t , t  t )  f (x, t )

伝熱

1.2
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
0.0

10

(9)

S
 M 1 S  m(x, t )  m eq (x, t )   M 1  I   M F t ,




2


(10)

m4    (e2x1  e2x2  e2x3 ) f   q,

-uv+


4
1  1 25 
4
 3  ur ur  
 2
 3  u u   O (  4 ), (8)24  r rr r r r 

Flatness3( 62  6  1) 
4 28 
1  
sin 4   3 
 2

2
6
4
r
r
r
r


r




+

E  D 3Q 27   2

(7)

Eii


1  1 2
5 
4
 
 2
 3  u u   O(  4 ),2
4  r rr r r Skewness

3   2
1  3 4  28 
4 2  u z  uz uz    3   urur
4  rr r 
r r 2 r 2 r r 3 +3( 62  6  1) 
3 2
 1 
sin 4    3 6r r 2 4  rEii

E  D 3Q19   2

Kim et al.
D3Q27 MRT

-9

10

-10

-1

10

10

-3

10

-2

-1

10

10
+

kz

チャネル乱流(Reτ=180)の乱流統計量

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展
n   c /  f  t c / t f を用い，散速度モデルは D3Q27 モデルが望ましいので，筆者らによって MRT 化が行われ，D3Q27 MRT モデ
ル[10]が提案されている．そこでは変換マトリク
ス M をグラム・シュミットの正規直交化法により与え，そのときの衝突マトリクス S は，動粘性係数 と直結している成分 s5 を除いては，計算の安定性等を加味しながら数値実験によって最適化されている．
なお，LES を行う場合は，
11
  δt , s5 2    sgs  cs2 f eq ,c  f eq , f , f neq ,c  n

f neq , f ,(12)
fの関係を使い情報を交換する．ここで肩文字 c,fはそれぞれ粗，密を意味し，neq は非平衡分布関数 ( f  feq ) を意味する．粗格子から密格子に情報を伝える際には，時間および空間において補間法を併用する．
LBM における LGR 法の問題点は，乱流解析を行うと，レイノルズ応力などの乱れの二次モーメ
ント以上の高次相関量において，図 6(a)(b)に示す
ように，破線で示す格子の継ぎ目で分布に捩れが出てしまうことである[15,16]．特に，格子間隔を比較的広くとれる LES を行うと顕著に現れ，図6(b)のように，有限体積法（FVM）による LES と比べても差は顕著で，LBM では格子間隔の変化に敏感であることがわかる．
従来の LGR 法におけるこの問題は，格子から格子に分布関数を受け渡す際，質量や運動量の保存(11)のように動粘性係数と SGS 渦粘性係数を足し合わせて s5 に関連付ける．また，格子ボルツマン方程式で扱う圧力は，式(4)からわかるように音速や密度と密接に関連した静圧であるので，ナビエ・ストークス ソルバー中の乱流圧力（静圧＋乱流
エネルギー）のように取り扱うことは原理的にできない．そこで，通常のフィルター操作に加え，テストフィルター
を用いて，
3
2
k sgs  Ckes  k 1 ( uk  uk ) のように乱流エネル
ギーを見積もり[11]，外力項 F に 2 / 3ksgs を反映させる．係数は通常 Ckes  1.0 とする．
図 4 は，チャネル乱流を D3Q27 MRT LBM で
DNS 解析した結果を，スペクトル法による結果[12,13]と比較したものである[10]．計算領域 4πδ
×2δ×2πδ を等間隔格子 1539×240×754 で分割し，解像度は全領域で Δ+=1.5，平均マッハ数は Ub/cs=0.01としている．平均速度分布，レイノルズ応力の各成分，乱流エネルギー輸送方程式の各収支項，
skewness，flatness，流れ方向とスパン方向の一次
元エネルギースペクトルの分布のどれを見ても，LBM の結果はスペクトル法の結果をほぼ完全に再現している．このことから，D3Q27 MRT LBM
による乱流の解析は十分に信頼できることがわか
る．
3.3

(a)

∆c

coarse grid

,
A

f→c

B fβ c

fβff

,
B

fβc→f
A
f
fine grid(b)coarse grid

c A

B

c
fβcf→c P

fβcc

f

IBC 局所細密格子法

2016 年 10 月

fβcc

∆f

LBM では規則格子を用いるので，境界層の解像度を上げようとすると，全領域で膨大な数の計算
格子点が必要になる．そこで大規模解析には局所細密格子（local grid refinement: LGR）法[14]を用いることが必須条件となる．LGR 法では図 5(a)に示すようにオーバーラップ部を設け，格子間隔比

伝熱

c

fβcf

A

c→f
P f fβcBfine grid

図 5 粗密格子の接続: (a) 標準的手法[14], (b)IBC 局所細密格子法[17]- 17 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展
条件  fneq  0,  ξ  fneq  0 が厳密には満たされないためであることを筆者らは明らかにし，それを修正する以下に概要を示す，imbalance-correction（IBC）LGR 法[17]を提案している．
図 5(b)の粗格子上の点Ｐc において密格子から与えられる左方向 c に向かう分布関数の成分の和とそれ以外の成分 c の和に分けて質量と運動量の残差を計算すると，
 

f

neq ,c

c

c

u 

f
ˆ c

ξ f

neq ,cc ccneq , f  c,
ˆ c
neq , f  c,
ˆ c ˆ cu*  u  u u ,(16)(a) 3.02s c 2  
 1
    ux  1  1  r   1  2   , cs  
 6
1

(a)

4.0
Reynolds stresses





1.0
vv

(21)

のように補正係数を見積もることができる．なお，
これを数回繰り返せば完全に収束させることも可能であるが，筆者らの経験上 1 回で十分である．
図 7 は，IBC LGR 法と従来法である Dupuis &Chopard(DC)法[14]によるチャネル乱流の LES 結果を比較したものである．SGS モデルは WALE モ
デ ル [18] を 用 い ， 密 格 子 ， 粗 格 子 の 大 き さ は

Kim et al.
Dupuis and Chopard method
IBC method

3.5
uu

(20)

 c 2  
 1
 uy   uz  1  1  r   2   , 18
 cs  

DNS
LBM(LES)

2.5



3.0
2.5
2.0
1.5
1.0
0.5

0.0
0

100

y+

200

300

0.0

400

(b)

3.0

DNS

2.5

FVM 
LES
LBM 

uu



2.0
1.5

ww



1.0
0.5

vv



0

100

200

300

-2
-4

Loss: 5%-6
Dupuisand ChopardDCmethodImproved
methodIBC
method

0

400

y+

5

1x10

5

2x10
Iteration

5

3x10

5

4x10

図 7 IBC LGR 法によるチャネル乱流 LES(Reτ=180) [17]: (a) レイノルズ応力, (b) 質量保存性

図 6 LGR 法による LES(Reτ=395): (a) オープン
チャネル乱流 [15], (b) チャネル乱流[16]2016 年 10 月

Loss: 0.2%

0

-8

0.0

20 40 60 80 100 120 140 160 180+y2
Total mass error [%](b) 3.5伝熱

2sが得られる．ここで式 (17)-(19) に現れる高次項
O (u 2 / cs2 , u u / cs2 ) を無視し，左辺をゼロとした連立方程式を解くと，と係数   ,  u を導入して補正し，再度残差を計算すると， r  nc /  f として，1 r
*  (1    )   

6
(17)
1  r cux
222
  ux 0

O
(

/
c
,

/
c
),
u
s
s
u
u
6 cs2

0.5

(19)

 O(u / c , u u / c ),
2

1

(15)

ww

1  r c u y , z
18 cs22(14)

*       ,

1.5

(18)

1  r c u x
 O(u 2 / cs2 , u u / cs2 ),6 cs2

*u y ,z  (1   u y ,z )u y ,z  u y ,zとなり， , u がゼロになるとき保存則が満たされる．そこで，これら残差を使って密度，速度を

2.0

1 r
c / 062

  ux

(13)

ξ f
ˆ c

*ux  (1   ux )ux  

- 18 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展 f   3.8 ， c  7.6 であり，緩和層内の y   17.7の位置で接合している．なお，計算領域，マッハ数などは図 4 の DNS と同様である．図 7(a)に明らかなように IBC LGR 法は捩れた点も無く，ほぼ完
全に DNS[12]の結果を再現するが，DC 法では破線の格子継ぎ目で分布が捩れるだけではなく，緩和層内に継ぎ目があることで大きく DNS から外れた結果となっている．また，図 7(b)は系全体の質量保存性を確認したものである．DC 法で 40 万
ステップ計算すると約 5%の誤差が生じるが，IBC
LGR 法ではこれを 0.2%まで改善できる．
3.4

いっぽうで，乱流温度場を解析する際，境界形状が複雑化すると，図 9(a)に示すように層流解析では問題なくても，通常の時間発展方程式(22)のままでは図 9(b)に示されるように解が不安定化し，いずれは発散してしまうことが報告されている[20]．これを解決するためには，regularization（正規化）が有効である．エネルギー密度分布関数の正規化した時間発展方程式は，(a) 1.0

+
+

0.4

0.0
0.0

(c)

0.2

0.4 0.6
y/H

Nu

0.8

1.0

0
0

50

100
+
y

150

50
A.K. Dhiman et al. 2005
D3Q19 with reg.
D3Q7 MRT

40
30
20

Pr=50
Pr=10

10 Pr=1

1[ g  ( x, t )  g eq ( x, t )],
g

0
A

B

C

D

図 8 離散速度モデルによる温度場の解析結果[19]: (a) チャネル乱流 LES(Reτ=180)の平均温度分布, (b) 同乱流熱流束, (c) 角柱周りの熱伝達(Re=45)(22)を速度場と合わせて解く．強制対流スカラー場のモデル[19]では，温度 T ，平衡分布関数 g eq は，Q 1 ξ u T   g  , g αeq  wT 1  α 2  ,
cs 
 0


4
2

0.2

温度場を解析する際，非圧縮性流体の LBM では，非圧縮性ナビエ・ストークス ソルバーと同様に，独立したエネルギー式をブシネスク近似で
速度場と連結する二分布関数法を用い，エネルギ
ー密度分布関数 g  の LBGK 表記の発展方程式
g  ( x  ξ α t , t  t )  g  ( x, t ) 

0.6

Kawamura et al. 2000
D3Q19 with reg.
D3Q7 MRT

6
u' '

T/

乱流温度場

(b) 8

Kawamura et al. 2000
D3Q19 with reg.
D3Q7 MRT

0.8(23)(a) Re=30と表す．また，緩和時間  g は温度拡散係数  と1    g   cs2 t ,
(24)
2

の関係があり，LES を行う際には    0   sgs / Prtのように分子温度拡散係数 0 に SGS 温度拡散係数を加えて緩和時間と対応付ける．
スカラー場の離散速度モデルに関しては，速度場のように D3Q27 モデルでないと不具合がおきることは報告されていない．図 8(a)(b)に示される
ように，D3Q7 モデルと D3Q19 モデルのいずれを用いてもチャネル乱流の解析結果には差がなく[20]，DNS データ[21]をよく再現する．しかし，角柱周りの熱伝達[22]を解析すると，層流であっても D3Q7 モデルでは図 8(c)に明らかなように角部の Nu 数に捩れが生じていることが分かり，離散速度成分が最少の D3Q7 モデルでは，複雑形状周りの解析には不適切であることがわかる[20]．

伝熱2016 年 10 月(b) Re=3000
T/ΔT
1

0

(c) Re=3000 with Reg.

図 9 正規化の効果: (a)(b) 正規化を用いない場合の角柱群周りの瞬時温度分布，(c)正規化を用いた
場合の瞬時温度分布- 19 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展
1 g   x  ξ α t , t  t   g eq  x, t   1   g   x, t  , (25)
  
g 

と書け，正規化した非平衡分布関数は，g  w

 i

Q 1

g  ,

cs2 0

neq
i


(26)

と計算される．これらを解くことで計算は安定化
し，図 9(c)のように問題は回避できる．
４．おわりに
以上まとめると，乱流速度場の解析には D3Q27
MRT LBM[10]が適しており，局所細密格子法にはIBC LGR 法[17]が有効である．また，二分布関数法を用いた乱流スカラー場の解析には速度場とは別 に ス カ ラ ー 場 に 正 規 化 LBGK モ デ ル [20] をD3Q19 モデル（または，筆者らは確認していない
が，D3Q15 モデル）とともに用いるのが有効である．なお，誌面の都合で離散速度モデルのパラメ
ーターについて本文中に詳細を記述しなかったが，

D3Q19, D3Q27 モデルについて，表 1, 2 に示す．
また，最初に述べたように LBM は，プログラ
ムが簡単で並列化が容易であることに加え，格子生成が不要で，複雑流路の形状表現に優れている．
したがって，図 10 に一例を示すように，多孔体材
料のような既存手法の適用が困難な不規則形状を有する複雑構造体周りの流れ場[23, 24]に対して，その威力を遺憾なく発揮できると考えられる．
本稿では，LBM による乱流の解析法についての筆者らの取り組みを紹介した．拙稿によって，適切に実行すれば，LBM が高い精度で乱流や乱流伝熱を解析できる有望な手法であることが，読者にご理解いただければ幸甚である．
表 1 D3Q19 離散速度モデルのパラメーターcsξα / cc

1
3(0, 0, 0)(1, 0, 0), (0, 1, 0), (0, 0, 1)(1, 1, 0), (1, 0, 1), (0, 1, 1)

表 2 D3Q27 離散速度モデルのパラメーターξα / c

(a)

cs
c

1
3

0
1~ 6
7 ~ 1819 ~ 26(0, 0, 0)(1, 0, 0), (0, 1, 0), (0, 0, 1)(1, 1, 0), (1, 0, 1), (0, 1, 1)(1, 1, 1)

1/ 3
1 / 181 / 36w
8 / 272 / 27
1 / 541 / 216

参考文献[1] Bhatnagar, P. L., Gross, E. P. and Krook, M., A
model for collision processes in gases. I. Small
amplitude processes in charged and neutral(b)

図 10 空隙率 0.8 のカーボンペーパー界面乱流のLES 解析例: (a) ひずみ速度テンソル第 2 不変量，(b) カーボンペーパー内速度コンター；Reb=3000

伝熱

0
1~ 6
7 ~ 18w

2016 年 10 月

- 20 -

one-component systems, Phys. Rev., 94 (1954),
511-525．
[2] Holdych, D.J., Noble, D.R., Georgiadis, J.G. andBuckius, R.O., Truncation error analysis of lattice
Boltzmann methods, J. Comput. Phys., 193 (2004),
595 – 619.
[3] Chun, B. and Ladd, A.J., Interpolated boundary
condition for lattice Boltzmann simulations offlows in narrow gaps, Phys. Rev. E, 75 (2007)
066705.
[4] Mayer, G. and Hazi, G., Direct numerical andlarge eddy simulation of longitudinal flow along
triangular array of rods using the latticeJ. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展
Dynamic subgrid scale modeling of turbulent
flows using lattice-Boltzmann method, Physica A,

Boltzmann method, Math. Comput. Simul., 72(2006), 173 – 178.
[5] Kang, S.K. and Hassan, Y.A., The effect of latticemodels within the lattice Boltzmann method in thesimulation of wall-bounded turbulent flows, J.

388 (2009), 2640–2658.
[16] Touil, H., Ricot, D. and Leveque, E., Direct andlarge-eddy simulation of turbulent flows on
composite multi-resolution grids by the latticeComput. Phys., 232 (2013), 100 – 117.
[6] Geller, S., Uphoff, S. and Krafczyk, M., Turbulent
jet computations based on MRT and cascaded
lattice boltzmann models, Comput. Math. Appl.,

Boltzmann method, J. Comput. Phys., 256 (2014)
220–233.
[17] Kuwata, Y. and Suga, K., Imbalance-correction
grid-refinement method for lattice Boltzmann

65 (2013), 1956 – 1966.
[7] Kuwata, Y. and Suga, K., Anomaly of the lattice
Boltzmannmethods
inthree-dimensional

flow simulations, J. Comput. Phys., 311 (2016),
348 –362.
[18] Nicoud, F. and Ducros, F., Subgrid-scale stressmodelling based on the square of the velocity

cylindrical flows, J. Comput. Phys., 280 (2015),
563 –569.
[8] Yu, H., Luo, L.-S. and Girimaji, S.S., LES of
turbulent square jet flow using an MRT lattice

gradient tensor, Flow Turbulence Combust., 62(1999) 183–200.
[19] Yoshino, M. and Inamuro, T., Lattice Boltzmann
simulations for flow and heat/mass transfer
problems in a three-dimensional porous structure,

Boltzmann model, Comput. Fluids, 35 (2006), 957– 965.
[9] d’Humi´eres, D., Ginzburg, I., Krafczyk, M.,
Lallemand, P. and Luo, L.-S., Multiplerelaxation- time lattice Boltzmann models in threeInt. J. Num. Meth. Fluids, 43 (2003) 183–198.
[20] Sakurai. Y., Kuwata, Y. and Suga, K., Lattice
Boltzmann LES of conjugate turbulent heattransfer in square rod arrays, Proc. ETMM10,
Marbella, Spain (2014) S20.4.
[21] Kawamura, H., Abe, H., and Shingai, H., DNS of
turbulence and heat transport in a channel flow
with different Reynolds and Prandtl numbers andboundary conditions, Proc. THMT3, Nagoya,Japan (2000) 15–32.
[22] Dhiman, A.K., Chhabra, R.P., and Eswaran, V.,
Flow and heat transfer across a confined squarecylinder in the steady flow regime: Effect ofdimensions, Phil. Trans. R. Soc. A, 360 (2002),
437–451.
[10] Suga, K., Kuwata, Y., Takashima, K. andChikasue, R., A D3Q27 multiple relaxation- timelattice Boltzmann method for turbulent flows,
Comput. Math. Appl., 69 (2015), 518–529.
[11] Inagaki, M., Nagaoka, M., Horinouchi, N. andSuga, K., Large eddy simulation analysis ofenginesteadyintake
flows
using
a
mixed-time-scale subgrid scale model, Int. J.
Engine Res. , 11 (2010) 229–241.
[12] Kim, J., Moin, P. and Moser, R., Turbulence
statistics in fully developed channel flow at low

Peclet number, Int. J. Heat Mass Transfer, 48(2005) 4598－4614.
[23] Kuwata, Y. and Suga, K., Large eddy simulationsof pore-scale turbulent flows in porous media by
the lattice Boltzmann method, Int. J. Heat Fluid

Reynolds number, J. Fluid Mech., 177 (1987),
133–166.
[13] Kasagi, N., Horiuti, K., Miyake, Y., Miyauchi, T.
and Nagano, Y., DNS database of Turbulence andHeat Transfer, Univ. Tokyo, Japan (1993).
[14] Dupuis, A. and Chopard, B., Theory andapplications of an alternative lattice Boltzmann

Flow, 55 (2015), 143 –157.
[24] Kuwata, Y. and Suga, K., Lattice Boltzmann
direct numerical simulation of interface
turbulence over porous and rough walls, Int. J.
Heat Fluid Flow, (2016), DOI:10.1016/
j.ijheatfluidflow. 2016.03.006.

grid refinement algorithm, Phys. Rev. E, 67(2003) 066707.
[15] Premnath, K.N., Pattison, M.J. and Banerjee, S.,

伝熱2016 年 10 月

- 21 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

GPU スパコンによる大規模 LBM 計算
Large-scale Computations of Lattice Boltzmann
Method on a GPU Supercomputer

青木 尊之（東京工業大学）
Takayuki AOKI (Tokyo Institute of Technology)
e-mail: taoki@gsic.titech.ac.jp
1. はじめにGPU（Graphics Processing Unit）は，パソコンのグ
ラフィクス・カードに搭載されて使われてきたが，最近では CPU に内蔵されたり，スマートフォンやタブレットにも搭載されている．GPU は画像を表示するために開発され，その後，より高速で，美しく，高精細な画像表示の要求により，描画機能が飛躍的に向上してきた．そして，画像表示だけではな
く汎用計算に使う「GPU コンピューティング」（GPGPU(General-Purpose computing on GPUs)と同義）の取組みが 2000 年頃から始まった．GPU はグ
ラフィクス処理を目的としているため，汎用計算に対する機能が限定されている反面，同世代の CPUに対して 10 倍以上高い演算処理能力を持つため，魅力的な演算デバイスとして次第に認識されるよ
うになっていった．2006 年には NVIDIA 社が自社の GPU に対し，GPU コンピューティングの統合開
発 環 境 と し て CUDA(Compute Unified Device
Architecture)をリリースした[1]．それまでは Cg (Cfor Graphics) 言 語 や HLSL(High Level Shader
Language) により画像処理の機能を汎用計算に置き換える必要があったが，CUDA の登場により通常の C 言語の拡張としてプログラミングすることができるようになり，GPU を汎用計算のアクセラレータとして利用する GPU コンピューティングが一気に広まり始めた．また，2009 年には CUDA のFORTRAN 版もリリースされ，GPU の性能向上に合わせた頻繁なバージョン・アップがなされている．
2012 年に NVIDIA 社，PGI 社，CAPS 社, Cray 社が協定し，CPU で実行する通常のソースコードに指示行を挿入するだけで GPU の機能を使えるようにする OpenACC[2]の開発が始まった．ホスト（CPU）
用コード，さまざまなアクセラレータ用コードを単一コードとして記述でき，メンテナンスが容易になるなどの利点がある．
並列化したいループに kernels,
loop ディレクティブを追加することにより，指定

伝熱2016 年 10 月

- 22 -された領域がアクセラレータで実行されるカーネ
ルへ変換され，GPU で実行される．CUDA でプロ
グラミングした場合と比較すると，まだ十分にGPU の本来の性能を引き出すことは難しいが，さまざまな不具合が解消されたバージョン 2.0 のOpenACC 対応のコンパイラが有償で提供されている．
スパコンへの GPU の導入は，東京工業大学・学
術国際情報センターが 2008 年に 680 個という当時としては大量の GPU を世界で初めてスパコン(TSUBAME 1.2)に導入した．2010 年には科学技術計算に適した倍精度浮動小数点演算や ECC メモリに対応する GPU がリリースされ，同年 11 月のスパ
コン Top500 のランキング（Linpack 計算のスコア）[3]では，
1 位と 3 位に中国の GPU スパコンが入り，4 位に NVIDIA Tesla M2050 を 4,224 個搭載した東京工業大学・学術国際情報センターのスパコン
TSUBAME 2.0（総合演算性能 2.4PFlops）がランクインし，GPU をアクセラレータとして搭載するス
パコンが上位を独占した．2016 年 6 月の Top10 の中では，GPU を利用するスパコンが 2 システムと減少しているが，これは次期の Pascal コア GPU のリリースが遅れているためであり，Top500 の中では GPU を搭載するスパコンが 86 システムと増え，GPU の利用が着実に普及している．
２．LBM 計算の GPU 実装
2.1 単一 GPU での実装
均一格子に対する LBM 計算を単一 GPU に対して CUDA を使って実装することは難しくない．GPUのアーキテクチャは汎用 CPU と大きく異なり，
NVIDIA の Tesla K20X を例に説明すると，192 個の演算ユニット（CUDA コア）が 1 つのストリーミン
グ・マルチプロセッサ（SMX と呼ばれている）を構成している．Tesla K20X には 14 個の SMX があり，合計で 2688 個の CUDA コアがある．1 つのJ. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展
SMX の中には 64k 個の 32bit レジスタファイルがあり，非常に高速な L1 キャッシュと共有メモリがある．1 つの SMX 内の共有メモリは 192 個の CUDA
コアからアクセスできるが，他の SMX 内の CUDA
コアからはアクセスできない．K20X を搭載したグ
ラフィクス・ボード上には 6GB のビデオメモリ（GDDR5）があり，GPU からはビデオメモリへ高速なアクセスが可能である．ビデオメモリには2688 個の CUDA コアから共有でメモリアクセスが可能である．CPU からビデオメモリには比較的低速な PCI Express バスを介したデータ通信が行われる．GPU スパコンの各ノードには，2 個程度の CPUと GPU を搭載したビデオカードが 1～4 枚程度搭載されていることが多く，さまざまなメモリが階層的構造になっていて，それぞれデータ転送の帯域が大きく異なるため，それらを考慮した実装を行うと高い実行性能が得られる．
GPU で計算する場合，1 コアが 1 スレッド以上を処理するが，通常の CPU と異なり 1 スレッドの処理は非常に軽量であり，2688 コアに対して 100 万スレッドを投入することも珍しくない．従来のプロ
グラミングでは，ループにより全格子点を順番に計算していたが，CUDA ではループの回数に相当する
スレッド数を生成し，GPU に投入する．1 スレッド
は LBM の 1～32（程度）の格子点を計算し，スレッドを投入する際に割り当てられる複数のスレッド ID（ビルトイン変数）を組み合わせ，格子点番号に対応させる（図 1）．GPU 計算は 32 個単位のInstruction で SIMT（Single Instruction Multiple Thread）
で実行されるので，if 文などの条件分岐をできるだ
け少なくするようなプログラミングを行うと高い実行性能が得られる．それ以外の部分は従来のプログラムとほぼ同じであり，LBM 計算の CUDA 実装
は比較的容易である．

K40 では 12GByte）であり，大規模計算を行うには複数の GPU を使った計算が必要となる．LBM 計算
では均一直交格子を前提とするので，通常は領域分割法が用いられ，分割された各領域の計算に 1 個のGPU が割り当てられる．分割された領域の境界近傍格子（袖領域）での計算は隣接領域の格子点にアクセスする必要があり，袖領域に対するデータ通信が発生する．グラフィクス・ボード上のデバイスメ
モリ上にあるデータの GPU 間通信は帯域の狭い
PCI Express バスを介して CPU 側のメモリに転送する必要があり，メモリアクセスと比較するとかな
り時間がかかる．さらに CPU 側のメモリから隣接の GPU に転送する必要があるため，大規模計算において GPU 間のデータ通信は大きなオーバーヘッ
ドになる．最新の NVIDIA の GPU Direct RDMA を利用することで異なるノードの GPU との通信の手
続きが簡単になり，通信速度も 2 倍近く速くなる
が，依然として大きなオーバーヘッドがあることには変わりがない．これから登場する Pascal アーキテ
クチャの Tesla P100 では NV-Link が導入され，各ノードのマザーボード上で NV-Link で連結された
GPU 間では非常に高速なデータ転送が可能になる．
しかし，依然として複数ノードをつかう大規模計算に対しての状況は変わらない．高い実行性能を得るためには，分割された各領域において境界格子を先に計算し，その終了とともに GPU 間データ通信を開始する．それと同時に内側の格子の計算を並行して実行する「計算と通信のオーバーラップ（図 2）」
により，GPU 間データ通信の時間を可能な限り隠蔽する必要がある．

図 2 計算と通信のオーバーラップ
３．1m 格子を用いた都市部 10km 四方の大規模 LES
気流シミュレーション LBM 計算の GPU 実装
東京都心部 10km 四方のエリアに対し，実際の建造物のデータを元に 1m 間隔の格子解像度で詳細な気流シミュレーションを行う．都市部の建造物は高図 1 ビルトイン変数と配列要素の対応
2.2 複数 GPU の利用
単一のグラフィクス・カードに搭載されるビデオ
メモリのサイズはせいぜい数 GByte（最新の Tesla

伝熱2016 年 10 月

- 23 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展のスパコン TSUBAME2.0 のほぼ全ノードの 4,032
個の GPU を用い，計算格子 10,080×10,240×512 に対して流入速度は対数則による速度分布関数で与え，高さ 100m において風速を 10m/s と設定して計算を行った．流入速度で規格化したレイノルズ数はおよそ 106 程度となる．
図 3 は多数の粒子を風速に乗せて移流させたときのスナップショットであり，動画にすると風速分布を把握することができる．高層ビルにより, 速度の速い風が上空へと押し出されている様子が確認できる. さらに, 中央部の高層ビル後方部では上空の風の巻き込みにより, 下向きに強い風が吹いている事などが確認された．また，幹線道路に沿って流れる「風の道」も再現されていることが分かった[5]．将来的には，台風の際の被害などが飛躍的な精度で予測できるようになったり，排ガス，事故やテロによる有毒ガスなどの汚染物質の拡散も詳細に予測可能となると期待できる．

層ビルが密集した複雑な形状をしており，それらに対する気流（風）は乱流状態となるため，汚染物質などの拡散過程に大きな影響を与える．これを精度良く予測するために，1m 間隔格子を用いて広範囲の大規模計算を行う．
3.1 ラージエディ・シミュレーション
都市部の高いレイノルズ数の計算には，ナヴィ
エ・ストークス方程式を解いて計算する場合と同じ
く，乱流モデルの導入が不可欠である．ラージエデ
ィ・シミュレーション（LES: Large-eddy simulation）
は格子解像度以下の変動成分を渦粘性で表し，高いレイノルズ数の流れに対しても安定かつ高精度に計算できる．しかし，乱流モデルとして広く用いられている動的スマゴリンスキー・モデルでは，格子で解像できる流れ場からモデル係数を動的に求めるため，テストフィルターとして広範囲の空間平均（メモリアクセス）が必要となり，大規模計算には適さないことが分かる．
コヒーレント構造スマゴリンスキー・モデル[4]は，乱流のコヒーレント構造を抽出するために，速度勾配テンソルの第 2 不変量と速度勾配テンソルからコヒーレント構造関数により局所的にモデル係数を決定可能な手法である．このモデル係数は明確な上限を持ち，その分散も小さい．また，層流ではコヒーレント構造が無くなるので自動的にモデ
ル定数はゼロとなり，壁面に対してもモデルが正しい振る舞いを示す．
LBM においても差分法と同様に渦粘性に基づ
く SGS モデルを用いることが可能であり，SGS モ
デルの渦粘性を分子粘性に足し合わせることで，変動成分をモデル化している．
3.2 スパコン TSUBAME の全ノードによる計算
東京都心部の 10km 四方のエリアに対し，実際の建造物のデータを元に 1m 間隔の格子解像度で詳細な気流のシミュレーションを行った．東京工業大学

４．局所細分化格子
4.1 再帰的な格子細分化
流体解析は物体近傍に高い解像度の計算格子を必要とするため，物体周りの流れを解析する場合な
どでは，計算領域全体で均一な格子を用いると計算効率が悪くなる．木データ構造に基づく再帰的な格子の細分化により，高解像度が必要な領域に細かい格子を割当てる．2 次元格子では 1 回の領域分割で元の領域が 4 つに分割され，3 次元格子では 8 つに分割される．GPU 計算ではメモリアクセスを効率化するため，木構造の末端であるリーフに 2×2×2
ではなく 16×16×16 や 32×32×32 のような大き

図 3 1m 格子による 10km 四方の範囲の都市部の気流計算

伝熱2016 年 10 月

図4

- 24 -

木データ構造に基づく細分化格子

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

い格子を割り当てる．木データ構造のリーフの再帰深度である細分化レベルが高いほどリーフの格子
解像度は高くなる．図 4 では，リーフ 5, 6, 7, 8 はリーフ 2, 3, 4 の 2 倍の格子解像度（格子間隔は 1/2）となる．
4.2 細分化レベルごとの時間刻み
格子ボルツマン法では，細分化レベルごとの格子幅に合わせて速度分布関数の 1 ステップ分の移動量を変更する必要がある．速度分布関数の速度 c =Δx/Δt は一定であるため，格子幅Δx に応じて時間刻み幅Δt を変化させる．最も細かい格子の細分化レベルを lm，時間刻み幅をΔt (lm)とすると，細分化レベル l の格子における時間刻み幅は，以下のよ
うに定められる．

動で記述することは，プログラムの生産性を考えると現実的でない．条件分岐の無い実装をするため，
全隣接パターンの GPU カーネル関数を記述し，複数カーネルの同時実行による高速化を図る．プログ
ラムの生産性を維持しながら記述するため，C++テンプレートを用いて 702 通りの GPU カーネル関数を自動生成する．すべての格子点を単一の GPU カーネルで計算する単純な実装に比べて 10 倍以上の高速化を達成している[7]．

Δt (l) = 2 lm−lΔt (lm)
細分化レベルが 1 段上がる毎に時間刻み幅は半分になり，同時に時間ステップ数は 2 倍になる．さらに，通常の細分化格子における空間補間および時間補間に加えて，偏差応力テンソルの連続性を保つために次式の細分化レベルの変換[6]を行う．
,

2

1
1

図 5

,

1
1
,
,
2
1
左辺はレベル変換により求められる速度分布関数である．上付き添字の1, ,1は細分化レベルを表し，チルド記号は時間・空間補間された値であることを表す．
4.3 局所細分化格子における計算の高速化
格子ボルツマン法に局所細分化格子を導入すると細分化の末端のリーフに割り当てる格子サイズが小さくなるため，GPU 計算の実行性能が大きく低下する．特に D3Q27 型では，リーフの外殻部分にある格子（外殻格子）の計算は，隣接格子点参照
に対して多数の条件分岐を含む．
外殻格子計算の高速化のためには，存在しうるア
クセスパターンをすべて列挙し，複数のカーネル関数に分割して記述する方法が考えられる．しかし，
D3Q27 格子ボルツマン法においては，外殻格子点の位置によって図 5 のようにパターンが 26 通りあり，そのそれぞれに 27 通りの隣接点参照が存在するため，それらの組合せによるアクセスパターンの総数は 702 に達する．これほど多数のパターンを手伝熱2016 年 10 月

- 25 -

隣接点参照のパターンに基づくリーフ内の格子の分類

4.4 走行中の自転車周りの流れの計算
本論文で提案する高速化手法を適用した例として，自転車競技の選手の周りの流れの計算を行った．
計算領域の中で自転車競技の選手と自転車は静止した物体として扱い，自転車の進行方向から一定速度で風が流入すると仮定している．このような空力計算では，計算境界において圧力波が反射し，計算精度を悪化させることが知られており，精度良く計算するためには計算領域を物体の遠方までとる必要がある．等間隔格子では広い計算領域をとるために物体近傍で格子解像度が不足するおそれがあるが，局所細分化格子を用いることで物体近傍を局所的に高解像度にすることができる．各リーフには32×32×32 の格子を割り当てる．局所細分化格子の最深の細分化レベルを 3 とし，物体近傍に高解像度格子を集中させた．最高の格子解像度は，計算領域全体に等間隔格子を割り当てた解像度（1024×
1024×1024）に相当し，格子点数は 5%にまで削減されている．
複数 GPU で計算するには，計算負荷が均等になるように領域分割する必要がある．リーフを空間充填曲線で辿り，各領域でリーフの数がほぼ等しくな

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

るように切断する．ヒルベルト曲線は交差しないため，モートン曲線などよりも分割領域が良好な形状となる．
本計算では，流体の計算と同時に速度場に沿って移流する煙状のトレーサ粒子を流体計算と伴に計算し，計算結果を可視化している．乱流が十分に発達した時刻のスナップショットを図 6 に示す．レイ
ノルズ数 40 万の流れになっているため，衝突ステップに Multiple Relaxation Time（MRT）モデル[8]を適用し，ここでも LES モデルとしてコヒーレント
構造スマゴリンスキーモデル（CSM）[4]を導入して
い る ． 自 転 車 の 境 界 条 件 に は Bouzidi ら のInterpolated bounce back スキーム[9]を適用している．
流れを可視化するためにトレーサ粒子を発生させ，自転車の背後に乱れた流れ場が生じていることが確認できる．

るため，GPU スパコンを使い高速な計算を行う．
固体粒子－流体間の直接相互作用計算を行うた
めには，粒子の運動を格子ボルツマン法に移動境界条件として組み込む必要がある．ここでは，各計算セルの固体体積占有率に基づいて流体と固体の重み付けをする Immersed Moving Boundary[12]を用いた．また，時間積分には安定性の高い MRT（Multiple
Relaxation Time）モデルを用い，移動境界条件を含む時間発展方程式は,,となる． は速度分布関数， はモーメント， は変
換行列， は緩和係数行列， は外力である．添字は速度分布関数の方向，eqは平衡状態の速度分布関数を表す． は固体体積率に基づく重み係数の行列，は単位行列， は移動境界の影響を与える追加項であり，を壁の速度として
Ω

図 6

走行中の自転車周りの流れのトレーサ粒子
による可視化

５．固液・固気の混相流シミュレーション
5.1 固体と流体の直接相互作用
これまでの殆どの固液・固気の混相流解析では，固体粒子は大きさを持たない点として扱われ，流体から固体粒子には理想化された抗力が働き，流体側には粒子運動をモデル化した力の項が運動方程式に加えられ，計算が行われてきた[10, 11]．流体を格子ベースの手法で計算する場合，1 格子内に複数個の粒子が存在することができる．しかし，このような流体と固体の相互作用の計算方法では，粒子が球形でない場合や高レイノルズ数の流れになった場
合，十分な精度で解析が行えなくなる．そこで，固体粒子のサイズより格子間隔を十分に細かく取り，
固体と流体が固体境界を介した直接相互作用を行う大規模計算を行う．必要とする計算資源が多くな

伝熱2016 年 10 月

- 26 -,,と与える．添字 ′は反対方向の速度分布関数を表す．
固体粒子の形状は，符号付きの距離関数で表現する．
固体粒子に作用する流体力は，固体粒子に覆われている格子点における運動量変化の総和を計算することで求める．
固体粒子の計算は，粉体計算でよく用いられている個別要素法[13]で行う．固体粒子間の接触力はバネとダンパでモデル化される．非球形の粒子や構造物の接触相互作用の計算には，粒子を剛体連結して形状を表現する剛体計算手法[14]を導入する．
5.2 流動層の大規模直接計算
床面から一様に流体を流入させて容器内の粉体を流動化する流動層について，粒子－流体間の直接相互作用を格子ボルツマン法と個別要素法を連成させ流動層解析を行う．
流動層の解析領域は 9.6 mm×256 mm×512mm とし，直径が 1 mm の球形粒子を 56 万個含む．粒子
直径に対して格子ボルツマン法の計算格子を 10
格子割り当てる．格子ボルツマン法の計算には 12.6
億 （ 96 × 2,560 × 5,120 ） 個 の 格 子 点 を 用 い ，TSUBAME 2.5 に搭載されている 128 台の GPU を利用した．初期に計算領域の下半分に粒子をほぼ稠密に配置し，床面全体に流入境界条件，上面全体に流出境界条件を課し，領域内の固体粒子群を流動化

J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

る．イチョウの葉の運動に合わせてレベルセット関数の座標変換を行い，格子ボルツマン法の移動境界条件における，各計算セルでの固体体積占有率の計算に用いる．レベルセット関数のゼロ等値面に球形粒子を配置することで，物体形状を表す粒子群デー
タを作成し，他のイチョウの葉との衝突を計算する．
TSUBAME2.5 の 128 台の GPU で 1024×1024×2048
格子を用い，512 枚のイチョウの葉が舞い落ちる過程を計算した[15]．図 10 に示すように，回転しながら落下するイチョウの葉もあり，特徴的な現象を再現できた．イチョウの葉同士の接触を考慮することで床に堆積する様子を計算することができた．

する過程の計算を行った．
図 7 は 128GPUで計算した領域全体のスナップ
ショットである．
粒子の初期位置に応じて層状に色付け（下から紺→紫→青→緑→黄→赤）して可視化している．図 8 の左側は図 7 の中の
1 つの GPU が計算を担当する領域の粒子と流れ場を可視化したもの，右図は粒子周りの流れ場を流
図 7 128GPU，12.6 億格子を用
速ベクトルで可いた 56 万個の球形粒子を視化したもので含む流動層の直接計算ある．粒子周りの微視的な流れも捉えられていて，高精度な解析を実現できている[15]．

図 8 1 つの GPU が担当する範囲である図 7 の 1/128の計算領域（左図）と，その一部分の粒子と流れの場の拡大図（右図）

2016 年 10 月

イチョウの葉の計算モデルの作成過程

図 10

512 枚のイチョウの葉の舞い落ちる計算

６．おわりに
最近，再び注目されるようになってきた格子ボルツマン法は，計算アルゴリズムが単純であり，高い演算性能を持つ GPU および GPU スパコンでの計算に非常に適している．大規模な流体計算の例を示したが，今後のスパコン性能の向上に伴い，さらなる大規模化の見通しは明るい．

5.3 舞い落ちるイチョウの葉のシミュレーション
小さな球形粒子を連結して複雑形状の物体を表現することにより，比較的容易に非球形粒子と流体の相互作用を計算することができる．物体から流体に及ぼす力は 4.1 節で述べた方法をここでも適用できる．例として，舞い落ちるイチョウの葉のシミュ
レーションを示す．図 9 のように，実際のイチョウの葉から 3 次元のポリゴンデータを作成し，ポリゴ
ンデータを符号付きのレベルセット関数に変換す

伝熱

図9

謝 辞
本研究では，国立研究開発法人海上技術安全研究所の小野寺直幸氏（東京工業大学在職当時），東京- 27 -J. HTSJ, Vol. 55, No. 233

特集：LBM の最近の進展

工業大学大学院生の渡辺勢也氏，長谷川雄太氏に大変協力して頂いた．慶應義塾大学の小林宏充氏にはLES モデルについて多くの助言を頂いた．本研究は科学研究費補助金・基盤研究(S) 課題番号 26220002「ものづくり HPC アプリケーションのエクサスケールへの進化」，科学技術振興機構 CREST「ポストペタスケール高性能計算に資するシステムソフト
ウェア技術の創出」，学際大規模情報基盤共同利用・共同研究拠点，および革新的ハイパフォーマンス・
コンピューティング・インフラから支援を頂いた．
記して謝意を表す．
参考文献[1] 青木尊之, 額田彰著, はじめての CUDA プログ
ラミング, 工学社 (2009)．
[2] http://www.openacc.org/[3] http://www.top500.org/[4] Kobayashi, H., The subgrid-scale models based oncoherent structures for rotating homogeneous
turbulence and turbulent channel flow, Phys. Fluids

boundaries, Physics of Fluids, 13-11 (2001) 3452–
3459.
[10]Yuu, S., Umekage, T., and Johno, Y., Numericalsimulation of air and particle motions in bubbling
fluidized bed of small particles. Powder Technology,

Powder technology, 77-1 (1993) 79-87.
[12]Noble, D. R. and Torczynski, J. R., A latticeBoltzmann method for partially saturated
computational cells, International Journal of Modern
Physics C, 9-08 (1998) 1189-1201.
[13]Cundall, P. A. and Strack, O. D., A discrete
numerical model for granular assemblies,
Geotechnique, 29-1 (1979) 47-65.
[14] 田中正幸, 酒井幹夫, 越塚誠一, 粒子ベース剛体 シ ミ ュ レ ー シ ョ ン と 流 体 と の 連
成,Transactions of JSCES, 20070007 (2007).
[15]Watanabe, S., Aoki, T., Hasegawa, Y., Large-scale
Simulations for Fluidization using Coupled Lattice
Boltzmann Method and Discrete Element Method on
a GPU Supercomputer. ECCOMAS Congress, Crete
Island, Greece (2016).

147-1 (1980) 219–228.
[7] 長谷川雄太，青木尊之, Octree 型 AMR を導入した格子ボルツマン法の C++テンプレートを用いたカーネル生成による GPU 計算の高速化, 情報2016 年 10 月

Mathematics with Applications, 69-6 (2015) 518–
529.
[9] Bouzidi, M., Firdaouss, M. and Lallemand, P.,
Momentum transfer of a Boltzmann-lattice fluid with

110-1 (2000) 158-168.
[11] Tsuji, Y., Kawaguchi, T. and Tanaka, T., Discrete
particle simulation of two-dimensional fluidized bed,

17 (2005) 045104.
[5] 小野寺直幸, 青木尊之, 下川辺隆史,小林宏充,
格子ボルツマン法による 1m 格子を用いた都市部 10km 四方の大規模 LES 気流シミュレーショ
ン. 情報処理学会ハイパフォーマンスコンピュ
ーティング研究会主催 HPCS シンポジウム(2013)[6] Olga, F. and Dieter, H., Grid Renement for LatticeBGK Models, Journal of Computational Physics,

伝熱処理学会論文誌コンピューティングシステム，
9-2 (2016) 34–45.
[8] Suga, K., Kuwata, Y., Takashima, K. and Chikasue,
R., A D3Q27 multiple-relaxation-time lattice
Boltzmann method for turbulent ﬂows, Computers &- 28 -J. HTSJ, Vol. 55, No. 233

行事カレンダー

本会主催行事開催日
2016 年
11 月
2(水)
～
4(金)
12

2017 年
3月

5月

行事名

申込締切

原稿締切

The Fourth International Forum on Heat
Transfer (IFHT2016)（開催場所：Sendai International Center）

2016.5.30
（延期後の
最終締切）

2016.7.29

2(金)講習会「計測技術」（開催場所：東京八重洲ホール 701 会議室）

2016.11.18

26(日)
～
30(木)

The First Asian Conference on Thermal
Sciences (ACTS2017)（開催場所：ICC Jeju, Korea（韓国済州島国際会議場））2016.9.21（Abstract
提出締切）

23(火)
～
25(木)

第 54 回日本伝熱シンポジウム（開催場所：大宮ソニックシティ）

本会共催，協賛，後援行事
開催日
行事名2016 年
11th Asian Thermophysical Properties
10 月
2(日)
Conference(ATPC2016)～（開催場所：パシフィコ横浜アネック
6(木)
スホール）
10 月
8(土)
可視化情報全国講演会（日立 2016）～（開催場所：茨城大学日立キャンパ
9(日)
ス）
10 月
13(木)
修善寺ワークショップ 2016～
（開催場所：ラフォーレ修善寺 研修14(金)センター（静岡県伊豆市））
10 月
18(火)
日本地熱学会平成 28 年学術講演会～（開催場所：群山市中央民館）
22(土)
10 月
22(土)
熱工学コンファレンス 2016～
（開催場所：愛媛大学）
23(日)
10 月
26(水)
第 44 回日本ガスタービン学会定期講
～
演会
27(木)
（ホテルリッチ＆ガーデン酒田）
The 5th International Conference on
10 月
29(土)
Human-Environment System(ICHES2016～
Nagoya)11 月
2（水） （開催場所：名古屋大学）
11 月

1(火)
～
4(金)

11 月

7(月)
～
10(木)
8（火）

11 月

11 月

16(水)
～
17(木)

11 月

23(水)
～
25(金)

問合先
Organizing Committee Chair Taku Ohara
URL:http://ifht2016.org/index.html

2016 年
10 月号

日本伝熱学会（担当 大澤）
Tel: 03-3259-7919, Fax: 03-5280-1616
E-mail: office@htsj.or.jp

2016 年
10 月号

2016.12.31

Conference Secretariat
Ms. Han Lee (KSME)
Tel: +82-2-501-5305，Fax: +82-2-501-3649
E-mail: info@acts2017.org
http://acts2017.org/
Conference Co-chair（日本語での問合せ）
高田 保之（九州大学）
Tel: 092-802-3100, Fax: 092-802-3098
E-mail: takata@mech.kyushu-u.ac.jp

2016 年
7 月号

2017.1.13

2017.3.10

実行委員会委員長 飛原 英治
東京大学大学院 新領域創成科学研究科
人間環境学専攻 内Fax: 04-7136-4647
E-mail: symp2017@htsj-conf.org

2016 年
10 月号

申込締切

原稿締切

2016.3.2

問合先

掲載号

問合せ先 田口良広（慶応大学）
Tel: 045-566-1809
E-mail: tag@sd.keio.ac.jp

2016.6.3

2016.7.29

2016.9.5

2016.9.5

2016.7.11

2016.9.11

2016.6.24

2016.9.1

実行委員長 稲垣照美（茨城大学）
Tel: 0294-38-5044
E-mail: terumi.inagaki.mech@vc.ibaraki.ac.jp
エレクトロニクス実装学会 乃万こずえ
Tel: 03-5310-2010, Fax：03-5310-2011
e-mail：noma@jiep.or.jp
日本地熱学会 実行委員 窪田 ひろみ
E-mail: hiromi-n@criepi.denken.or.jp
熱工学コンファレンス 2016 実行委員会事務局
Tel＆Fax: (089)927-9718
E-mail: tedconf16@me.ehime-u.ac.jp
公益社団法人日本ガスタービン学会 事務局
Tel：03-3365-0095, Fax：03-3365-0387
E-mail：gtsj-office@gtsj.org
ICHES2016 Nagoya 事務局
E-mail:iches2016office@davinci.nuac.nagoya-u.ac.jp

The 11th International Symposium on
Advance Science and Technology
inExperimental
Mechanics(11th
ISEM'16-Ho Chi Minh, Vietnam)（開催場所：Saigon Vissai Hotel）
第 1 回 ICMC アジア会議―CSSJ50 周年記念合同国際会議

日本実験力学会事務局 小林
Tel: 025-368-9310, Fax: 025-368-9309
E-mail: office-jsem@clg.niigata-u.ac.jp

ステップアップ・セミナー2016
計算機の大規模化と数値計算法（開催場所：大阪科学技術センター）
第 346 回講演会 実務者のための流体解析技術の基礎と応用（各種シミュレーション技術の適用
事例紹介及びデモ展示付き）（開催場所：大阪科学技術センター）
第 54 回燃焼シンポジウム（開催場所：仙台国際センター）

日本機械学会 関西支部事務局
Tel: (06)6443-2073, Fax(06)6443-6049
E-mail:info@kansai.jsme.or.jp

伝熱 2016 年 10 月

掲載号

低温工学・超電導学会
Tel:(03)3818-4539, Fax:(03)3818-4573
E-mail:LDJ04246@nifty.com

日本機械学会 関西支部事務局
Tel: (06)6443-2073, Fax(06)6443-6049
E-mail:info@kansai.jsme.or.jp

2016.7.15- 29 -

2016.9.9

東北大学 流体科学研究所 小林研究室
Tel: 022-217-5273, Fax: 022-217-5323
E-mail: sympo53@combustionsociety.jp

J. HTSJ, Vol. 55, No. 233

行事カレンダー
11 月

28(月)
～
30(水)

第 37 回日本熱物性シンポジウム（開催場所：岡山国際交流センター）

2016.7.8

2016.8.19

12 月

3(土)

2016.7.1

2016.10.21

12 月

8(木)
～
9(金)
10(土)

第 19 回スターリングサイクルシンポ
ジウム（開催場所：宇都宮大学峰キャンパス）
第 6 回潜熱工学シンポジウム（開催場所：京都工芸繊維大学センターホール）
2016 年度 計算力学技術者（CAE 技術者）資格認定「1・2 級」認定試験
第 3０回数値流体力学シンポジウム（開催場所：タワーホール船堀）

2016.9.16

2016.10.21

2016.9.23

2016.10.31

日本流体力学学会 事務局
Tel:(03)3714-0427, Fax: (03)3714-0434
E-mail jsfm@rf7.so-net.ne.jp

2016.9.2

2016.10.24

講演論文担当事務局 榎本啓士（金沢大学）
Tel: 076-234-4370
E-mail：eno@t.kanazawa-u.ac.jp
参加担当事務局 小橋好充（北海道大学）
Tel: 011-706-6384
E-mail：kobashi@eng.hokudai.ac.jp

12 月

12 月

12(月)
～
14(水)

12 月

19(月)
～
20(火)

第２５回微粒化シンポジウム（富山国際会議場 大手町フォーラ
ム）

3月

19(日)
～
22(水)

7月

10(日)
～
2019 年
12 月
末日

The 6th International Symposium on
Micro and Nano Technology inExperimental Mechanics(ISMNT-6)（開催場所：Fukuoka (Japan)）
EU Green Gateway to Japan（開催場所：東京）

日本熱物性学会事務局（東北大学流体科学研究所分子熱流動研究分
野内）Tel・Fax: 022-217-5277
E-mail: jstp@microheat.ifs.tohoku.ac.jp
実行委員長 戸田富士夫（宇都宮大学）
Tel: 028-649-5330, Fax: 028-649-5330
E-mail: toda@cc.utsunomiya-u.ac.jp
問合せ先 日出間るり（神戸大学）
Tel: 078-803-6657
E-mail: hidema@port.kobe-u.ac.jp
日本機械学会 事業企画Ｇ 石澤
Tel: (03)5360-3506, Fax: (03)5360-3509
E-mail caenintei@jsme.or.jp

2017 年

伝熱 2016 年 10 月

実行委員長 高橋厚史（九大）
Tel: 092-802-3015
E-mail: takahashi@aero.kyushu-u.ac.jp
日本能率協会 産業振興センター
Tel: (03)3434-0998, Fax:(03)3434-8076
E-mail: noriyasu_kamioki@jma.or.jp- 30 -J. HTSJ, Vol. 55, No. 233

お知らせ

平成 28 年度 日本伝熱学会 学術賞・技術賞・貢献賞・奨励賞候補者推薦募集のお知らせ
日本伝熱学会では本年度も各賞の表彰を実施いたします．つきましては，下記の要領に従って本年度の学術賞，技術賞，貢献賞および奨励賞の候補者推薦募集を行いますので，自薦，他薦を問わず，多数ご応募下
さい．
１．対象となる業績（１） 学術賞の対象は，原則として，最近 3 年間の Thermal Science and Engineering 誌または Journal of
Thermal Science and Technology 誌 に掲載された，あるいは，最近 5 回の日本伝熱シンポジウムにおいて発表し，国内外で審査のある学術論文集に掲載された伝熱に関する優秀な研究を主体的に行った研究者とその共同研究者とします．なお，受賞対象研究課題名は，必ずしも論文題目と厳密に一致する必要はありません．
（２） 技術賞の対象は，公表された優秀な伝熱技術を開発した者とします．
（３） 貢献賞の対象は，日本伝熱学会が主催する事業の実施において，顕著な貢献を行った会員とします．
（４） 奨励賞の対象は，原則として，最近 2 回の日本伝熱シンポジウムにおいて優秀な論文を発表した若手研究者で，発表時に大学院生，またはこれに準じる者（大学卒業後 5 年以内の者）とします．
（５） 学術賞および奨励賞の対象資格は，原則として本会会員に限ります．
（６） 贈賞数は，学術賞 2 件程度，技術賞 1 件程度，貢献賞 2 名程度，奨励賞 4 件程度とします．
２．選考方法（１） 各賞の選考は，「表彰選考委員会」が「日本伝熱学会賞審査・選考方法内規」によって行います．
（２） 表彰選考委員会は，公募の他に，各賞の候補を推薦することが出来るものとします．
３．提出書類（１） 申請書・推薦書 1 通・ 申請書・推薦書の書式は，学会ホームページからダウンロードしてご使用下さい．
・ 学術賞，技術賞，貢献賞，奨励賞それぞれに申請書・推薦書の書式が異なりますので，ご注意く
ださい．
（２） 論文抜刷または技術内容参考資料
6 部 （貢献賞以外の賞）（３） 日本伝熱シンポジウム講演論文集抜刷 6 部 （学術賞，奨励賞のみ）
４．提出期限平成 28 年 12 月 16 日（金）５．提出先および問い合わせ先：〒223-8522 横浜市港北区日吉 3-14-1
慶應義塾大学理工学部 システムデザイン工学科 長坂E-mail： nagasaka@sd.keio.ac.jp
Tel & Fax: 045-566-1735

伝熱2016 年 10 月

- 31 -

雄次 宛

J. HTSJ, Vol. 55, No. 233

お知らせ

日本伝熱学会主催 第 4 回 国際伝熱フォーラム
International Forum on Heat Transfer (IFHT2016)
日本伝熱学会は，世界最大かつ最もアクティブな伝熱の研究者・技術者集団として，世界中の研究者が集い，新しい伝熱研究のトレンドを生み出す国際伝熱フォーラムを開催します．2004 年に京都でフォーラムがスタートし，2008 年には東京，前回，2012 年には長崎で第 3 回のフォーラムが開催され，毎回 150 件を超える発表で盛り上がりを見せています．第 4 回の仙台では，欧米アジア各国から，およそ 223 件の発表が予定
され，伝熱分野の先端的研究と今後の研究動向を把握できるよい機会と期待されています．本フォーラムでは，日本伝熱学会創立 50 周年を記念して設立された抜山記念国際賞（第 3 回）を受賞された店橋護 氏（東京工業大学 教授）の受賞講演と授与式が予定されている他，伝熱分野の最先端を走る 6 名の研究者によるキーノート講演も予定されております．多くの方々に是非，国際伝熱フォーラムに参加いただけますよう，お願い申し上げます．
開催日時

2016 年 11 月 2 日（水）〜 4 日（金）参 加 費

事前申込み（9 月 30 日まで）
一般：50,000 円，学生：20,000 円
事前申込み（10 月 1 日以降）
一般：60,000 円，学生：30,000 円※参加登録 Web サイト http://ifht2016.org/registration/index.html（当日会場でもお申込みを受け付けております．）

会

仙台国際センター（http://www.aobayama.jp/）〒980-0856 仙台市青葉区青葉山無番地（仙台駅から地下鉄東西線で 5 分）

場

抜山賞受賞講演
Prof. Mamoru Tanahashi (Tokyo Institute of Technology, Japan)「Impacts of DNS and Advanced Laser Diagnostics on the Next-Generation IC Engine Development」
キーノート講演
Prof. Yoshinori Itaya (Gifu University, Japan)「Non-homogeneous Heat and Mass Transfer Problems in Upgrading Processes of Low Rank Energy」
Prof. Gary Rosengarten (RMIT University, Australia)「Towards Next Generation Solar Thermal Collectors: Keeping Things Hot and Cold」
Prof. Young-Hoon Song (Korea Institute of Machinery & Materials, Korea)「Thermal Management Technologies for Control of Diesel Emission」
Prof. W. Q. Tao (Xian Jiao Tong University, China)「On The Beauty of Field Synergy Principle for Enhancing Convective Ｈeat Transfer」
Prof. Shohji Tsushima (Osaka University, Japan)「Understanding, Designing and Fabricating Fuel Cells and Batteries: Coupled Processes of Transport
Phenomena and Electrochemical Reactions」
Prof. Evelyn N. Wang (MIT, USA)
「Liquid-Vapor Phase-Change Via Nanoengineered Surfaces」
Web ページ
問い合わせ先最新情報は，http://ifht2016.org でご確認ください．
組織委員会委員長
小原 拓〒980-8577 仙台市青葉区片平 2-1-1
東北大学 流体科学研究所
Tel&FAX: 022-217-5277 E-mail: ifht2016@microheat.ifs.tohoku.ac.jp
実行委員会委員長
麓 耕二〒036-8561 弘前市文京町 3
弘前大学 大学院理工学研究科 知能機械工学専攻
Tel&FAX: 0172-39-3676 E-mail: kfumoto@hirosaki-u.ac.jp

伝熱2016 年 10 月

- 32 -J. HTSJ, Vol. 55, No. 233

お知らせ

日本伝熱学会主催講習会「計測技術」
開催日 2016 年 12 月 2 日（金）13.00～17.00
会 場 東京八重洲ホール 701 会議室
地図（http://yaesuhall.co.jp/accessmap/）
要旨

今後センサ社会が益々発展されると言われているため，これを機に『温度計測』を基礎から学ぶ必要がある方々，計測分野で著名な専門家の方々とディスカッションをすることでより広い知見を得たい方々を対象に講習会を開催する運びとなりました．計測の基礎知識（原理や不確かさ）を習得せずに，実測結果を見誤ってしまうことで設計検証に時間を要することはありませんか？ 計測ばらつきを小さくするために，どのような工夫をしていますか？ 本講習会では，講演者と参加者が一体となった参加型イベントにするため，設計業務や研究で悩んでいることを共有するディスカッションの時間も長く設
けております．

題目・講師
13:00 ～ 13:10／開催の挨拶
日本伝熱学会/産学交流委員長

近藤 義広（（株）日立製作所）

13:10 ～ 14:00／熱電対の正しい使い方（測定誤差の要因と精度向上への鍵）
田川 正人 教授（名古屋工業大学）
14:10 ～ 15:00／放射温度計の正しい使い方
中村 元 教授（防衛大学校）
15:10 ～ 16:00／有効熱伝導率の測定手法
大串 哲朗 技術顧問（（株）アドバンスドナレッジ研究所）
16:20 ～ 17:00／パネルディスカッション
パネラー：講演講師， 司会：近藤 義広
17:30 ～ 19:00／講師との技術交流会
定先着 70 名員

参加費

（資料集 1 冊，技術交流会参加費込み）
学会員 25,000 円，非会員 35,000 円，特別賛助会員企業の方 20,000 円，学生 15,000 円

資料集

本セミナーに参加された方 1 名につき 1 冊準備します．
資料集のみの販売はいたしません．

申込方法

伝熱学会の講習会案内ページ（http://www.htsj.or.jp/lecture20161202）をご覧ください．
また，参加費は申し込みフォームに記載の口座に事前にお振り込みください．
申し込み期限 2016 年 11 月 18 日

問合せ先

伝熱

（公社）日本伝熱学会（担当 大澤）／電話（03）3259-7919／FAX（03）5280-1616／E-mail：office@htsj.or.jp

2016 年 10 月

- 33 -J. HTSJ, Vol. 55, No. 233

お知らせ

第 54 回日本伝熱シンポジウム研究発表論文募集
第 54 回日本伝熱シンポジウム実行委員会
委員長 飛原英治
幹事 党 超鋲
開催日： 平成 29 年 5 月 24 日（水）～5 月 26 日（金）会 場： 大宮ソニックシティ（https://www.sonic-city.or.jp/）
所在地〒330-8669 埼玉県さいたま市大宮区桜木町 1-7-5 電話
アクセス JR・東武野田線大宮駅より徒歩約 3 分
講演申込締切：講演論文原稿提出締切：
事前参加申込締切：

048-647-4558

平成 29 年 1 月 13 日（金）平成 29 年 3 月 10 日（金）平成 29 年 4 月 7 日（金）ホームページ URL： http://htsj-conf.org/symp2017/【シンポジウムの形式】
 講演発表形式としてa) 一般セッション（口頭発表）b) オーガナイズドセッション（口頭発表）c) 学生および若手研究者を対象とする優秀プレゼンテーション賞セッションを実施します．
 １講演あたりの割当時間は，一般セッションでは 15 分（発表 10 分，個別討論 5 分）で，各セッションの最後に総合討論の時間（5 分×セッション内の講演件数）を設ける予定です．オーガナイズドセッ
ションについては，オーガナイザーの指示に従ってください．
 優秀プレゼンテーション賞セッションについては，本号掲載のお知らせ「優秀プレゼンテーション賞（第 54 回日本伝熱シンポジウム）について」をご参照ください．
【研究発表申込方法】
 講演者（登壇者）は日本伝熱学会会員（学生会員，賛助会員，特別賛助会員を含む）に限ります．講
演者が日本伝熱学会会員でない場合は，会員申し込みと会費納付を行ってください．
 原則としてウェブによる申込に限らせていただきます．本シンポジウムで準備するウェブサイトから必要なデータを入力してください．ウェブによる申込ができない場合には，実行委員会事務局にお問い合わせください．
 申込の際に，一般セッション，オーガナイズドセッション，優秀プレゼンテーション賞セッションのいずれで発表するかを選択してください．優秀プレゼンテーション賞セッションにお申込みの場合には，本号掲載のお知らせ「優秀プレゼンテーション賞（第 54 回日本伝熱シンポジウム）について」をご参照ください．
 発表申込時に，論文要旨（日本語で 200～250 字）を入力してください．講演論文集の抄録として科学技術振興機構（JST）のデータベースに登録します．
 講演発表申込は，講演者 1 名につき 1 題目とさせていただきます．
 講演の採否およびセッションへの振り分けは，実行委員会にご一任ください．
【講演論文集電子版】
 講演論文集として電子版のみを発行し，冊子版は発行しません．講演論文はウェブ上に一定期間公開
します．日本伝熱学会会員は講演論文電子版を事前にウェブから読むことができます．
 講演論文集電子版を参加者に配布します．なお，講演論文集電子版は参加できなかった日本伝熱学会
会員に希望があれば配布しますので，シンポジウム終了後に日本伝熱学会事務局にお申込みください．
 講演論文集の長さは 1 題目あたり A4 サイズで 4～6 ページです．講演論文の作成様式は，1 段組み 50
字×46 行です．
 講演論文原稿の執筆要綱は会誌 1 月号の「お知らせ」および本シンポジウムのホームページに掲載します．

伝熱2016 年 10 月

- 34 -J. HTSJ, Vol. 55, No. 233

お知らせ 講演論文原稿は，原則として PDF ファイルで提出していただきます．
 講演論文締切日までに提出されなかった講演論文は，講演論文集電子版およびウェブに掲載されません．十分にご注意ください．
【講演論文の公開日】
 講演論文は日本伝熱学会会員を対象にパスワードを配布し，シンポジウム開催日の 1 週間程前からウェブ上で公開します．従って，特許に関わる公知日はシンポジウム開催日よりも早くなりますので，ご注意ください．
 会員が当該講演論文を後日原著論文として学術論文誌等に投稿される場合に配慮し，本講演論文集が限られた範囲に配布されたものであることを明確にすべく，シンポジウム終了後，一定期間を経て講演
論文の公開を終了します．
【参加費等】
 シンポジウム参加費
会員一般
事前申込：10,000 円，会場申込：12,000 円
非会員一般
事前申込：14,000 円，会場申込：16,000 円
会員学生
事前申込： 5,000 円，会場申込： 6,000 円
非会員学生
事前申込： 6,000 円，会場申込： 7,000 円※講演をするために日本伝熱学会学生会員になる場合，参加費は 7,000 円の割引価格となります．これには，平成 29 年度の学生会員の会費が含まれます．
※特別賛助会員は 1 口につき 3 名，賛助会員は 1 口につき 1 名，参加費が無料になります．
※会員とは，日本伝熱学会会員のほか，協賛学協会の会員を含みます．
※講演論文集電子版は参加者全員に配布されます．
 講演論文集電子版
講演論文集電子版のみの販売：5,000 円（シンポジウム後の販売となります．）【懇親会】
 開催日：平成 28 年 5 月 25 日（木）
 会 場：大宮ソニックシティ
 参加費：
一般
事前申込：8,000 円，会場申込：10,000 円，同伴配偶者無料
学生
事前申込：4,000 円，会場申込： 5,000 円【ご注意】 講演申込の取消および講演論文原稿の差し替えは，シンポジウムの準備と運営に支障をきたしますのでご遠慮ください．
 論文題目と著者名が，講演申込時と論文提出時で相違しないように特にご注意ください．
 口頭発表用として実行委員会事務局が準備する機器は，原則としてプロジェクタのみとさせていただ
きます．パーソナルコンピュータは各自ご持参ください．
 参加費，懇親会費等は参加取消の場合でも返金いたしません．
 本シンポジウムに関する最新情報については，随時更新するホームページでご確認ください．
 その他，ご不明の点がありましたら，実行委員会事務局まで e-mail または FAX でお問い合わせください．
【お問い合わせ先】
第 54 回日本伝熱シンポジウム実行委員会事務局
東京大学 大学院新領域創成科学研究科 人間環境学専攻内Email: symp2017@htsj-conf.org
FAX: 04-7136-4647

伝熱2016 年 10 月

- 35 -J. HTSJ, Vol. 55, No. 233

お知らせ

優秀プレゼンテーション賞（第 54 回日本伝熱シンポジウム）について
日本伝熱学会 学生会委員会
委員長 桃木 悟
第 54 回日本伝熱シンポジウムでは，下記の要領で，学生および若手研究者を対象とした優秀プレゼンテーシ
ョン賞セッションを設けます．日頃の研鑽の成果を披露するチャンスとして，奮ってご応募下さい．
開催
日： 平成 29 年 5 月 24 日（水） シンポジウム第 1 日
発 表 形 式： 発表者 1 名あたり，発表内容をまとめた 1 枚のアピールスライド提出とポスタープレゼン
テーションを行う形式をとる予定です．アピールスライドは冊子体として配布されると共にスライドショーとして会場にて案内される予定です．詳細については，決定し次第，シンポジウムのホームページに掲載いたします．
対象： 企業・大学・研究機関等の技術者・研究者で，平成 29 年 3 月 31 日現在で 28 歳以下の方，
または，申込当日に学生である方（ただし，社会人大学院生を除く）．
応 募 資 格： 発表者は日本伝熱学会の会員（正員・学生員）に限ります．発表者が日本伝熱学会会員でない場合は，講演論文原稿提出までに，会員申し込みを行ってください．なお，本セッシ
ョンで発表する方は，応募資格を必ず満たす必要があります．また，過去に本賞を受賞された方は応募することはできません．
応 募 件 数： 指導教員または研究グループ長等あたり，1 名（1 件）とします．
応 募 方 法： 第 54 回日本伝熱シンポジウム発表申込時に，本誌掲載の研究発表申込方法に従って，“優秀プレゼンテーション賞”を選択し，“指導教員または研究グループ長等の氏名”を入力
してください．なお，講演論文原稿の様式については一般セッションと同様のものとします．
審査・選考方法： 複数名の審査員による評価に基づいて選考し，受賞者を決定します．
表彰： 受賞者はシンポジウム第 2 日の学会総会で表彰されます．

伝熱2016 年 10 月

- 36 -J. HTSJ, Vol. 55, No. 233

事務局からの連絡

編集出版部会からのお知らせ

一各種行事・広告などの掲載についてー
インターネットの普及により情報発信・交換能力の比類ない進展がもたらされました．一方，ハードコピーとしての学会誌には，アーカイブ的な価値のある内容を手にとって熟読できる点や，一連のページを眺めて全貌が容易に理解できる点など，いくら電子媒体が発達してもかなわない長所があるのではないかと思います．ただし，学会誌の印刷・発送には多額の経費も伴いますので，当部会ではこのほど，密度のより高い誌面，すなわちハードコピーとしてぜひとも残すべき内容を厳選し，イン
ターネット（HP：ホームページ，ML：メーリングリスト）で扱う情報との棲み分けをした編集方針を検討いたしました．
この結果，これまで会告ページで取り扱ってきた各種行事・広告などの掲載につき，以下のような方針で対応させていただきたく，ご理解とご協力をお願いする次第です．
対 象

対 応

具体的な手続き（電子メールでの連絡を前提としています）

本会（支部）主 無条件で詳細を，会誌と HP
催による行事に掲載，ML でも配信

申 込 者 は ，記 事 を 総 務担 当 副 会 長補 佐 協 議 員（ ML 担当），広報委員会委員長（HP 担当）あるいは編集出版部
会長（会誌担当）へ送信してください．

条件付き掲載会誌：1 件当たり 4 分の l ページ程度で掲載（無料）
HP：行事カレンダーに掲載
してリンク形成（無料）ML：条件付き配信（無料）申込者は，まず内容を説明する資料を総務担当副会長補佐
協議員に送信してください．審議の結果，掲載可となった
場合には総務担当副会長補佐協議員より申込者にその旨通知しますので，申込者は記事を編集出版部会長（会誌担当）と広報委員会委員長（HP 担当）に送信してください．

会誌：掲載せず
HP：条件付き掲載 （無料）
ML：条件付き配信（無料）申込者は，公募内容を説明する資料を総務担当副会長補佐
協議員に送信してください．審議の結果，掲載可となった
場合には総務担当副会長補佐協議員より申込者にその旨通知しますので，申込者は記事を広報委員会委員長（HP 担当）に送信してください．

関係学会や本会
会員が関係する組織による国内
外の会議・シンポジウム・セミナー
大学や公的研究機関の人事公募（伝熱に関係のある分野に限
る）

一般広告
求人広告

申込者は，編集出版部会長（会誌担当）または広報委員会
委員長（HP バナー広告担当）に広告内容を送信してくだ
さい．掲載可となった場合には編集出版部会長または広報
会誌：条件付き掲載（有料）委員会委員長より申込者にその旨通知しますので，申込者
HP：条件付き掲載
は原稿を編集出版部会長または広報委員会委員長に送信し（バナー広告のみ，有料）てください．掲載料支払い手続きについては事務局からご連絡いたします．バナー広告の取り扱いについては
http://www.htsj.or.jp/wp/media/36banner.pdf をご参照下さい．

【連絡先】
・総務部会長：大竹浩靖（工学院大学）：ohtake@cc.kogakuin.ac.jp
・編集出版部会長：廣田真史（三重大学）：hirota@mach.mie-u.ac.jp
・広報委員会委員長：元祐昌廣（東京理科大学）：mot@rs.tus.ac.jp
・総務担当副会長補佐協議員：伏信一慶（東京工業大学）：fushinobu@htsj.or.jp
・事務局：大澤毅士・村松佳子：office@htsj.or.jp【注意】・原稿は Word ファイルまたは Text ファイルでお願いします．
・HP はメンテナンスの都合上，掲載は申込月の翌月，また削除も希望掲載期限の翌月程度の時間遅れがあることをご了承願います．
・ML では，原則としてテキスト文の送信となります．pdf 等の添付ファイルで送信を希望される場合 はご相談ください．
伝熱 2016 年 10 月

- 37 -J. HTSJ, Vol. 55, No. 233

事務局からの連絡
55 期入会（2016.6.9～2016.9.7）正 6 名
資
正
正
正

氏 名
杉井 泰介
齊藤 克弘
濱﨑 元貴

所 属
（株）日立製作所 研究開発グループ
三菱重工業株式会社
東芝キヤリア株式会社

資
正
正
正

氏 名
小倉 正大
坂本 明洋
竹山 亮

所 属
マイクロ・イクイップメント（株）
新日鐵住金株式会社 技術開発本部
株式会社フロロコート「伝熱」2016 年 4 月号の入会者リスト（p. 69）に，下記の方のお名前が抜けていました．
ここに謹んでお詫びし訂正申し上げます．
資
学生

氏 名
曽我 宏輔

伝熱 2016 年 10 月

所 属
東京工業大学 大学院 理工学研究科- 38 -J. HTSJ, Vol. 55, No. 233

編集出版部会ノート
Note from the JHTSJ Editorial Board
10 月号の特集は「LBM の最近の進展」と言う題目で，格子ボルツマン法（LBM）の手法や応用計算についての最新の研究動向について取り上げました．LBM は一般的に計算格子生成とは無縁
で，並列計算に適しているなど大変魅力的な流体の数値計算法ですが，1990 年ごろに登場したまだ新しい手法であるため，有限差分法や有限体積法，有限要素法などに比べると研究の蓄積量は格段に少ないといえます．したがって，モデルやプログラミング技術などに多くの研究テーマがあり，国内外の研究者が盛んに取り組んでいます．本号では，LBM による乱流や乱流伝熱，気液混相流の解析，さらに大規模並列計算機を活用した最新の研究内容について，第一線で活躍されている研究者の方々に執筆をお願いしました．伝熱研究に従事されている読者の方々に，LBM の魅力と可能性，そして何がどこまで計算できるのかについて現状を知って頂ければ幸いです．
最後に，本特集の内容に快く賛同頂き，多忙な時期にもかかわらず，締切日を遵守してご寄稿頂
きました著者の皆様に深く御礼申し上げます．

須賀一彦（大阪府立大学）

Kazuhiko Suga (Osaka Prefecture University)e-mail: suga@me.osakafu-u.ac.jp

副会長

長坂 雄次（編集出版担当，慶應義塾大学）部会長

廣田 真史（三重大学）

委員
（理事）

梅川 尚嗣（関西大学），桑原 不二朗（静岡大学），星 朗（東北学院大学）（協議員） 須賀 一彦（大阪府立大学），多田 幸生（金沢大学），土屋 敏章（富士電機（株））
中山 浩（中部電力（株）），西村 顕（三重大学），西村 伸也（大阪市立大学）
服部 博文（名古屋工業大学），村上 陽一（東京工業大学），
TSE チーフエディター

花村 克悟（東京工業大学）

編集幹事

伏信 一慶（東京工業大学）

編集出版部会 連絡先：

〒514-8507 三重県津市栗真町屋町 1577
三重大学大学院工学研究科機械工学専攻
廣田 真史
Phone & Fax: 059-231-9385 / E-mail: hirota@mach.mie-u.ac.jp

伝熱2016 年 10 月

- 41 -J. HTSJ, Vol. 55, No. 233