渦法および粒子追跡法の内部流れ混相流解析への適用
Application of Vortex Method and Particle Trajectory Tracking Method into Numerical Simulation of Internal Multiphase Flows
2007 年 3 月
磯 良 行
目次
第1章 序論··· 1
1.1 緒言··· 1
1.2 流体－粒子の混相流解析に関する従来の研究··· 3
1.2.1 流体－粒子の混相流解析の概要とモデル化の視点··· 3
1.2.2 オイラー・ラグランジュ型の混相流解析··· 5
1.2.3 渦法と粒子追跡法によるラグランジュ・ラグランジュ型の混相流解析··· 7
1.3 本研究の目的··· 8
1.4 本論文の概要··· 9
第2章 渦法と粒子追跡法による混相流解析手法と数値モデル··· 12
2.1 はじめに··· 12
2.2 渦法と粒子追跡法によるラグランジュ・ラグランジュ型の混相流解析の概要·· 12
2.2.1 解析手法の特徴··· 13
2.2.2 解析のアルゴリズム··· 13
2.3 渦法による流れ解析··· 17
2.3.1 渦法による流れ解析の概要··· 17
2.3.2 基礎方程式··· 19
2.3.3 Biot-Savart法による速度場の計算方法··· 20
2.3.4 境界要素法による内部流れの定式化（間接法および直接法）··· 21
2.3.5 渦度場の離散化（渦要素モデルと壁面渦導入法）··· 34
2.3.6 渦度場の時間変化の計算方法··· 43
2.3.7 圧力場の計算方法··· 52
2.3.8 高速化··· 56
2.4 粒子追跡法による粒子挙動解析··· 60
2.4.1 粒子追跡法による粒子挙動解析の概要··· 60
2.4.2 粒子運動のモデル化と計算方法··· 62
2.4.3 並進運動··· 63
2.4.4 回転運動··· 66
2.4.5 衝突運動··· 67
2.5 相間相互作用の計算方法··· 70
2.5.1 粒子による流れ場の変化とモデル化の概要··· 70
2.5.2 粒子の並進運動による流れ場の変化の計算方法··· 71
2.5.3 粒子の回転運動による流れ場の変化の計算方法··· 75
2.5.4 粒子周りのミクロな流れの計算方法··· 76
2.6 まとめ··· 78
第3章 混相内部流への適用と検証··· 80
3.1 はじめに··· 80
3.2 二次元解析（One-way解析）··· 81
3.2.1 固液二相噴流の解析··· 81
3.2.2 混合層を伴うチャネル内での固液二相流の解析··· 88
3.2.3 鉛直チャネル内での固液二相流の解析··· 95
3.2.4 T字型合流路内での固液二相流の解析と実験··· 99
3.2.5 高速化手法の有効性と妥当性の検証···112
3.3 三次元解析（One-way解析）··· 120
3.3.1 一様流中への固液合流噴流の解析··· 120
3.3.2 固気二相噴流の解析··· 125
3.3.3 T字型合流路内での固液二相流の解析··· 131
3.4 相間相互作用を考慮した解析（Two-way解析）··· 137
3.4.1 固体粒子群の自由沈降流の解析（二次元解析）··· 137
3.4.2 固体粒子群の自由沈降流の解析（三次元解析）··· 142
3.4.3 堆積を伴う高濃度粒子流の解析と実験（二次元解析）··· 151
3.4.4 堆積を伴う高濃度粒子流の解析と実験（三次元解析）··· 163
3.5 まとめ··· 169
第4章 実用流体機器への応用··· 170
4.1 はじめに··· 170
4.2 氷蓄熱装置における熱交換器内での氷水二相流の解析··· 170
4.2.1 対象流れ場と解析目的··· 170
4.2.2 モデル化··· 172
4.2.3 解析条件··· 173
4.2.4 解析結果と考察··· 174
4.3 石炭焚き発電用ボイラにおける伝熱管群周りの灰粒子挙動の解析··· 178
4.3.1 対象流れ場と解析目的··· 178
4.3.2 モデル化··· 179
4.3.3 解析条件··· 181
4.3.4 解析結果と考察··· 182
4.4 まとめ··· 191
第5章 結論··· 192
5.1 本論文のまとめ··· 192
5.2 今後の展望··· 193
謝辞··· 195
参考文献··· 196
本論文を構成する公表論文··· 203
主要記号
w , v , u
u ：速度ベクトル
z y
x, ,
ω ：渦度ベクトル または 角速度ベクトル
z y
x, ,
α ：渦強度ベクトル
, ：循環量
：渦要素の粘性核半径（コア半径）
p ：圧力
：密度
：粘性係数 または 二重吹き出し または 摩擦係数
：動粘性係数
G ：Green関数
H ：Bernoulli関数
：デルタ関数 または 微小な長さ n ：法線方向の単位ベクトル
k ：渦要素の回転軸方向の単位ベクトル
：速度ポテンシャル
n ：速度ポテンシャル流束（ の法線方向微分値）
：吹き出し
V , v ：速度 または 体積
L , l ：長さ
T , t ：時間
S ：面積
R , r ：位置のベクトル または 距離（二点間距離）のベクトル
D , d ：直径
M, m ：質量 または 流体と粒子との質量混合比
I ：慣性モーメント
F , f ：力のベクトル
T ：回転トルクベクトル
J ：衝撃力（力積）のベクトル
g ：重力加速度ベクトル
e ：反発係数
CD ：抗力係数
CL ：揚力係数
CT ：流体粘性による回転トルク係数
CV ：体積濃度
lbp ：境界パネルの長さ（一辺の長さ）
Nbp ：境界パネル枚数
h ：壁面上の速度参照点高さ（壁面上の渦層高さ）
t ：時間刻み（タイムステップ）
Np ：一個の標本粒子（計算上扱う粒子）が代表する実際の粒子数
way
ltwo ：Two-wayカップリング用の計算格子幅
ν L V
Ref ：流れ場のレイノルズ数（慣性力と粘性力の比）
ν d V Vf p p
Rep ：粒子のレイノルズ数（慣性力と粘性力の比）
f f
f p
p
V L d 18 St
2
：ストークス数（粒子の応答時間と流れ場の特性時間の比）
添え字
下付き添え字 f ：流体相
下付き添え字 p ：粒子相（分散相）
下付き添え字 ave ：時間平均 上付き添え字 ’ ：時間変動
上付き添え字 * ：境界面上 または 衝突後
その他
単位 ：値が表示される都度に記す
第 1 章 序論
1.1 緒言
流れ場中に固体粒子，気泡，液滴などの分散した粒子が混じった流体－粒子系の混相 流は，輸送，攪拌，分離，反応，燃焼を伴う多くの工業装置内で観察され，工学的に非 常に重要である．例えば，固気混相流の工業装置としては，石炭焚きボイラ，流動層，
微粉炭燃焼器，サイクロン，集塵器，ガスタービン翼，固体ロケットエンジン，サンド ブラスト装置，などがある．固液混相流の工業装置では，スラリー輸送（固体粒子の水 力搬送装置），固形物汚泥の処理・沈殿装置，脱硫吸収塔，バイオリアクタ，氷蓄熱，高 レベル放射性廃液固化用のガラス溶融炉，などがある．気液混相流では，火力・原子力 発電，気泡塔，オゾン処理装置，バブリング攪拌槽，噴霧燃焼器，などが挙げられる．
上記のような混相流の中で，気相や液相の流体場に固相である固体粒子が混じった固 気・固液混相流では，固体粒子の非定常な挙動によって，摩耗・堆積や，粒子濃度の空 間・時間的な不均一による性能低下等がしばしば問題となる．最悪の場合，摩耗による 流路壁の損傷・破損や，堆積による流れの閉塞など，運転を停止せざるを得ない重大な トラブルに発展する可能性がある．そのため，上記のような混相流装置を設計・開発す る際には，初期段階において，粒子フローパターンや，粒子と壁面との干渉等の非定常 な流動特性を把握する必要がある．
これらの混相流装置の設計・開発には，理論・実験式，経験的なノウハウだけでなく，
実物スケール・サブスケールでの実証試験を行う場合が多い．ただし，ニーズの拡がり・
高性能化・低コスト化が強まる中，コストや期間を短縮する目的で，実証試験に代わり，
数値流体解析技術（Computational Fluid Dynamics, CFD）の適用が強く求められている．
また，数値解析を用いることで，計測では得ることが難しい粒子の濃度・速度分布など の定量的な情報を得ることができる．さらに，粒子解析にラグランジュ的手法（粒子追 跡法）を用いるならば，粒子が壁面へ衝突する際の位置・角度・速度等を個々の粒子に 対して算出できるため，摩耗や堆積を予測する上で有益なデータを得ることが可能にな る．産業界の現場においても，新規装置の設計や開発に数値解析を活用する事例が増加 している．近年，定常的な流れ現象については，市販の汎用CFDツールで概ね対応でき る環境が整いつつある．しかしながら，産業機器の多くは，高レイノルズ数で非定常な 渦流れを伴うため，この非定常現象に起因した性能低下やトラブルを，予め設計段階で 予測することが非常に重要となっている．
一般に，工学上問題となる流れの多くは高レイノルズ数の乱流であり，分散相である 粒子の挙動は，流れ場に存在する様々なスケールの乱流渦に支配される[1] - [5]．特に，噴 流，急拡大，合流流れ，鈍体物体後流などの様に，非定常性の強い渦流れが伴う混相流 では，乱流渦と粒子との干渉を解析することが重要となる．Croweら[1], [2], [3]は，粒子運
動が流れ場の組織渦への粒子の追従性により分類できると考え，粒子運動を評価するパ ラメータとして，組織渦の時間スケールに対する粒子の緩和時間の比であるストークス 数を用いることを提案した．その後，二次元混合層，噴流，物体後流などの様々な流れ における計測 [4], [5]によって，Croweらの評価方法の有用性が確認された．
このように，流れ場の渦運動と粒子運動とは，密接に関係していることが知られてい る．例えば，粒子の乱流拡散現象を数値解析によって再現する場合，工業上実用的で汎 用性のあるReynolds Averaged Navier-Stokes（RANS）型の手法では，乱れの統計的性質 に基づき，時空間スケール・等方性・確率密度関数等の仮定を導入することで，粒子に 及ぼす乱れの影響をモデル化する方法が用いられる．
また，流れ場に対してLarge Eddy Simulation（LES）や直接数値計算（Direct Numerical
Simulation, DNS）を用いて，乱流中の粒子運動を解析する研究が行われている．湯ら[6] -
[8]は，固気二相噴流を対象にして，乱流場と粒子との干渉を調べている．解析結果から，
流れ場の大規模渦やそれに付随する不規則運動によって，粒子がクラスタを形成するこ とを示している．さらに最近では，粒子周りの流れまで解像して解析する本来の意味で の混相流直接数値計算を行う研究例が報告されている．梶島ら[9], [10]は，静止流体中を落 下する単独粒子の挙動や，一様な流れ場中を多数の粒子が落下する現象について，乱流・
粒子運動の両方とも実験・理論式や物理モデルを用いない直接数値計算を実施している．
粒子を質点モデルとして扱う手法と，粒子周りの流れまで解析する手法とを比較するこ とで，質点モデルでは，粒子からの非定常渦放出の影響，粒子の回転とせん断による揚 力，の二つの問題が顕在化することを示している．ただし，現状ではLESやDNS を適 用できるのは乱流状態を維持する下限に近い低レイノルズ数で，比較的単純な形状の流 れ場で，解析領域も限られる．そのため，これらの解析法を実機スケールの工業装置へ 適用するには，計算コストの増大，経験定数の使用，解析領域外の影響を反映させるた めの複雑な境界条件，などの面で，現状では実用することが難しい．
以上説明した解析は，流体場にオイラー型の計算格子を用いる領域型の解析手法であ る．これに対して本研究では，流れ場を渦法で，粒子運動を粒子追跡法で，流体相と粒 子相ともラグランジュ的に解析する手法を用いることによって，流れ場に存在する流体 塊の渦運動と，その渦から誘起される流体力によって運動する粒子を時々刻々追跡する．
粒子運動を支配する様々なスケールを持つ乱流渦は，非定常で非線形運動する渦要素で 表現されるため，粒子に及ぼす乱れのモデル化を必要とせずに，高レイノルズ数で非定 常性の強い渦流れに対しても比較的簡易な解析が可能となる．このように，渦法と粒子 追跡法とのカップリング手法は，①RANS型乱流モデルが不要である，②渦構造の発達・
散逸過程を精度良く捉えられ，乱流渦と粒子との干渉を時々刻々追跡できる，③格子形 成が不要であり，境界の移動や変形を伴う問題を比較的簡易に扱える，④渦同士の非線 形運動により，自動的に高解像度な解析ができる，等の利点を持つ[11] - [14], [50], [51]．
渦法と粒子追跡法を組み合わせた混相流解析手法がもつ上記のような特徴を活かして，
渦運動と粒子運動との関係を調べる研究が多く行われている．二次元混合層，円形噴流，
物体後流などのせん断乱流中での粒子の乱流拡散解析において有用性が数多く報告され
ている[11] - [14]．また最近では，内山ら[15] - [17]，Joiaら[18] ，Waltherら[19]により，流体相と
粒子相との二相間相互作用を考慮したTwo-wayモデルも提案されている．流体相および 粒子相ともにラグランジュ解析を用いると，Two-way couplingの計算の際に離散点数の べき乗で演算量が増大する問題があるが，上記研究例ではオイラー的手法を援用するこ とで計算負荷を低減させる試みが行われている．さらに，内山らは，Two-wayモデルを 導入した三次元解析に発展させており，固気二相同軸円形噴流や自由落下粒子群が形成 する粒子噴流の解析を行っている[20], [21]．
このように，渦法と粒子追跡法を用いた混相流解析は，ラグランジュ解析が得意とす る混合層や噴流等の外部流れに適用されている．ただし，混相流を対象としたこれまで
の研究例[11] - [21]では，はく離せん断層位置が固定された解析がほとんどであり，このよ
うな解法では解析領域の全域において粘性が無視できないような内部流れを解析するこ とはできない．そのため，渦法と粒子追跡法による混相流解析を内部流れに適用する試 みはこれまで行われていなかった．しかしながら，管内輸送に代表される混相流装置の 多くが内部流れであり，壁面摩耗や堆積等が問題となる．このような工学上重要な諸問 題を数値解析によって予測するためには，流路の複雑な構造や三次元形状を離散的にモ デル化すると共に，流路壁面からのはく離や逆流を考慮した内部流れの解析が必要とな る．そこで本研究では，渦法と粒子追跡法によるラグランジュ・ラグランジュ型
（Lagrangian-Lagrangian）の解析法を固気・固液の内部流れ混相流解析に適用する．
1.2 流体－粒子の混相流解析に関する従来の研究
1.2.1 流体－粒子の混相流解析の概要とモデル化の視点
流れ場中に，固体粒子，気泡，液滴などの分散した粒子が混じった流体－粒子系混相 流に対する数値解析には，基礎となる数理モデルによっていくつかの種類がある．空間 スケールを基準にすると，図1.2.1[23]のように整理することができる[23], [24]．
一つ目に，流体相および分散相（粒子）共にミクロな視点に基づく解析は，個々の粒 子周りの流れと粒子個々の運動を解く手法であり，最も厳密で信頼性が高い．このよう に，流体相を微視的な視点で解くアプローチでは，粒子間隙の流れに対してNavier-Stokes 方程式が解かれ，粒子に働く流体力はモデル化無しに求められる．粒子がつくる流体乱 れの影響が強く，粒子個々の運動が問題となる流れでは，この視点からのアプローチが 有効となる．ただし，この手法は大きな計算負荷を要するため，取り扱える粒子数が限 られる．最近では，数千個の球形粒子の運動を扱った研究例が報告されている[9], [10]．
二つ目に，粒子の運動はミクロな視点（粒子追跡法）で，流体相に対しては粒子スケ ールよりも大きな空間スケールで平均化されたマクロな流れとして解く手法がある（オ イラー・ラグランジュ法，ラグランジュ・ラグランジュ法）．この手法では，粒子に働く
流体相
ミクロな視点
・粒子径より小さい解像度
・粒子に働く流体力にモデル化不用
マクロな視点
・粒子径より大きい解像度
・相間相互作用をモデル化
分散相（粒子）
ミクロな視点
・個々の粒子について運動を記述
・並進，回転，衝突など
マクロな視点
・粒子群を連続体としてモデル化
・相間相互作用をモデル化
分散相（粒子群） 流体相 モデル化無し
の直接解析
オイラー・
ラグランジュ
ラグランジュ・
ラグランジュ
オイラー・オイラー
流体力およびその反力についてはモデリングが必要となる．この視点に基づく数値解析 では，数十万個以上の粒子を取り扱うことも可能であり，流動層における気泡の形成や，
粒子クラスタのメゾスケール構造が解析で再現されることが報告されている[25], [26]． 三つ目に，流体相および分散相ともマクロな視点でモデル化する解析法では，上記の ような流体相に対するモデルに加えて，粒子群も一種の連続体として取り扱う連続体モ デルを用いる（オイラー・オイラー法）．この手法は，二流体モデル[27]とも呼ばれてい る．粒子個数による解析の制限が無く，計算負荷が小さい利点がある．しかしながら，
二流体モデルでは，粒子間の相互作用のモデル化が特に問題となる．これまでに，分散 相にニュートン粘性を仮定したモデルや，粒子に対する気体運動論に基づくモデルなど が提案されている．また，連続的な粒子径分布を考慮することが困難である事や，数値 拡散の影響が大きい事などの短所がある．
以上のように，流体－粒子の混相流解析手法には，空間スケールによって分類すると 三つの手法がある．ここで，緒言にて述べたように，本研究で対象とする固気・固液の 混相流において，特に問題となる摩耗や堆積などを数値解析で予測する場合には，固体
図1.2.1. 流体－粒子系混相流解析のモデル化の視点
粒子と壁面との衝突現象を精度良く解析することが必要となる．このような場合には，
個々の粒子について並進や回転の運動だけでなく，衝突運動を精度良く取り扱うことが できるラグランジュ手法（粒子追跡法）が有効となる．例えば，数値シミュレーション によって固体粒子の挙動をラグランジュ的に追跡することで，固体粒子と壁面との衝突 時の速度，角度などから，摩耗位置や摩耗量を推定する研究が行われている[28] - [31]．
そこで，粒子運動の解析にラグランジュ的なアプローチを行うオイラー・ラグランジ
ュ法（Eulerian-Lagrangian），および本研究で用いるラグランジュ・ラグランジュ法
（Lagrangian-Lagrangian）について，以降で詳しく解説する．
1.2.2 オイラー・ラグランジュ型の混相流解析
(1) オイラー型流れ解析による乱流の解析法
緒言にて述べたように，工学上問題となる流れの多くは高レイノルズ数の乱流であり，
分散相である粒子の挙動は，流れ場の様々なスケールの乱流渦に支配される[1] - [5]．そこ で，オイラー型流れ解析による乱流の解析法のおおまかな分類を図1.2.2[23]に示す．
DNSは，流れの基礎方程式であるNavier-Stokes方程式を乱流のモデル化無しに直接数 値計算する．エネルギの大半を有する大規模な渦から，エネルギ散逸を担う最小規模の 渦まで全て計算する．計算領域は最大の渦（流路や物体の寸法）よりも大きく，計算格 子は最小の渦（コルモゴロフスケール）よりも小さくしなければならない．したがって，
計算に必要な格子分割数は，少なくとも両者の比であり，例えばレイノルズ数 Re に対 してRe 3/4程度と見積もられている．乱流は本質的に三次元現象であることを考慮すれば，
少なくともRe 9/4程度の格子数を要することになる．これに対して，現在の計算機性能で
図1.2.2. オイラー型流れ解析による乱流の解析法
基礎研究レベル
RANS (Reynolds Averaged Navier-Stokes) simulation ほとんどの乱れをモデル化
LES (Large Eddy Simulation)
グリッドスケール以上の乱れは計算
DNS (Direct Numerical Simulation) 乱れエネルギの散逸まで計算
モデル依存大
計算負荷大 実用レベル
取り扱える格子数は 1003～10003程度が現実的である．以上の見積もりから明らかなよ うに，現状でDNSを適用できるのは，乱流状態を維持できる下限に近い低レイノルズ数 で，非常に単純な形状の流れ場に限られる．
一方，工業的にはレイノルズ平均Navier-Stokes方程式（RANS）を用いた解析が実用 的である．しかし，単相乱流でさえRANSに対する普遍的な乱流モデルは確立されてお らず，粒子による様々な効果のモデリングはさらに困難となる．本質的な問題として，
流れ場の乱流に対する平均操作と混相流に対する平均操作の関係を普遍的に定義できな い点がある．したがって，混相乱流に対するRANSでは，流れ場ごとに固有のモデルを 構築しなければならず，多くの経験的パラメータに依存せざるを得ない．
これに対して，DNSとRANSの間に位置する手法としてLESがある．LESの利点は，
平均化する空間スケールが定義されていることである．LESでは，流れの基礎方程式に 空間的なフィルタをかける．フィルタの幅としては，計算の格子幅程度が選ばれる．そ れよりも小さなスケール（Sub-Grid Scale，SGS）だけに乱流モデルが適用され，格子で とらえられるスケール（Grid Scale，GS）の変動のみが計算の対象となる．単相乱流に対 するSGSモデルが確立されているとは言い難いが，計算格子サイズより小さなスケール の乱流渦だけをモデル化すればよいので，RANS に比べると普遍性がある．ただし，混 相乱流に対するSGSモデリングは，まだ研究が始まったばかりである．
(2) オイラー・ラグランジュ法による流体－粒子の混相流解析法
上述のように，自然現象や工業装置内での混相乱流の解析には，現状ではRANSかLES を用いるのが現実的である．いずれにおいても，乱流中の粒子運動方程式および粒子に よる乱流変調のモデルが必要となるが，十分に確立されていない．
先ず，個々の粒子の運動方程式について述べる．球形の剛体粒子に関しては，
Basset-Boussinesq-Oseen formulation（BBO 方程式）を基礎とし，様々な効果を表現する
項の総和として質点モデル[32], [33]が提案されている．粒子の運動量変化は以下の式で表さ れる．
dt m d p up
＝抗力（定常力＋非定常力）＋揚力＋重力・浮力＋圧力勾配力
＋履歴力＋接触力＋その他の力
ここで，mpは粒子の質量，upは粒子重心の速度である．
質点モデルは，比較的静かで，広い流体中を運動する単独の粒子についてはよく研究 されている．しかし，非球形粒子，高濃度粒子群，壁近傍，乱流中での挙動を表現でき るかについては，ほとんど検証されていないのが現状である．また，粒子が高濃度で存 在する場合には，接触力（運動量・角運動量などの変化）のモデル化が必要となる．こ れについては，バネとダッシュポットなどを用いたモデルが提案されている[22] - [24]．
次に，粒子による乱流変調モデルについて述べる．乱流変調は，粒子の混入による乱
れの増減を意味することが多い．これに関してはいくつかのパラメータが提案されてい
る[34], [35]．複雑な乱流モデルに粒子効果を取り入れるためには，粒子による乱れの増減だ
けでなく，非等方性やエネルギ注入スケールが問題となるが，普遍的なモデルは見いだ されていない．LESによって乱流変動が解析される場合には，粒子に作用する流体力の 反力が流れの運動方程式に返されるが，SGS乱れに対する効果はモデルで与えられなけ ればならない．したがって，粒子の乱流への影響が妥当に与えられる空間スケールには 下限がある．
以上のように，オイラー・ラグランジュ法による数値解析を用いて，混相乱流場を高 精度で再現するためには，上述したモデル化や計算誤差を極力削減することが必要であ る．したがって，その手法はDNSあるいはLESであることが望まれる．ただし，この ような混相流計算を実機スケールの工業装置へ適用するには計算コストが増大するため，
現状では実用化することが困難である．
1.2.3 渦法と粒子追跡法によるラグランジュ・ラグランジュ型の混相流解析
本研究で用いるラグランジュ・ラグランジュ型の混相流解析の特徴を明確に示すため，
オイラー・ラグランジュ法とラグランジュ・ラグランジュ法との解析方法の違いを図
1.2.3に示す．前項で記したオイラー・ラグランジュ型の混相流解析は，図1.2.3 (a)に示
すように，流れ場の計算にオイラー的な計算格子を用いる手法である．上述したように，
オイラー・ラグランジュ法では以下のような課題が挙げられる（一部はDNSを除く）．
・時間や空間によって平均化された乱流モデルが必要となる．
・流体－粒子間での乱流モデルが未だ研究段階である．
図1.2.3. オイラー・ラグランジュ法とラグランジュ・ラグランジュ法
Vortex Particle
Vortex
Particle Computational grid
(a) オイラー・ラグランジュ法 Fluid
flow
(b) ラグランジュ・ラグランジュ法
・計算格子内で流体相の速度を平均化あるいは内挿する必要がある．
・流れを予め推測し，必要とされる空間解像度と計算負荷に応じて，最適な計算格子 サイズの粗密を決定する必要がある．
これに対して本研究では，流れ場を渦法で，粒子運動を粒子追跡法で解析する．図1.2.3 (b)に示すように，流体相および粒子相ともにオイラー的な計算格子を必要としないラグ ランジュ・ラグランジュ型の混相流解析手法である．主な利点を以下に挙げる．
・乱流モデルなしに，高レイノルズ数の非定常渦流れを解析できる．
・渦構造の発達・散逸過程を精度良く捉えられ，乱流渦と粒子との干渉を時々刻々追 跡できる．
・格子形成が不要であり，境界の移動や変形を伴う問題を比較的簡易に扱える．
・渦同士の非線形運動により，自動的に高解像度な解析ができる．
・粒子位置（粒子の重心）での流体相の速度や速度勾配を直接算出できる．
このようなラグランジュ・ラグランジュ型混相流解析の特徴を活かして，渦運動と粒 子運動との関係を調べる研究が多く行われており，ラグランジュ解析が得意とする外部 流れに適用されている[11] - [21]．ただし，緒言にて述べたように，混相流を対象とした従 来の研究では，内部流れに適用する試みは行われていない．そこで本研究では，混相内 部流を対象とした数値解析法を提案する．本研究の目的を次節に記す．
1.3 本研究の目的
緒言および従来の研究例にて述べたように，本研究で用いる渦法と粒子追跡法による ラグランジュ・ラグランジュ型の混相流解析法を，実用産業機器にみられる多種多様な 混相流場に対して有効に利用するためには，複雑形状や三次元構造を考慮した内部流れ に適用する必要がある．一般に，二次元を対象とした外部流れとして円柱・平板等の比 較的単純な形状周りの流れ解析については，等角写像を用いて境界をモデル化すること が可能である．しかし，内部流れや三次元複雑形状に対しては，境界を離散的にモデル 化することによって，数値的に解析する必要がある．また，壁面からのはく離や逆流が 発生する内部流れを解析するためには，解析領域を囲む全ての固体壁において，粘性に より生成される渦度を離散化して導入する必要がある．
これらの問題を解決するため，本研究では混相内部流を対象とした境界要素法と壁面 渦導入法を組み合わせた二次元および三次元の解析法を提案する．また，解析法の妥当 性を検証するため，既存文献による比較だけでなく，検証実験を行うことで，様々な混 相内部流場に対して，本解析の精度と適用範囲を定量的に明らかにする．さらに，相間 相互作用を考慮した内部流れ混相流解析まで発展させると共に，実際に利用されている
工業装置を対象とした解析を試みることで，一般産業機器に応用可能な数値解析技術と してまとめる．
1.4 本論文の概要
本研究は，輸送や攪拌などを伴う多くの工業装置内で観察される固気・固液混相流の 数値解析に関するものである．流れ場の解析には乱流モデルや格子生成を必要としない 非定常流れ解析法として応用の進んでいる渦法を適用し，粒子運動の解析には個々の粒 子について並進・回転・衝突の運動を追跡する粒子追跡法を適用することによって，新 たに開発した工業的に汎用性の高いラグランジュ・ラグランジュ型の内部流れ混相流解 析法について論じた論文であり，5章から構成される．
第1章では，内部流れを対象にした渦法と粒子追跡法による混相流解析の研究を行う 背景と目的について述べる．流れ場中に固体粒子，気泡，液滴などの分散した粒子が混 じった流体－粒子系混相流の数値解析法として，直接解析法，オイラー・ラグランジュ 法，オイラー・オイラー法，そして本研究で開発するラグランジュ・ラグランジュ法の それぞれについて，モデル化の概要や特徴などを述べる．
一般に，工学上問題となる流れの多くは高レイノルズ数の乱流であり，分散相である 粒子の挙動は，流れ場に存在する様々なスケールの乱流渦に支配されるため，乱流渦と 粒子との干渉を精度良く解析することが重要となる．渦法と粒子追跡法による混相流解 析法は，乱流モデルなしに高レイノルズ数の非定常流れを解析でき，乱流渦と粒子との 干渉を時々刻々追跡できる特徴を持つ．この特徴を活かし，ラグランジュ解析が得意と する混合層や噴流等の外部流れに適用されている．ただし，従来の研究では，はく離せ ん断層位置が固定された解析がほとんどであり，内部流れに適用する試みはこれまで行 われていなかった．しかしながら，粉体やスラリーの管内輸送に代表される混相流装置 の多くが内部流れであり，粒子による壁面摩耗や堆積等が問題となる．そこで本研究で は，上記のような工学上重要な諸問題を数値解析によって予測するため，産業機器に多 くみられるような複雑な流路構造や三次元形状を有する混相内部流に適用できるラグラ ンジュ・ラグランジュ型の混相流解析法を開発することを目的とした．
第2章では，渦法と粒子追跡法によるラグランジュ・ラグランジュ型の内部流れ混相 流解析の手法と数値モデルを解説し，内部流れ解析へ拡張するために新たに提案した解 析手法の詳細を述べる．渦法と粒子追跡法による混相流解析は，流体相および粒子相と もラグランジュ的に解析する手法であり，流れ場に存在する流体塊の渦運動と，その渦 から誘起される流体力によって運動する粒子を時々刻々追跡する．粒子運動を支配する 様々なスケールの乱流渦は，非定常で非線形に運動する渦要素で表現されるため，粒子 に及ぼす乱れのモデル化を必要とせずに，高レイノルズ数で非定常性の強い渦流れに対 しても比較的簡易で計算負荷の小さい解析が可能となる．本研究で開発した渦法と粒子 追跡法による内部流れ混相流解析において主体となる数値計算方法には，内部流れの計
算方法，粒子運動の計算方法，相間相互作用の計算方法の三つがある．
第一に，外部流れの非定常解析において有用性が数多く報告されている渦法による流 体解析法を内部流れ解析へ拡張するため，内部流れ場を境界要素法によって解析する手 法と壁面渦導入法とを組み合わせた二次元および三次元の解析手法を提案した．境界要 素法を用いた内部流れ解析では，流路内壁および流入出面などの境界を，ポテンシャル 値 およびポテンシャル流束値 nの両方を未知数として離散化し，境界積分方程式 を直接法によって定式化する方法を提案した．流入面には速度流入境界条件，流出面に は自由流出境界条件および連続の条件，壁面には法線方向速度がゼロの境界条件を与え，
および nを解いた．これに加えて，粘性による壁面上のNo-slip条件を満足させる ため，壁面上に発生する渦度を渦要素で全て置き換え，それらの挙動をラグランジュ的 に追跡する手法を用いることで，流路内で発生するはく離や逆流などを捉えられるよう にした．ここで，流路壁近傍では速度勾配が非常に大きくなるため，渦要素の急激な伸 縮が起こる．これに対しては，壁面を離散化したパネルサイズに応じて，導入渦要素を 多数の渦要素に予め分割して導入する方法を提案し，計算精度の低下を抑制した．また，
壁面近傍に限らず，渦要素の伸びが激しい場合に対しては，対象渦要素を多数の小球形 要素に再分割する方法を用い，渦度の発散と計算精度の低下を回避した．
第二に，粒子運動の解析には粒子追跡法を用いて，個々の粒子の並進・回転・衝突の 三つの運動を考慮してラグランジュ的に追跡した．粒子の並進運動に作用する外力は，
流体抗力（定常力および非定常力），Magnus揚力，Saffman揚力，重力・浮力を考慮し，
三次元の固体球に対する実験式や理論式から算出した．粒子の回転運動は，流体粘性 および衝突による回転を計算した．粒子の衝突運動に関しては，剛体球粒子を仮定して 摩擦のある非弾性衝突モデルを用いた．本研究では，混相流装置内で発生する摩耗や堆 積などの現象を精度良く予測するため，決定論的手法によって衝突判定を行い，個々の 粒子に対して衝撃運動方程式を解いた．
第三に，流体と粒子の間の相間相互作用の解析法には，渦法と粒子追跡法による外部 流れを対象とした混相流解析で有効性が報告されているオイラー的手法を援用した
Two-wayモデルを用いた．粒子の混入によって生じる流れ場の変化を，個々の粒子が受
ける流体力の反作用としてモデル化することによって，流れ場解析における渦要素の渦 度変化を時々刻々計算した．
第3章では，第2章で提案した渦法と粒子追跡法を用いた混相流解析法を二次元およ び三次元の内部流れに適用し，解析手法の妥当性，精度や適用範囲を検証した結果を述 べる．二次元解析では，固液二相噴流，混合層を伴うチャネル内での固液二相流，鉛直 チャネル内での固液二相流，T 字型合流路内での固液二相流を対象とした．三次元解析 では，一様流中への固液合流噴流，固気二相噴流，T 字型合流路内での固液二相流を対 象とした．また，本解析において計算負荷が高いBiot-Savart法による誘起速度の算出に 対して，高速化手法を提案し，計算負荷と精度の最適化を図った．
以上のような様々な混相内部流を対象に解析を行い，流体相および粒子相に関する時
間平均および時間変動特性などについて，既存文献による比較検証だけでなく，水とガ ラスビーズを用いた固液二相流の検証実験を実施し，粒子フローパターンの可視化，
Laser Doppler Velocimetry（LDV）やParticle Image Velocimetry （PIV）を用いた速度計測 を行うことによって，解析手法の精度と適用範囲を定量的に明らかにした（検証した条 件範囲は，Re=5000~42000, St=0.05~50）．また，他手法（オイラー・ラグランジュ法お よびオイラー・オイラー法）と本手法との比較を行うことで，本研究で用いたラグラン ジュ・ラグランジュ法の優位性を示した．本研究で開発した内部流れ混相流解析法を用 いることによって，高レイノルズ数で非定常性の強い渦流れによる粒子の乱流混合，摩 耗，堆積などの現象に対して，粒子のフローパターンだけでなく，壁面衝突や堆積の有 無およびそれらの発生位置等を予測し，装置の設計・開発に有益なデータが得られるこ とを示した．
さらに，流体と粒子の相間相互作用を考慮した内部流れ混相流解析として，粒子群の 自由沈降流や，堆積を伴う高濃度粒子流を対象にTwo-way解析を行った．実験との比較 から本解析の妥当性を確認し，低濃度だけでなく，相間相互作用の影響が無視できない 中・高濃度の粒子を含む混相流場の内部流解析が可能であることを示した．
第4章では，工学的に利用されている実用流体機器への応用について述べる．氷蓄熱 装置における熱交換器内での氷水流れ，および石炭焚き発電用ボイラにおける伝熱管群 周りの灰粒子流れを対象に，本解析技術を応用した．非定常な渦流れを伴う実際の装置 での混相内部流に対して本解析の有用性を示すだけでなく，装置内で発生する摩耗や堆 積など混相流特有の現象を解明し，機器の改良に資する適切な知見が得られた．このよ うに，解析技術の開発や検証に止まらず，工学的な実用問題を取り扱うことで，一般産 業機器に応用可能な技術にまで進展させた．
第5章は本研究をまとめるものであり，本研究で提案した渦法と粒子追跡法によるラ グランジュ・ラグランジュ型の内部流れ混相流解析法の有効性について総括し，今後の 展望を述べる．
最後に，本論文を構成する公表論文を本論文末尾に挙げる．
第 2 章 渦法と粒子追跡法による混相流解析手法と数値モデル
2.1 はじめに
本章では，渦法と粒子追跡法を用いたラグランジュ・ラグランジュ型の内部流れ混相 流解析の手法と数値モデルについて述べる．特に，渦法と粒子追跡法による混相流解析 において，これまで適用されていなかった内部流れへ適用することを目的に，本研究で 新たに提案する数値解析法について詳しく述べる．主体となる数値解析法には，境界要 素法と渦法による内部流れの計算方法，粒子追跡法による粒子運動の計算方法，相間相 互作用の計算方法の三つがある．本章では，渦法と粒子追跡法による混相流解析の特徴 や解析のアルゴリズムから解説し，上記した三つの手法についてそれぞれ詳しく述べる．
2.2 渦法と粒子追跡法によるラグランジュ・ラグランジュ型の混相流解析 の概要
本研究では，流れ場を渦法で，粒子運動を粒子追跡法で解析する．解析手法の概略を
図2.2.1に示す．流体相および粒子相ともラグランジュ的に解析する手法を用いることに
よって，流れ場に存在する流体塊の渦運動と，その渦から誘起される流体力によって運 動する粒子を時々刻々追跡する．粒子運動を支配する乱流渦は非定常で非線形に運動す る渦要素で表現されるため，粒子に及ぼす乱れのモデル化を必要とせずに，高レイノル
completely grid-free multiphase flow
(fluid and solid particles) wall
inlet
outlet
solid particles vortex elements (fluid)
図2.2.1. 渦法と粒子追跡法によるラグランジュ・ラグランジュ型の混相流解析法
ズ数で非定常性の強い渦流れに対しても比較的簡易な解析が可能となる．
以降に，渦法と粒子追跡法を用いたラグランジュ・ラグランジュ型解析の特徴を記す．
2.2.1 解析手法の特徴
緒言にて述べたように，本研究で用いる渦法と粒子追跡法とのカップリング手法は，
①RANS 型乱流モデルが不要である，②渦構造の発達・散逸過程を精度良く捉えられ，
乱流渦と粒子との干渉を時々刻々追跡できる，③格子形成が不要であり，境界の移動や 変形を伴う問題を比較的簡易に扱える，④渦同士の非線形運動により，自動的に高解像 度な解析ができる，等の利点を持つ[11] - [21], [50], [51]．また，流体相および粒子相ともオイ ラー的なグリッドを用いないため，グリッドベースで規定される格子解像度や時間刻み 等の制限がなく，解析精度と計算負荷を比較的自由に決定できる．そのため，パラメー タスタディや試行錯誤が伴う設計・開発に有用な手法である．ただし，渦要素数の増加 に伴って計算時間が増大すること，三次元流れでは渦度（渦要素）の急激な伸縮に対応 する工夫が必要なこと，圧縮性流れへの拡張が難しいこと，などの課題も挙げられてお り，最近ではこれらの解決に向けた研究も行われている．
上記のような特徴を活かして，渦運動と粒子運動との関係を調べる研究が多く行われ ており，せん断乱流中における粒子の乱流拡散解析への有用性が報告されている[11] - [14]． また最近では，内山ら[15] - [17], [20], [21] ，Joiaら[18]やWaltherら[19]により，二相間相互作用 を考慮したTwo-wayモデルも提案されている．流体相および粒子相にラグランジュ解析 を用いるとTwo-way couplingの計算の際に，離散点数のべき乗で演算量が増大する問題 があるが，上記研究例ではオイラー的手法を援用することで計算負荷を低減させる試み が行われている．このように，渦法を用いた混相流解析は，ラグランジュ解析が得意と する混合層や噴流等の外部流れに適用されている．しかし，これまでの研究例では，内 部流れ混相流に適用されていない．
本研究では，混相内部流を対象とし，境界要素法と壁面渦導入法を組み合わせた手法 を提案する．境界要素法を用いた内部流解析では，ポテンシャル値 およびポテンシャ ル流束値 nの両方を未知数として境界を離散化し，境界積分方程式を直接法によっ て定式化する．加えて，粘性による壁面上のNo-slip条件を満足させるため，壁面上に発 生する渦度を渦要素で全て置き換え，それらの挙動をラグランジュ的に追跡する．
2.2.2 解析のアルゴリズム
本研究で提案するラグランジュ・ラグランジュ型の内部流れ混相流解析のアルゴリズ ムの概略を図2.2.2に示す．内部流れの解析を行うため，境界要素法（直接法）と渦法を 用いる．ポテンシャル値 およびポテンシャル流束値 nの両方で境界をモデル化し，
流入出・壁面などの境界条件を与え， および nを解く．加えて，粘性による壁面
上のNo-slip条件を満足させるため，壁面上から渦要素を導入し，渦度の変化や移動をラ グランジュ的に追跡する．粒子運動は粒子追跡法を用い，並進・回転・衝突を考慮し，
個々の粒子について運動方程式を解く．また，流体と粒子との相互作用を考慮し，粒子 による渦度の変化を計算する．解析のアルゴリズムは，以下のような計算手順で構成さ れる．
① 解析対象のモデル化，解析の開始
対象とする混相流場に対して，解析領域を設定する．一般に，計算時間が制限 されるため，解析領域を限定する必要がある．また，解析すべき物理現象を分析 して，必要十分でかつ効率の良いモデル化を行う．以上から，設定した解析領域 の境界形状（流入出，壁面など）をパネル状に離散化する．
② 境界要素法による解析
離散化した境界パネル上で，境界条件（流入出，壁面など）を与える．境界条 件を満足するように，境界要素法を用いて境界積分方程式を解き（直接法または 間接法），解析領域内のポテンシャル流れ（渦無し流れの流速upot）を求める．
上記の①，②のステップにて，ポテンシャル流れの解析結果が得られる．次に，粘性 による渦度のある流れと粒子運動（分散相）をカップリングして，非定常解析を行う．
③ 渦要素の導入
壁面上においてNo-slip条件を満足するように渦要素を流れ場に導入する．また，
流入境界面において速度不連続面がある場合には，速度勾配に相当する渦要素を 流れ場に導入する．ここで，流路壁近傍では速度勾配が非常に大きくなるため，
壁面を離散化したパネルサイズに応じて，導入渦要素を多数の渦要素に予め分割 して導入する．
④ 粒子（分散相）の導入
分散相である粒子に対して，物性・濃度・速度などの解析条件に基づき，粒子 を流れ場に導入する．
⑤ 粒子の運動
粒子に対して，並進・回転・衝突の運動方程式を解き，粒子を移動させる．こ のとき，粒子重心位置での流体相の速度や速度勾配を算出する必要があるため，
Biot-Savartの法則を用いる．
⑥ 渦度の変化
渦要素の渦度変化は，渦度輸送方程式に基づき計算する．移流による伸縮およ び粘性による拡散から，渦度の変化を計算する．また，流体相と粒子相との相互 作用による渦度の変化をTwo-wayモデルにより計算する．以上のような渦度の変 化によって渦要素が変形するため，必要に応じて複数の渦要素に分割し，再配置
を行う．
⑦ 渦要素の速度計算と移動
渦要素を代表する位置（渦要素の中心点）における速度をBiot-Savartの式より 算出し，その時間積分量から渦要素を移動する．
⑧ 境界要素法による解析
渦要素（渦流れ）とポテンシャル流れとを重ね合わせた粘性流れに対して，再 度，境界条件（流入出，壁面など）を与える．境界条件を満足するように，境界 要素法を用いて境界積分方程式を解き，解析領域内のポテンシャル流れを求める．
これにより，粘性流れの流速uは，渦無し流れの流速upotと，渦度が存在する流 れの流速uvorの重ね合わせで表される．
⑨ 解析領域内の速度場・圧力場の算出
必要に応じて，現在の計算時刻における速度場・圧力場を算出する．圧力場は，
同時刻における速度場によって決定される圧力 Poisson 方程式よって，従属的に 求めることができる．
上記の③～⑨のステップにて，粘性による渦度のある流れ，および分散相の粒子運動 をカップリングした非定常解析を行うことができる．必要とする解析時間に到るまで，
③～⑨のステップを繰り返す．
⑩ 解析終了
上記ステップで得られた解析結果から，流体相および粒子相に対する時間平均 値・変動値などを算出し，必要とする情報を抽出して評価・考察する．
粒子と渦要素の干渉
境界要素法による解析 渦要素の導入
渦度の変化
渦要素の速度計算・移動
解析領域内の速度場・圧力場の算出 境界条件：流入，流出，壁面等 壁面上のNo-slip条件
渦度の伸縮，拡散，再配置
Biot-Savart，ラグランジュ移動
Biot-Savart，圧力ポアソン 渦法
粒子（分散相）の導入
粒子の運動
物性，濃度，速度
並進，回転，衝突
Two-wayカップリング
粒子追跡法 境界要素法による解析
境界条件：流入，流出，壁面等 解析開始
対象場（境界形状）の離散化
解析終了
解析結果の評価・考察
図2.2.2. 解析のアルゴリズム
繰り返し
2.3 渦法による流れ解析
2.3.1 渦法による流れ解析の概要
渦法は，速度場の計算方法によりBiot-Savart法とVortex in Cell （VIC） 法の二つの 手法に大きく分けられる．渦要素の持つ渦度の変化は粘性・非粘性，あるいは層流・乱 流に関わらず，一般的に非圧縮性流れに対して成り立つ渦度輸送方程式によって求めら れる．また，速度場については，VIC法では差分法と同様に流れ場を計算格子で分割し，
渦度分布をソース項とするPoisson方程式から流れ関数を求める．各格子点での速度は，
流れ関数の微分係数を差分により近似することで求める．一方，Biot-Savart 法において は，渦度の定義式から積分形として導かれる一般化されたBiot-Savartの式を用い，各要 素からの誘起速度の和として速度場が計算される．このように，渦度輸送方程式および 速度場の計算において，渦度を媒介変数として考えるならば，未知数は流れ関数あるい は速度である．これらの方程式が与えられた境界条件のもとに解くことができるならば，
流れ場が速度もしくは渦度によって規定されることを意味し，圧力は直接的な関わりを 持たない．つまり非圧縮粘性流れにおいて，流れ場は渦度の生成，輸送および消散の過 程によって決定され，圧力は従属的な量と見なすことができる．これは，渦法を理解す る上で重要な事柄である．このように，渦法は渦度輸送方程式および速度場の計算によ って流れを解析する計算法であり，流れ場の圧力計算を本質的に必要としないという特 徴を持つ．
渦法においては，流れ場に存在する渦度分布は，渦要素の持つ渦度分布の線形和によ って離散化される．渦要素のもつ渦度は，渦度輸送方程式に従い時間的に変化する．渦 度輸送方程式において，渦度は渦度ベクトル方向の伸縮を伴う移流および粘性による拡 散によって変化する．粘性拡散の考慮法としてChorin[36]は，渦度の拡散方程式の解が確 率的には分子のブラウン運動で表現できることを利用し，渦要素の移流に乱数による
Random walkを加える手法を提案した．またLeonard[37]は，時間の経過と共に渦要素モデ
ルの粘性核半径を拡大することにより渦度の減衰を表現する Core-spreading 法を提案し た．大上ら[38]は，渦度輸送方程式の粘性拡散項から速度の次元をもつ拡散速度を定義し，
渦要素同士の相互作用として粘性拡散効果を考慮する手法を提案した．一方，三次元流 れ場においては，粘性効果に加え渦度ベクトル方向に伸縮を伴う移流の効果も考慮する 必要があり，渦糸，Stick，VortonおよびBlobなどの渦要素モデルを用いて，それぞれ異 なる方法で考慮されている．また，三次元の高レイノルズ数乱流への応用を考慮した渦 法における乱流モデルの開発も近年積極的に行われている．微小渦やサブグリッドスケ ールの渦の効果をモデル化することによるLESモデルに関しては，Leonard & Chua[40]は
Smagorinskyの渦粘性係数の概念に基づく渦核成長モデルを提案した．伊澤・木谷[41]は，
Leonard & Chua[40]の渦核成長モデルをやや変形した簡便な乱流モデルを提案し，円形ノ
ズルからの非定常噴流解析に適用した結果を報告している．Mansfieldら[42], [43]は，単一 渦輪の挙動および二つの渦輪の衝突に対してLESモデルを導入して解析した結果を報告
している．
現在のところ，三次元流れへの適用性および拡張性の観点から，Blob モデルと
Core-spreading法を組み合わせたアルゴリズムが最もシンプルであり，実用的であると考
えられる．Core-spreading法を三次元流れへ適用するため，移流による伸縮と粘性拡散に よる粘性核半径の拡大を同時に考慮した三次元 Core-spreading 法が中西・亀本ら[39]によ り提案されており，本研究ではこの手法を適用することにし，二次元および三次元流れ
におけるCore-spreading法についてそれぞれ詳しく説明する（2.3.5項，2.3.6項）．
渦法で非常に重要とされるのが，渦度の発生に伴う渦要素の流れ場への導入法である．
渦度は，流体粘性の作用によって形成される速度せん断層に分布することから，固体壁 面上の境界層や噴流が形成する自由せん断層などは，離散渦要素を分布させることによ って表現することができる．渦法においては，これら速度せん断層をどの程度のスケー ルの渦要素を用いて表現するかによって解析精度を決定することになる．例えば，微細 な渦要素を多数導入すれば，それに応じて高精度な解析を行うことが可能となる．また，
はく離せん断層のみを離散渦要素の分布で置き換えるような大胆な近似的扱いをするな らば，大規模な渦構造が支配的な非定常流れなどは計算負荷を大幅に軽減して簡略的に 解析することができる．このように，渦法は求める解の解像度に応じてフレキシブルに 適用することが可能である．亀本ら[44] - [46], [50]は，物体壁面の極近傍の速度勾配を直線勾 配で近似することによって渦層を仮定し，渦層内の渦度の移流および拡散を計算するこ とによって，渦要素を流れ場へ導入する手法を提案した．また，小島・亀本ら[52] - [54]は，
ブラフボディ，風車，トラクター・トレーラー周りの流れ等の三次元の非定常外部流れ に対して，亀本らが提案した渦要素導入法を適用し有効性を示した．そこで本研究では，
渦法による流れ解析を内部流へ適用するため，上記の亀本らにより提案された壁面から の渦要素導入法を用いた（2.3.5項）．
圧力に関しては，圧力 Poisson 方程式に時間項が含まれていないことから分かるよう に，参照する時刻以前の速度場・圧力場の状態とは無関係に，その時刻における速度場 から従属的に求めることができる．渦法においては，任意の位置における速度が
Biot-Savartの式より得られることから，圧力Poisson方程式から陽的に求めることができ
る．一般に圧力 Poisson 方程式を解くためには，空間に分布する未知数を反復法によっ て求める方法が用いられるが，計算格子が必要となるため渦法に適しているとは言えな い．そこで中西ら[47]は，圧力Poisson方程式に含まれる速度の空間的な勾配をBiot-Savart の式を空間微分することにより解析的に求め，境界要素法を用いて圧力 Poisson 方程式 を解く手法を提案した．ただし中西らの手法においても，差分法で用いられるような領 域分割型の格子ではないまでも，空間に速度参照点を設ける必要がある．これに対して，
Uhlman[49]は，圧力に関するPoisson方程式を積分方程式で定式化することによって，空
間に計算格子を設けることなく圧力を解析する手法を提案した．2.3.7項では，圧力計算 の基礎式として，速度場と圧力場の関係式である圧力に関する Poisson 方程式を示し，
圧力に関する積分方程式の導出方法について説明する．
2.3.2 基礎方程式
本研究では，非圧縮粘性流れを対象にする．非圧縮性流れでは密度 は一定であり，
基礎方程式は，次に示すような非圧縮 Navier-Stokes方程式と連続の式である．
u F
u
u u 2
1grad
grad p ν
t (2.3.1)
0
divu (2.3.2)
ただし，u，p， ， ，t，Fはそれぞれ，速度ベクトル，圧力，動粘性係数，密度，時 間，外力を表す．ここで，外力 Fを F = 0とすれば，式（2.3.1）は以下のようになる．
u u u 2u
1grad
grad p ν
t (2.3.3)
上式のrot（回転）をとれば，rot grad は常にゼロであるので圧力項は消え，次式のよ うに渦度ωと速度 uに関する式となる．
ω ω
ω u 2
t rot (2.3.4)
ここで，渦度ωは速度 uによって次式のように定義される u
ω rot (2.3.5)
さらに，式（2.3.4）を変形し，div u = 0 （連続の式）及び div ω=div rot u = 0 を代入 すれば，以下のような渦度輸送方程式（Vorticity transport equation）が得られる．
ω u
ω ω
ω u 2
grad
t grad (2.3.6)
この式において，右辺第一項は伸縮項（ストレッチ項，Stretch term）とも呼ばれ，渦 度の移流による渦度の回転軸方向の変化，及びその方向での渦度の大きさの変化を表す．
第二項は粘性による渦度の拡散を表す拡散項（Diffusion term）である．このように，非 圧縮性流体に関しては，流れ場の時間変化は渦度 ω と速度 u を未知数とする渦度輸送 方程式によって決定される．ここで，
grad t u
Dt
D (2.3.7)
の関係より，式（2.3.6）は以下のように表記できる．
ω u
ω ω 2
Dt grad
D (2.3.8)
特に，二次元流れ場の場合においては渦度の移流を表す右辺第一項は常にゼロになる ので，渦度輸送方程式は以下のようになる．
ω 2ω Dt
D (2.3.9)
したがって，二次元場を対象とした渦度輸送方程式は，流れ場中で渦度の存在する空 間を十分微小な体積に分割し，微小体積内の渦度の拡散を表す拡散項だけの上式を計算 することによって解くことができる．
一方，圧力pは，非圧縮 Navier-Stokes方程式の発散をとり，連続の式と任意ベクトル の回転の発散が恒等的にゼロであることを考慮して得られる次式の圧力ポアソン方程
式（Pressure Poisson equation）から求めることができる．
u u grad
2p div (2.3.10)
ここで，上式に時間項が含まれていないことからわかるように，参照する時刻以前の 速度場・圧力場の状態とは無関係に，その時刻における速度場から従属的に圧力場を求 めることができる．このように，渦法による流れ解析では，圧力場の計算を本質的に必 要としない特徴をもつ．
以上のように，渦法は渦度 ω と速度 u について，渦度輸送方程式と渦度の定義式
（ω rot u）を基礎式として流れを解析する数値計算法である．
2.3.3 Biot-Savart法による速度場の計算方法
流速uをベクトルポテンシャルAによって次のように表す．
A
u rot (2.3.11)
上式を渦度の定義式に代入すると，次式となる．
A A
u ω
div 2
grad rot rot rot
(2.3.12)
ベクトルポテンシャルAには任意のスカラ関数の勾配を付け加えてもよい性質がある
ため，div A = 0 の条件（Coulombゲージの条件式）を付け加えれば次式が得られる．
ω
2A (2.3.13)
上式を，ベクトルグリーン関数を用いて積分し，uについて解くと次式を得る[55]． ds
G G
dv
G S i i i i i i
Vωi i n u n u
u (2.3.14)
ここで，添字 iは離散化した多数の渦要素や境界の位置riでの変数を表し，niはriを 含む境界面上における法線方向単位ベクトルである．Gはスカラー・ラプラス方程式 の基本解である．Gはディラックのデルタ関数 を用いて，
2G 0 (2.3.15)
であり，二次元および三次元流れに対して以下のようになる．
G R1
2 log
1 ：二次元 (2.3.16)
G R 4
1 ：三次元 (2.3.17)
ここで，Rは任意の位置rに対して，R r ri である．
式（2.3.14）の右辺第一項の空間積分は，流れ場に存在する渦度から誘起される速度を
示す．第二項の境界積分は，境界上の速度で決定される渦度のないポテンシャル流れの
速度を示している（ポテンシャル流による誘起速度の計算方法の詳細は2.3.6項に記す）． このことから第二項を速度ポテンシャルの勾配 で置き換えると，次式を得る．
dv
Vωi iG
u (2.3.18)
上式はBiot-Savartの法則として知られる式[55]であり，流れ場に存在する渦度と速度場
の関係を表す．この式により，流れ場に存在する渦度の分布，および境界から誘起され る速度ポテンシャルの勾配が既知であれば，任意の位置における速度も一意的に決定す ることができる．
2.3.4 境界要素法による内部流れの定式化（間接法および直接法）
(1) 境界要素法の概要
非圧縮粘性流れの流速uは，渦無し流れ（Irrotational flow）の流速upotと，渦（渦度）
が存在する流れ（Rotational flow）の流速uvorの重ね合わせで表される．ここで，渦無し 流れ（ポテンシャル流れ）の解析に関して，円柱・平板等の比較的簡単な形状周りの二 次元を対象とした外部流れに限っては，等角写像を用い，理論的に解析することが可能 である．しかし，工業上利用される一般的な流体機器を想定した複雑形状や三次元構造 に対応するためには，解析技術の三次元拡張が必要不可欠となる．そのため，三次元複 雑形状や内部流れに対しては，境界や解析領域を離散化して，数値的に解析を行うこと が必要となる．そこで本研究では，三次元の境界要素法と渦法とを組み合わせ，三次元 の内部流れへ適用できる解析手法を開発する．
境界要素法（Boundary Element Method，BEM）では，解析対象の領域内部および境界 上での未知数を表す微分方程式（ポテンシャル流れ解析ではラプラス方程式）を，境界 値のみを関係付ける積分方程式に変換し，境界条件（ディリクレ，ノイマン）を満たす ように境界に分布させた未知数を数値解析的に解く．本研究で参考にした境界要素法に 関する書籍[57] - [62]および論文[63] - [70]を参考文献一覧に�
