解説◊◊◊◊◊◊解説◊◇◇◇◇◇

固気混相流の数値シミュレーション—九州工業大学情報科学センター研究用計算機システムを利用した研究事例紹介—

湯 晋 ー 1．梅景俊彦 2r

1 はじめに
著者らの研究室では，固気混相流すなわち気流中に各種温度で固体粒子が混在した流れ場において粒子群と気流の力学的相互干渉によって生起される複雑現象のメカニズム解明を目指して，数値シミュレーションを中心とした研究（粒子を含まない流体のみの単相の流れ場や粒子濃度が非常に高い粉粒体層に関する研究も含まれる）を行っている．数年前まで著者らは主に九州大学大型計算機センターのベクトル計算機を利用してそれらの大規模計算を実行してきた．
ター研究用計算機システム更新を機に，しかし，平成 8年度に行われた前回の本学情報科学センPVM(ParnllelVirtualMachine)を用いた計算プログラムの並列処理化に取り組み，現在ではかなりの数値シミュレーションを本学の研究用計算機システムで実行してい(）年度）の CPU利用実韻は約 2(）（）（）（）時間であり，本年度（平成 11年度）もほぼ同程度のる．昨年度（平成 1利用を見込んでいる．本稿では平成 8年度以降を対象として，著者らが行っている研究で本学情報科学センターの研究用計算槻システムを利用した事例の中から，その一部を簡単に紹介する．

r、

2 著者らの研究室における研究用計算機システムの利用状況
最近の著者らの研究室における種々の研究テーマの中から，本学情報科学センター研究用計算機システ
ムを利用した代表的なものとして以下の①から◎を挙げることができる．
①自由噴流の乱流遷移過程の数値解析と実験による検証② Largじ EddySimulationを用いた高レイノルズ数固気混柳翡流中の気流と粒子の運動の数値解析と実験による検証③衝突による粒子間相互作用を考慮した高濃度固気混相噴流中の気流と粒子の運動の数値解析と実験による検証
l工学部機械知能工学科，

工学郭愕械矧能丁学科，

2yuuii:mech.
kyutech.
ac.
jpumekage(•.mech.k~·utech.ac.jp‑ 17← -~~“•

九州工業大学情報科学センター

広報第 12号 2000.
3

解説① DirectSimulationofMonte‑Carlo法を用いた乱流及び循環流動層の数値解析と実験による検証⑤離散要素法 (Di:‑;t.
inctElementMetho<l)を用いた微小粒子気泡流動層の数値解析と実験による検証

R Smoothe<lParticle法を用いた微小粒子気泡流動層の数値解析と実験による検証⑦ Smoothe<lParticle法を用いた流動粉粒体層の数値解析と実験による検証
R離散要素法を用いた粉体層の流動化現象（フラッシング現象）の数値解析と実験による検証⑨離散要素法を用いた粉体層の摩擦現象の数値解析と実験による検証◎液中粒子群の沈降及び圧密機構の数値解析と実験による検証◎撹拌槽内の流体運動の数値解析と実験による検証

上記 1
1テーマの中で①から＠の数値計算にはいずれも本学情報科学センターの研究用計算機システム
を 使 用 し て い る ． こ の う ち ① か ら ④ の 4テーマが，PVMを用いた複数の CPUによる並列計算を行う.
3年の間に九州大学大型計算機センターのベクトル計算機から計算処理を移行させたことによって，最近 2テーマである．

゜

離散要素法を用いた計算（上記の研究テーマでは⑤．R及び⑨）については平成 8年度以前から一貫して本学情報科学センターの研究用計算機システムを利用してきた．離散要素法は高濃度の粒子群が多体接触している場合の粒子間相互作用を 1個 1個の粒子に着目して計算する手法であるが，離散要素法は構造的にプログラムの高ベクトル化率の達成が困難でベクトル計算機を使用する利点がほとんど無いことがその理由である．
⑥．⑦及び＠の各テーマは平成 8年度の前回の計算機システム更新後に着手した研究テーマで，当初から本学の計算機システムの利用を前提とした数値解析コードの開発を行ってきた．
＠のテーマは現在でも九朴1
大学大型計算機センターを利用している．撹拌槽内の流体運動の数値解析では複雑な形状をした撹拌翼周りの境界条件とその移動の取り扱い等にやや複雑な点があるため上述の①カ~ら①の各テーマに比べると取り組みは遅れているが，現在 PVMを用いてプログラムの並列処理化を行うことによって本学の計算機システムヘの移行を進めている．

3 固気混相流の数値シミュレーション事例の紹介
前節で挙げた①から①の各テーマにおける固気混相流の数値シミュレーションをその粒子濃度によって大別すると次の三つに分類できる．
(1) 粒子濃度が希薄で，粒子どうしの衝突が無視できる場合，すなわち流体の運動が支配的で，粒子・流体間の相互干渉項を含んだ Nav
ier‑Stoker‑,式と粒子の運動方程式を解く問題．流体のみの流れを対象としている①と◎を含めると，①．②及び◎がこの範疇に含まれる．

‑ 18‑

九州工業大学情報科学センター
広報第 12号 2000.
3

゜解説

(2)

粒子濃度が高濃度となり，流休の運動とならんで粒子どうしの衝突が運動を支配する問題．前述のテーマの中では③と①がこの範疇に含まれる．なお，気流の運動が大きな役割をしており，かつ粒子の濃度がより高濃度になってもはや 2体衝突では表せない場合として気泡流動層の流れ（＠及び＠）や粉体の流動化現象(R)があるが，それは次の (3)の場合に含めて記述する．
(3)

さらに粒子濃度が濃くなり，3個以上の粒子どうしの多体衝突やすれあいが現象を支配し，流体の運動は副次的となり，現象を支配する新しい基礎方程式を導出して，それを解く問題．⑤から＠がこの範疇に含まれる．
これらの分類に従って，それぞれの数値シミュレーションを簡単に見ていくと以下のようになる．
分類 (1)の低濃度固気混相流の研究の発展に最も大きく寄与したのは Na
vier‑Stoke
H式の直接数値計算（粒子・流体間の相互作用項を抗カ・揚力としてモデル化している場合を含む）と考えられる．ほとんどの低濃度固気混相流の運動を支配するのは乱流である．直接数値計算では現在でも低レイノルズ数乱流しかr計算できないが，それでも従来の不完全な乱流モデルに基づいて乱流を計算するよりもはるかに正確かつ詳細に計算できる．なお固気混相流の高レイノルズ数の流れについても著者らの研究室において現在すで

argeE<l<lySimulationの計算モデルが導出さにサブグリッドスケールの流れに粒子の存在を考慮した Lれており，高レイノルズ数固気混相乱流の計算も着実に進んでいる．図 1に高レイノルズ数 3次元噴流 (Re=1000(））の流れを LargeEddySimulationを用いて計算した結果（等渦度線図）と可視化実験結果の比較を示す．計算結果と実験結果はいずれも噴出してから約 9(）無次元時間後の流れであり，計算結果は噴出先端の形状をも含めて実験結果を良く表現しているといえる．図 2に噴流の横断面から見た気流のグリッドスケールとサブグリッドスケールの運動エネルギー分布の単相噴流と混相噴流の計算結果の比較を示す．グ
リッドスケール，サプグリッドスケールとも乱流発達領域では粒子は気流の運動エネルギーを減少させ，初期領域では増加させる．
次に (2)の場合，すなわち高濃度固気混相流で流体と粒子どうしの衝突の両者が運動を支配している場

1)の場合と同様に Na.
vier合であるが，この分野においても現実の流れの説明に最も大きく寄与したのは (‑Stoke
H式と粒子のラグラジアン形運動方程式を連成して解く方法であろう．このとき粒子どうしの衝突

rは粒子を 1訊 r<lparticleと考えて，すなわち 2体衝突と考えて，衝撃方程式によって衝突後の粒子の速度を求める．衝突時の粒子どうしの摩擦は考慮する．このとき，1)の場合と違って，最も大きな問題となるの（は粒子の個数の問題である．粒子の運動をラグラジアン的に 1個 1個計算するので，計算できる粒子の個(1)の場合は粒子濃度が低く，衝突が無視できるので粒子個数はほとんど問題とならな数には限界がある (い）．この問題を解決する方法として研究テーマ①で用いられているのが，統計的な代表粒子の衝突等を考えて計算する DirectSimulationofMonte‑Carlo法である．不均質な流れ場中の粒子軌跡の計算に代表粒子を用いることにはかなり大きな誤差が見込まれるが，気流の一計算メッシュ内に存在する代表粒子数を数百とすることにより，近似度を増すことは可能である．図 3に 3次元乱流流動層中の粒子の瞬時ベクトル線図 (a)及び瞬時粒子位置図 (b)の計算結果及び同一条件の実験による瞬時粒子位置図（（:)を示す．実験では白い所が粒子の多く存在している部分で，計算結果の図とは白黒が反転している．図 3中の (b)．
（（：)の比較から，実際の乱流流動層中の現象を計算結果は表現していることがわかる．
最後に (3)の場合，すなわち粒子濃度が非常に濃くなって粒子の衝突がもはや 2体衝突として取り扱え‑ 19‑

九州工業大学情報科学センター
広報第 12号 2000.
3

解説s・91.
.
-~-.主­言．ーコ—
ミ= ‑:‑ ． て゜

GArこ―‑＝ーェニfぎ

゜

g.9II

¥/10.
020.
030.
0(a)

4.
040.
0 4X/D

g9Ie

←‑

¥︒

5.9IR

10.
020）(
43000(b)り、 44040X/D

図 1: 高レイノルズ数 (Re=10000)3次元単相噴流の計算結果と実験結果の比較，（叫 L孔rgeE<l<lySinrn‑

0=37.
5m,js.D =4.0mm,t=9.
39x1(）
―:1s)
.
(b)可
lat.
ionを用いた等渦度線図の計算結果（空気． U
視化実験結果（水． U
o=2.0m/s.D =5.0mn．
，)t=3.
14X 10‑1s)‑ 2(）

九州工業大学情報科学センター
広報第 12号 2000.
.
3/

解説？
こ7

吋くJ:e。
.l
直マ‑l ――—_― Iol00ず￨

20.
030,0X/D(a)単梱噴流7
i0•F30.F¢

i'‑・・・Q5し9,D︑,1^9fI
4x[︒

40.
010、0r(b)混相噴流

瞬時グリッドスケール連動エネルギー線図

8．

0念．
． 8}．
}
s.‑
,．
．b・:

︱︱︱

N
a︑!O

0.Voo.<i

10.
0、 も.oto戸賃ミ4;
,. ？.

な1,

20.
0l I
40.
0 440

30.
0

X/D

(c)扉相噴流．

︱︱︱

Gnt

‑・←,

9.
.

￠a.

心＿，
ー、、
ヽ炉た↓，．
心ダざと．•・a: ．

＂へ⇒~心9

ヽ

:・7~.

.
・‑
←︒

3^rootv,r

ヤ

O,"

9 .

な ⇔ 一

l100

20.
0

30.
0

40.
04
4.
0
X/D

(d)混相嗜流

瞬時サブグリッドスケール運動 J．こネルギー線図図2: 粒子の存在を考慮した LargeEddySimulationを用いた高レイノルズ数 (Re=10000)
3次元単相噴流及び固気混相噴流における瞬時グリッドスケール及び瞬時サブグリッドスケール運動エネルギー分

7.
5rn/s,D = 4.
0rnm,t= 1.
28x10‑
2s)
布の計算結果（空気． U0= 3‑‑21‑

九州工業大学情報科学センター
広報第 12号 2000.
3

解説

2525

202020exー

5GAX

5iuaーex゜io

10u15ょYIDYID(c)(b)

図3:DirectSimulationofMonte‑Carlo法を用いた 3次 元 乱 流 流 動 層 の 計 算 結 果 と 実 験 結 果 の 比 較(Uo =6.0m/s,Dp=31011,m),(a.
)瞬時粒子速度ベクトル線図の計算結果 (t =0.
651s),(b)瞬時粒t=0.
651s)．（(：）可視化実験結果
子位置図の計算結果 (‑‑22‑九州工業大学情報科学センター
広報第 12号 2000.
3︒

18︒

(a)

48゜︒

A3︒

゜4YID

解説

なくなり，3休以上の多体衝突として考慮しなければならない流れの場合である．対象となる流れは各種粉粒体の流れ，気泡流動層，高混合比空気輸送等が考えられる．このような多体衝突の力学機構を，近似的にではあるが，簡単に計算可能にしたのが離散要素法 (DistinctElementMetho<l,以下 D.E.M.と略す）である．

D.E.M.を用いると，従来の連続体モデルでは到底計算できないと思われていた粉粒体特有の諸現象を，はとんどの場合は定性的にではあるが，容易に計算できるようになった．一例として，スラギン
グ現象が現われる気泡流動層中の粒子の瞬時位置図の数値計算結果と実験結果の比較を図 4に示す．粒子間相互作用は D.E.M.を用いて計算している．複雑な気泡の形態を計算結果は良く表現している．

D.E.M.の出現は Nav
ier‑Stokes式の直接数値計算と並んで固気混相流の計算力学の発展に大きく貢献してきたが，やはりそこには計算可能な粒子個数の限界という大きな問題点が残る．現在，

D.E.M.で計算できる最高の粒子数は球形 3次元粒子で通常約 10万個であるが，これは粒径が少し小さくなるとスプーン一杯分ほどの粒子数に過ぎない．

2体衝突の場合のように統計的代表粒子の衝突で多体衝突を表現することには多大の困難があり，連続体近似を再考することも一つの大きな道と思われる．連続体近似のrarticle法，略して S.
P.法であ有望な方法として現在著者らの研究室で取り組んでいる方法が SmoothedPる．連続体近似であるから，構成方程式，状態方程式が必要であるという基本的な困難は当然のことながら存在するが，S.P．法では運動方程式を常微分方程式に変換して，ラグラジアン的に解くので，数値解を簡単に得ることができる．粉体には粒子どうしの摩擦が入るので構成方程式はかなり複雑となる．従来の偏微分方程式で解く場合だとここで解けなくなる場合が非常に多かったが，ことができる．例えば，色々な構成方程式を導出し，S.
P.法だと比較的簡単に解く

S.
P．法によってどれが最も実際の現象をよく表現するかを見ることができるわけである．図 5にホッパーから流出する粉体の瞬時粒子位置図及び瞬時のベクト.P．法による計算結果と同一条件における実験による可視化図を示す．図から分かるように計算ル線図の S
結果は実験結果を良く表わしている．

4 おわりにr以上，著者らの研究室に於ける本学情報科学センター研究用計算機システムの利用状況と研究テーマについて簡単に紹介した．本稿の記述は平成 8年度のシステム更新以後を対象としたが，その期間だけに限定しても著者らは本学情報科学センター研究用計算機システムを利用して行った研究成果を 8編の論文 [1]‑[8]（投稿中及び投稿準備中の論文を除く）として公表している．平成 1
2年 4月から稼働する新システムではさらに研究用計算機の能力向上が図られるが，今後もそれを有効に活用して研究をさらに進展させていきたぃど考えている．

参考文献[1] 湯晋一，林昭久，脇昌裕．梅景俊彦：流動粉粒体（粗大粒子）の応力と歪み速度の関係一離散要素法と4巻． 4号， pp.212‑220(1997).
実験による導出ー．粉体工学会誌．第 3‑ 23―‑九州工業大学情報科学センター

広報第 12号 2000.
3

解説

3.
0[3.
0‑2.
5卜 12.0￨

2.5f，

i2.
0~}i,・

g1.5￨l。「><1.
5[ll・&

1.
0ト

1.0[￨

0.
5‑

0.
5,
̲̲

,~
3鼠軍 §:〉゜贔惑︒

゜︒﹂

0.o[0.0-•
I II I 9 9 10.
Lぷ

L0.
51.
0

Y/D(a)(b)h••, l I I L 」」

0:
51.
0

Y/D。=

図4: 離 散 要 素 法 を 用 い た ス ラ ギ ン グ 状 態 が 現 れ る 気 泡 流 動 層 の 計 算 結 果 と 実 験 結 果 の 比 較 (U=0.
4rn/s,Dp=31011,m),(a)瞬時粒子位置図の計算結果 (t 0.
48s),(b)可視化実験結果

24‑

九什I
工業大学情報科学センター
広報第 12号 2000.
3

解説

コ

s・ss9

g80L

麟鴫s.
oooos
g̀g8<s.88C'[aモ lA＼StooN.
.
.
.
̲・.
・;
.
.:
̲.
.
.
・̲
・.
.
'
令•.: ・•• ‑
•..
．
・・.
‑.
.
.
.
.
.
.
.
•,...
.
．
, ．．

.
‑.
,・,:・.
‑・.
:::.
.
.
.
'9.

;,•*．、; :,-:•.

':

魯．：：•'・ら‑・・~．
:・・．
，；．
.
.
.
‑,・:・・.
.•，ょ．． u ,，．．
疇｀

.,.,.

●●’•嗜’嶋・•、,．、,•̲:;;•'；・ ・•“←・、•',_.,..

・.
:.
.
.
』
,•;
••一·.,.',•・:;..
．
:.
.
.
‑.
―
:<・:
•に”•:・‑.
::.
^'1• :·~. ・•、・ー・ ・こ．
＇＂奪',,•r.
.'l.：錢忍玲．，・、•:＂'..,.....'..'{_:~_·
・；さ牡i.
1t;;·:.~_’’;；:t'滋硲．名ペ••心•,．．ュ:•．
．
•-噌．・疇 -t,•

•み・脅 9 ふ』 •9•し＇ー・●12(a)

ーx[mmlx!
mm]

xlm叫(b)︒

ー2l210OIO︒

IOし゜S(c)図5:SmoothedP
article法を用いたホッパーから排出される粉体の流れの計算結果と実験結果の比較（い=ずれも t 0.
20s)．（叫瞬時粒子位置図の計算結果，（ b)瞬時粒子速度ベクトル線図の計算結果．（（・）可視化実験結果‑ 25‑-•

九州工業大学情報科学センター
広報第 12号 2000.
3

解説[2] 湯晋ー．離散要素法による粉体流動の数値解析．日本機械学会 RC14(）混相流動解析プログラムの高度化
p.
55‑60(1997).と応用に関する研究分科会研究報告書． p[3] 湯晋ー．脇昌裕．岩政明．梅景俊彦． Smoothe<lParticle(S.
P.）法を用いた流動粉粒体層の数値解析と5巻． 3号． pp.
174‑
182(1998).
実験による検証．粉体工学会誌．第 3[4]Yun
1S.andT.
U
11
1ekage1Nu
mer
icalSimulationofA
ir叫ParticleMotionsinS
mallP
articleBub‑blingFluidizedBed,ProceedingsofPTF'98,TheM
ianliBeach1
998AIChEAnnualMeeting.
pp.
115‑
120(1998).
[5]Umekage,
T.
,S.
Yuu,T.
ShinkaiandT.Abe,Nu
mer
icalS
imu
lationforBlockageo
fCohesiv
ePa.
1
‑
ti
‑
desi
naHopperU
singt
heDistictElementM叫 10dandi
tsCorrelationwithExperimentalRe‑
sultsofRealCohesiveGranularMaterials,AdvancedPowderTechnology,Vol.
9,No.
4,pp.
331‑
344(1998).

[6]湯晋ー．河野浩幸，梅景俊彦高レイノルズ数スリットノズル固気混相噴流の流れに与える粒子の影響°

日本機械学会論文集 (B編)‑第 66巻． 641号． pp.
57‑
66(2000).

[7]Yuu,S.
,M.Waki,A.
IwarnasaandT.Umekage,NumericalSimulationoftheVelocityai
1dStressFieldsforaFlowingPonderU
singtheSmoothedP
articleMethodandE
xperimentalVerification.
AdvancedPowderTechnology,Vol.
11,No.
1(2000)［掲載決定］．
[8]Yuu,S.
,T.UmekageandY.
Johno,NumericalSimulationofA
irai
1dParticleMotionsinBubblingF
luidizedBedo
fS
mallP
articleH.PowderTechnology,(2000)［掲載決定］．

︒

‑26‑

九州工業大学情報科学センター
2号 2000.
3
広報第 1