博 士 論 文
エキゾーストマニホールド内の 非定常流れに関する研究
2020 年 3 月
群馬大学大学院工学研究科 先端生産システム工学領域
篠﨑 勇二
エキゾーストマニホールド内の非定常流れの研究
第1章 序論
1.1 本研究の背景 1
1.2 従来の研究
1.2.1 車体周りの流れの研究 4
1.2.2 エンジン筒内流れの研究 4
1.2.3 吸排気系脈動流の研究 5
1.3 本研究の目的 6
第2章 CFD理論と流速計測法について 2.1 緒言 9 2.2 CFD理論 2.2.1 支配方程式 9 2.2.2 エネルギ方程式 10
2.2.3 運動方程式 11 2.2.4 質量保存式 11 2.2.5 一般形方程式 11 2.2.6 メッシュの離散化 12 2.2.7 非定常一次元熱伝導の支配方程式と離散化 12 2.2.8 流れおよび熱伝導の数値計算法 16
2.2.9 時間平均方程式の誘導 16
2.2.10 乱流モデル 18 2.2.11 RANS 18
2.2.12 k-εモデル 19 2.2.13 k-ωモデル 19 2.2.14 LES 19 2.3 流速計測手法 2.3.1 熱線流速計 2.3.1.1 熱線流速計の特徴 21 2.3.1.2 熱線流速計の計測原理 21
2.3.1.3 定温度型熱線流速計 22
2.3.1.4 定電流型熱線流速計 22
2.3.1.5 温度補償 23
2.3.1.6 熱線プローブ 24
2.3.2 レーザドップラ流速計
2.3.2.1 LDAの特徴 25
2.3.2.2 LDAの原理 25
2.3.2.3 差動光方式の原理 28
2.3.3 粒子画像流速測定法
2.3.3.1 PIVの特徴 30
2.3.3.2 PIVの原理 30
2.3.3.3 ステレオPIVの概要 31
2.3.3.4 ステレオPIVの測定原理 32
2.3.3.5 撮影時間間隔の算出方法 33
第3章 燃焼を伴わない（モータリング）エキゾーストマニホールド内流れの 計測と数値シミュレーション（CFD）
3.1 緒言 35
3.2 計測方法
3.2.1 エキゾーストマニホールド入り口の流速測定 36
3.2.2 エキゾーストマニホールド入り口の可視化及び
PIV計測 37
3.2.3 触媒前集合部の可視化及び流速計測 38
3.3 実験結果及び考察
3.3.1 エキゾーストマニホールド入り口測定結果 40
3.3.2 触媒前集合部可視化結果 43
3.4 数値シミュレーション（CFD）
3.4.1 CFD手法 50
3.4.2 乱流モデル 50
3.5 CFD結果
3.5.1 CFD結果と相関係数によるPIVとの比較 51
3.5.2 相関係数によるCFD結果とPIV結果の比較 57
3.5.3 乱流モデル及びメッシュサイズの影響 58
3.5.4 触媒前の評価 61
3.6 まとめ
3.6.1 本実験手法からわかったこと 61
3.6.2 逆流に関する考察 62
3.6.3 本実験結果からわかった 62
第4章 燃焼を伴う（ファイヤリング）エキゾーストマニホールド内流れの計 測と数値シミュレーション（CFD）
4.1 緒言 65
4.2 計測方法
4.2.1 エキゾーストマニホールド入り口の流速測定 65
4.2.2 エキゾーストマニホールド集合部の流速計測 66
4.2.3 エキゾーストマニホールド入り口流速計測時の
散乱用粒子 68
4.2.4 エキゾーストマニホールド集合部流速計測時の
散乱用粒子 68
4.2.5 LDAのデータ処理 68
4.2.6 排ガス温度計測 69
4.3 計測結果
4.3.1 排ガス温度計測結果 69
4.3.2 LDA計測結果 70
4.3.3 PIV計測結果 72
4.4 CFD手法及びモデル 77
4.4.1 CFD手法及びモデル 77
4.4.2 相関係数について
4.5 CFD結果とPIV結果の比較 78
4.6 まとめ 84
第5章 商品化に対する知見の利用
5.1 緒言 87
5.2 CFDの商品化への応用
5.2.1 形状の異なるエキゾーストマニホールド 87
5.2.2 形状の異なるエキゾーストマニホールド内の
流れ 88
5.2.3 形状の異なるエキゾーストマニホールド
圧力損出 106
5.2.4 センサー検出可否 108
5.3 商品化検討 109
第6章 まとめ
6.1 まとめ 111
6.2 燃焼を伴わない（モータリング）エキゾーストマニホールド内 流れの計測と数値シミュレーション(CFD) 111 6.3 燃焼を伴う（ファイヤリング）エキゾーストマニホールド内
流れの計測と数値シミュレーション(CFD) 111 6.4 商品化に対する知見の利用 112
謝辞 113
- 1 - 第1章
序論
1-1 本研究の背景
近年,自動車開発の目指すところはめまぐるしく変化している．つながる (Connectivity)，自動化(Autonomous)，利活用(Shared & Service)，電動化(Electric) のいわゆる CASE の波が産業構造を大きく変える勢いとなっている．特に技術 革新が進んでいるのは自動運転と電動化である．自動運転は，ハンドルを握ら ずとも安全に目的の場所に到着できるという夢の自動車である．昨今では，移 動時間を利用したオフィス，ショッピング，食事が可能な乗り物が提案されて おり(1)，自動車という言葉ではなく，モビリティ(移動可能性)として提案されて いる．自動運転は色々な課題を解決する手段として，1)交通事故，交通事故死 亡者数の低減，2)交通渋滞削減によるCO2削減，3)運転者の高齢化問題解決，
4)輸送効率化などの問題解決に有効な手段と考えられている．自動運転にはレ
ベルが設定されおり，レベル 1:システムがアクセル，ブレーキまたはハンドル 操作のどちらかを部分的に行う．レベル 2:システムがアクセル，ブレーキまた はハンドル操作の両方を部分的に行う．レベル 3:決められた条件下で，全ての 運転操作を自動化する．ただし，運転自動化システム作動中も，システムから の要求でドライバーはいつでも運転に戻れなければならない．レベル 4:決めら れた条件下で，全ての運転操作を自動化する．レベル 5:条件なく，すべての運 転操作を自動化する．となっている．現在はいち早くアウディ(ドイツ)が 2018 年にレベル3の量産車を発売，国内では日産が2016年にレベル2の量産車の販 売をスタートしているが，自動ブレーキを国内に広めたのは，ぶつからない車 として評価されているスバルのアイサイトが果たした役割は大きい．自動運転 車の開発は自動車会社に留まらず，海外ではMicrosoft，Apple，Uber Technologies，
GoogleなどのIT企業が参入しており，国内でも，DeNAやSoftBankどが参入し
ている．また，自動運転には地図情報も非常に重要で，使用されている地図は 高精度地図と呼ばれている．これらは，ZENRINなどの地図メーカがいち早く手 がけているが，Toyota Research Institute – Advanced Development(TRI-AD)や，
KDDI，ZENRIN，富士通 3 社などが高精度地図を作成するための情報収集基盤
の開発を行っている．自動運転開発を取り巻く環境の国内での大きな変化点は，
2020 年だと言われている．政府は東京オリンピック開催に合わせて自動運転の 実用化を目指しており，また，次世代通信システム5Gの通信がスタートするこ
- 2 -
とにより，より一層レベル5への開発が進むと考えられる(2)〜(5)．
次世代モビリティに関する開発の大きな潮流のもう一つは，電動化である．
フォルクスワーゲン社(VW)のディーゼル排出ガス規制不正問題を発端に世界 は急速にEVシフトを始めた．フランスのマクロン大統領は，「2040年までにガ ソリン車とディーゼル車の新たな販売を禁止する」と発表している．また，ノ ルウェーが「2025 年以降にガゾリン車およびディーゼル車の新たな販売を禁止 する方針で主要政党が合意」．オランダも「2025年までに内燃機関を動力とする 車の新たな販売を禁止する案件を議会が賛成多数で同意」と発表した．これら の発言は，現段階では拘束力はなく，法規制の段階に進んではいないが今後も 世界は排出ガス規制を厳しくする方針であることが伺える(6)(7)．
すでにスタートした規制として，2018 年米国カルフォルニア州で採用された 厳しいZEV規制が挙げられる.本規制は新車販売台数の4.5％を電気自動車(EV)，
プラグインハイブリット車(PHV)，燃料電池自動車(FCV)のいずれかにするよ うにメーカに求めるものである．この規制において上記販売台数の割合は 2019 年に7％，2020年に9.5％へと厳しくなる．他の国でも同様な動きになっており．
中国ではNEV規制がスタートする．2019年から新車販売台数の10％をEV，PHV， FCV にする必要があり，2020 年には 12％となり，上記カリフォルニアの ZEV 規制よりも厳しい規制内容である．EUではCAFÉ(企業別平均燃費)という規制 がある．バッテリ充電式電気自動車(EV)は排出ガスゼロの ZEVとして優遇し，
PHEVは電気動力だけで走行できる距離に応じて CO2排出数値を優遇する．日 本国内では，経済産業省が自動車新時代戦略会議で HEV やBEV,FCEV を電動 車としてCO2排出削減の方針が話し合われており，HEVを電動車として明確に 定義しているが，欧州ではまだ電動車の定義がはっきりしない．このため，電 動車を増やす目標や規制はあるものの全ての自動車が EV になる未来はまだま だ先で，当面は内燃機関を積んだ自動車も混在する時代は続く．実際，イギリ スのエネルギー関連事業会社British Petroleum(BP p.l.c.)はEV保有台数を2040 年には3億2400万台で全体の16％と予想している(8)．また経済産業省の自動車 新時代戦略会議，パワートレイン別長期見通し(図 1)によると 2020 年に 5％の EVは2040年に15％と10％程度の伸びに対し，PHV,HV系は2020年10％だっ たものが 2040 年には 35％と 25％上昇するとの見方もしており，今後排出ガス 規制は確実に厳しくなるが，全ての自動車がEVに変わるわけではなく，一定数
がPHV,HV系へシフトされると予想される．排出ガス規制に対する技術として
は，現在までの基本的な三元触媒を使用した浄化技術，ディーゼルエンジンや 直噴エンジンから出るNOxを集める技術，DPF，GPF触媒を使用した浄化技術 以外に，低排気量ながら高出力が得られるダウンサイジングターボ，排ガス浄 化向上を狙った EGR，早期暖気によるフリクション低減を目的とした熱交換機
- 3 -
の採用など，いろいろな技術が研究開発，市販車へ採用されているが，いまだ に最も有効な手段は三元触媒を使用した浄化技術である．三元触媒を最も有効 に利用できる早期暖気が今後も重要なことは間違いない．
Fig.1 Power train another long perspective.
※国土交通省資料より引用
三元触媒の早期暖気が難しいエンジンとして上げられるのが水平対向エンジ ンである．低振動，低重心の強みがあるものの，車両レイアウト上，排気ポー トをエンジン下面から出す必要があり，直列エンジンや，V型エンジンのように 排気ポート直下に三元触媒を配置することが難しく，触媒までに長いエキゾー ストマニホールドを使用する必要があり，早期暖気にとっては不利なレイアウ トとなっている．このとき，単純にエンジンからのパイプを集合させると，パ イプ同士が向き合うためパイプの曲げ角度を大きくすることをはじめとした集 合形状を工夫する必要がある．パイプ長や集合部形状によっては，触媒の早期 暖気が不利となる.水平対向エンジンにおける早期暖気では，できるだけ早く高 温のガスを触媒まで供給することが課題となる．この課題解決のためには，エ キゾーストマニホールド内の高温時の流れの挙動を把握する必要がある．エキ ゾーストマニホールド内の流れは，流速が大きい脈動流れであり，かつ高温で あることから計測するのが難しく，数値流体解析（CFD）を利用した解析が有 効である．これまでにエキゾーストマニホールド内の流れを計測し，CFD の結 果と比較している研究例がいくつか報告されている(9)〜(16)が，その数は多くなく，
早期暖気が難しい水平対向エンジンを使用した計算例はない．さらに，CFD 検 証用の実験データは皆無である．
- 4 -
本研究は，水平対向エンジンの排気系の流動特性を明らかにするために，(1) 水平対向エンジンのエキゾーストマニホールド入口と集合部の流れの計測，
(2)CFD の確からしさを計測結果と比較し確認，(3)集合部の形状違いによる流 れの変化を明らかにすることを目標にした．これらの知見は，CFD を利用した 今後のエキゾーストマニホールドの開発に使用できるようになる．
1-2 従来の研究
自動車開発における流れの研究は主に空力という車体周りの流れ研究とエ ンジン筒内の流れ研究が主流である．下記に主な研究内容を記す．
1.2.1 車体周りの流れ研究
車体周りの流れでは主にセダンタイプやスポーツタイプの車両での研究が多
かったが(16)(17)，近年は非常に厳しい燃費規制をクリアーするため，SUVタイプ(18)
や，ワンボクックスタイプ(19)の車両開発でも盛んに研究が進んでいる．車体周 りの流れは，非常に複雑な非定常流れであるため，計測が難しかったが，近年 PIV計測(粒子画像流速計)など空間分解能が高い計測技術発展により，実験的研 究が進んでいる(16)～(19)．数値流体解析（CFD）を活用した研究開発も行われてい る(16)～(19)．空力問題は複雑な非定常流れのため，CFD における乱流モデルの選 定も非常に重要である．鬼頭等はASMO と呼ばれる 1/5 スケールの簡易車両モ デルを使用し，標準k-εモデルとLarge Eddy Simulation(LES)モデルの比較を行 い，LES が車体背後の大規模渦を含む車体周りの複雑な非定常構造を正しく再 現していることをした(20)．MenterがSASは非定常流においてLES的な挙動を再 現する(21)結果を示したことから，坂口等は同じくASMOを使用し，標準k-εモ デルとLESより計算コストの少ないScale-Adaptive Simulation(SAS)モデルの比 較を行い， SASの優位性を示している(22)．これは車両開発において非常に重要 な考え方で，LESモデルが他のモデルに比べ優れていることは理論上分かるが，
車両開発時は多種類の形状を短期間で計算，比較，選択を行う必要があるため，
CFD における精度と同時に計算時間，計算コストも考慮しながら活用する必要 がある．
1.2.2 エンジン筒内流れの研究
エンジン筒内流れの研究は，ガス流動，混合気形成などを PIV，LDA などに よる計測や，筒内流れの可視化により，タンブル流の生成や崩壊，制御の研究 を進めている(23)．筒内流動は燃焼制御にも利用できるため現在も盛んに研究が 進んでいる(24).現在市販されている自動車の熱効率は 40％程度まで上がってい るが，ここまで熱効率を上げるのには40年程度の時間を有した．しかし，2018
- 5 -
年に一旦終結した内閣府総合科学技術・イノベーション会議の戦略的イノベー ション創造プログラム(SIP)「革新的燃焼技術」においてガソリンエンジン及び，
ディーゼルエンジン共に，正味最高熱効率50％を上回ることに成功した(18)．こ のプロジェクトは開始から5 年で熱効率10％向上するという非常に短期間で目 標を達成した素晴らしい取り組みだったと言える．その中でも筒内可視化が重 要な役割を示している．横森等は，ハイスピードカメラ 2 台を用いてクランク 角度1degごとでの高速連続PIVを行っており，これによりピストン上面の境界 層は十分に発達した乱流境界層の状態に至っていないことが示された（19）．従来 のエンジン研究の多くが筒内流動は極めて強い流れを形成していることから壁 面近傍での境界層は十分に発達した乱流境界層であると考えられてきたが，そ うではない可能性を示唆した．このように，筒内流動に関しても，今後のさら なる計測，検証が必要となっている．
1.2.3 吸排気系脈動流の研究
吸排気系脈導流の研究は，出力向上や，燃費低減，排気エミッション低減な どに寄与し，極めて重要である．
吸排気系のシミュレーションには，吸気から排気までを 1 次元流体解析ソフ トを使って解析するのが一般的だが，野村等は，エキゾーストマニホールド部 の 3 次元メッシュを作成，筒内とエキゾーストマニホールドの解析を非練成と し，エキゾーストマニホールド入り口条件を一定としてエキゾーストマニホー ルドのみの評価を行っている（20）．この手法では，エキゾーストマニホールド内 の脈導流の干渉を確認できる上，k-ε乱流モデルを使用し，排圧履歴は実機とほ ぼ一致と報告している(9)．森等はエキゾーストマニホールドだけでなく，吸気マ ニホールドも 3 次元メッシュを使い吸気と排気系の 3 次元シミュレーションを 行い，EGRガス量とA/Fセンサーの空燃比について実測値と相関があったと報告 している(10)．吉沢等はエキゾーストマニホールド内にインクを添加し，触媒の インク付着範囲と CFD の結果を比較する手法で CFD の確からしさを確認した
（15）．どの報告もCFDとの比較に流速以外の値を使用して比較を行っている．直 接的にエキゾーストマニホールド内の流速を計測することは，作動流体である 空気が高温であり，流れは流速の大きい脈導流のため，困難となる．小保方等 は，2サイクルエンジンのシリンダ直前，直後の吸排気系流速測定をLDA用いて 行い，その特性を特性曲線法により評価した．LDA 計測をモータリングとファ イヤリングで実施しており，トレーサ粒子選定，計測窓の取り付け位置，計測 窓の汚れなどについて，丁寧に言及している．LDA 結果と特性曲線法は非常に 良く一致しており，単気筒の 2 サイクルエンジンに対し特性曲線法が有効であ ることを示した(11)．しかしながら，この手法も，単気筒のパイプ中心を一点の
- 6 -
みで計測にした値との比較にとどまっている．現在主に市販されているエンジ ンは他気筒エンジンであり，パイプ集合部や，触媒前の流れ場は非常に複雑な 流れとなっている．実際の開発現場ではこのようなパイプ内の複雑な流れ場の 挙動を理解し，より正確なCFDを行う必要がある
1.3 本研究の目的
本論文では，水平対向エンジンの排気系の流動特性を明らかにするために，(1) 燃焼を伴わない水平対向エンジンのエキゾーストマニホールド入口と集合部流 れの計測，(2) 燃焼を伴う水平対向エンジンのエキゾーストマニホールド入口 と集合部流れの計測，(3)計測結果と CFD の比較とエキゾーストマニホールド 開発に適した乱流モデルの選定を目的とする．本論文により得られた知見は CFD を利用した今後のエキゾーストマニホールドの開発に使用できると考えら れる．本論文では得られた結果を基に量産品開発に適用した事例も紹介する．
実験では水平対向エンジンを用い，燃焼を伴わないエンジン運転条件（以下モ ータリング），燃焼を伴うエンジン運転条件（以下ファイヤリング）の2条件と した．本研究は排ガス浄化技術の中でも特に触媒早期暖気技術向上に関する知 見を得ることを目的としているため，エンジン始動時のアイドリング条件を模 擬し，低回転域である600 rpmの定常運転を計測条件とした．エキゾーストマニ ホールド入り口は，エンジンとエキゾーストマニホールドの間に全長 100 mm のジョイントを設定し，エキゾーストマニホールド入り口の流速計測を行った．
モータリング時は，熱線流速計により計測し，ファイヤリング時には，石英ガ ラスパイプを設置し，レーザドップラ流速計 （Laser Doppler Anemometer: LDA） により計測した．またエキゾーストマニホールド集合部の流れを，可視化部を 有したエキゾーストマニホールドを用いて PIV（Particle Image Velocimetry: PIV） により流速成分計測を行った．
CFD では熱線流速計，LDA で計測を行ったエキゾーストマニホールド入り口 流速を入力条件としエキゾーストマニホールド全体の三次元CFDを行った．乱 流モデルはk-ε，SST k-ω，DES乱流モデルの合計3モデルによる比較を行った．
集合部の流れと PIV 結果を比較し、エキゾーストマニホールド内の流動特性を 最も評価できる乱流モデルとその計算コストも合わせて確認した．この手法で 確認できた内容を使って開発を行ったエキゾーストマニホールドの開発事例も 合わせて紹介する．
- 7 - 参考文献
1. トヨタ自動車株式会社ホームページ https://global.toyota/jp/company/profile/
2. 自動車技術会，“自動車用運転自動化システムのレベル分類及び定義”JASO TP 18004(2018)
3. 国土交通省，“自動運転に関する各局の取組進捗状況について”，第1会会合 資料2(2016)
4. 国土交通省，“自動運転に関する各局の取組進捗状況について”，第5会会合 資料2(2018)
5. 内閣官房 IT 総合戦略室，“ITS・自動運転を巡る最近の動向”，第１回 道路 交通ワーキンググループ 第37回 SIP第１期自動走行システム推進委員会 第 ４回 SIP第2期自動運転推進委員会 合同会議，参考資料5(2018)
6. 自動車技術会，“自動車技術 新たなステージに入った水素社会実現への取り 組み”，2018 Vol.72 No.1，p.24-30(2018)
7. 自動車技術会，“自動車技術 年鑑”，2015 Vol.69 No.8，p.2-114(2015)
8. 独立行政法人 石油天然ガス・金属鉱物資源機構，“石油・天然ガスレビュー”，
2018.9 Vol.52 No.5(2018)
9. 野村佳洋，大澤克幸，山田昌史，縄田英次，“3 次元流れ計算を用いた排気 マ ニ ホ ー ル ド の 圧 力 損 出 の 低 減 ”， 自 動 車 技 術 会 学 術 前 刷 集 No.971,p.349-352(1997)
10. 森光司，吉澤幸大，“エンジン吸排気系開発への数値シミュレーションの適 用”，自動車技術会会誌自動車技術Vol.55,No6,p.57-61(2001)
11. 小保方富夫，松尾典孝，平野嘉男，“小形二サイクル火花点火機関における 管内流の LDV による測定”，日本機械学会論文集(B 編)50 巻 450 号，
p.504-512(1983)
12. 和佐田信，遠藤正樹，“排気脈動流による熱伝達に関する研究”，日本機械学 会論文集Vol.80 No.818，p.1-12(2014)
13. 桜井雅人，田中健一，”排気脈動流れと排気吐出音シミュレーション技術”，
自動車技術会振動騒音シンポジューム
14. Atsushi Nishiyama: Application of Endoscopic Stereo PIV to 3D Exhaust Gas Flow Measurements in a Practial SI Engine，16th Int Symp on Applications of Laser Techniques to Fluid Mechanics，p.1-12 (2012)
15. 吉沢：排気マニホールド内の非定常流れ解析技術とその応用，自動車技術会 学術講演会前刷集No.975，p.129-132 (1997)倉谷尚志，川村哲裕，”高い空力
- 8 -
性能を有するフルモデルチェンジ・コンパクトセダンの開発”， 自動車技術 会2015年春季大会学術講演会講演予稿集,p.570-574(2015)
16. 農沢隆秀，岡田義浩，大平洋樹，岡本哲，中村貴樹，“自動車の空気抵抗を 増大させる車体周りの流れ構造(第2報，セダン車体の特徴的な流れ構造)”， 日本機械学会論文集(B編)，75巻757号，p83-89(2009).
17. 星田良光, 町田健太郎，金子宗嗣，小川厚，“量産SUV開発における空力抵 抗低減技術の紹介”，自動車技術会 2015 年春季大会学術講演会講演予稿 集,p.575-578(2015)
18. 館山裕之,金子宗嗣，小川厚，“広い室内空間を有するミニバンにおける空力 開発”，自動車技術会2015年春季大会学術講演会講演予稿集,p.579-582(2015) 19. 鬼頭,坪倉,”車体まわりの流れの大規模 LES -格子依存性の評価-“, 自動
車技術会論文集 Vol.38,No.6,p.297-303(2007)
20. Menter.F.R et al”A scaie-adaptive simnlation model for turbulent flow predictions,41st Aerospace Science Meeting,AIAA Paper 2003-0767(2003)
21. 坂口裕,近藤良樹,鬼頭幸三,加藤進,山本誠,”車両まわり流れ場における SASモデルの評価”,自動車技術会学術講演会前刷集No.138-08,p.5-8(2008) 22. 金子誠，“内燃機関の筒内「ながれ」解析”，日本流体力学会 ながれ 第 35
巻，p.353-358(2016)
23. 横森剛，松田昌祥，“PIVを用いたエンジン筒内の流動計測技術”，計測と 制御 第57巻 第5号，p.313-317(2018)
24. 自動車技術会，“自動車技術 燃焼技術”，2018 Vol.72 No.4，p.98-116(2018)
- 9 - 第2章
CFD理論と流速計測手法について
2.1 緒言
エキゾーストマニホールド内の流れを理解するには数値流体力学(Computational fluid dynamics :CFD)を使用した理論的解放での可視化と流速計測による実験的な可視化が理 解しやすい．本章ではCFD基礎理論と流速計測手法について基礎的な内容をまとめる．
2.2 CFD理論(1)〜(4) 2.2.1 支配方程式
パイプを流れる化学種 i（化学種とは流体を意味し，今回は排ガスを考えるため，以 後ガスとする）を考えるとき，全体を考えるとわからなくなるため，極小の一部を切り 出して考えるとわかりやすい．
Fig 2-1. Exhaust gas balance in mesh.
図2-1のような直方体（以下メッシュと呼ぶ）を流体中に考えてガスの出入りを考える．
そうすると流れによる流入量，拡散流束による流入量，反応による生成，消滅量が考え られる．時間変化のない定常流であれば，これらは三者間につり合いが成立する．
時間変化がある場合は，このメッシュ内のガスの質量の時間変化率も含まれる．図よ り，質量分率をmiとするとx方向に流れる正味流出量は
𝜕
𝜕𝑥(𝜌𝑢𝑚𝑖)𝑑𝑥 × 𝑑𝑦𝑑𝑧 x方向の拡散流束Jixによる正味拡散流出量は
𝜕𝐽𝑖𝑥
𝜕𝑥 𝑑𝑥 × 𝑑𝑦𝑑𝑧 反応による生成量は
𝑅𝑖 × 𝑑𝑥𝑑𝑦𝑑𝑧 体積内に含まれる時間変化量は
𝜕
𝜕𝑡𝜌𝑚𝑖 × 𝑑𝑥𝑑𝑦𝑑𝑧
で表せる．他の2方向y,z軸方向についても正味流出量と正味拡散流出量を考え，体積
𝐽𝑥 𝐽𝑥+𝜕𝐽𝑖
𝜕𝑥 𝑑𝑥 𝜌𝑢𝑚𝑖+ 𝜕
𝜕𝑥(𝜌𝑢𝑚𝑖)𝑑𝑥
dz dx
𝜌𝑚𝑖 dy 𝑅𝑖
- 10 - で割ると単位体積について表せる．
𝜕
𝜕𝑡(𝜌𝑚𝑖) + 𝜕
𝜕𝑥(𝜌𝑢𝑚𝑖) + 𝜕
𝜕𝑦(𝜌𝑣𝑚𝑖) + 𝜕
𝜕𝑧(𝜌𝑤𝑚𝑖) +𝜕𝐽𝑖𝑥
𝜕𝑥 +𝜕𝐽𝑖𝑦
𝜕𝑦 +𝜕𝐽𝑖𝑧
𝜕𝑧 = 𝑅𝑖・・・(1.1) ベクトル表示を用いると
𝜕
𝜕𝑡(𝜌𝑚𝑖) + 𝑑𝑖𝑣(𝜌𝑢𝑚𝑖) + 𝑑𝑖𝑣(𝐽𝑖) = 𝑅𝑖・・・(1.2)
となる．ρは密度，u,v,wはx,y,z方向速度成分である．さらに，拡散流束Jiは普通miの 勾配によって生じ，フィック（Fick）の法則に従って，
𝐽𝑖 = −Γ𝑖𝑔𝑟𝑎𝑑 𝑚𝑖・・・(1.3)
と表せる．gradはスカラー場勾配を考えた時，ベクトル微分演算子∇＝( 𝜕
𝜕𝑥1, 𝜕
𝜕𝑥2, 𝜕
𝜕𝑥3)を 導入するがこのとき，スカラー関数に作用させたものを grad と呼ぶ．これを式（1.2）
に代入すれば，
𝜕
𝜕𝑡(𝜌𝑚𝑖) + 𝑑𝑖𝑣(𝜌𝑢𝑚𝑖) = 𝑑𝑖𝑣(Γ𝑖𝑔𝑟𝑎𝑑𝑚𝑖) + 𝑅𝑖・・・(1.4)
となる．Γiは拡散係数である．divはgradをベクトル関数と内積をとる形で作用させた ものを発散といい英語表記でdiver genceと表し，𝑑𝑖𝑣𝐴 = ∇・𝐴 =𝜕𝐴1
𝜕𝑥1,𝜕𝐴2
𝜕𝑥2,𝜕𝐴3
𝜕𝑥3のことを示 す．この式は化学種i（今はガス）の保存式であり，miが濃度を示す．
2.2.2 エネルギ保存式
１メッシュについて単位面積当たりの比エンタルピhの釣り合いを考える．厳密なエ ネルギ方程式には多くの微分項が含まれるが，ここでは基本的な保存式の形にのみ注目 する．粘性摩擦係数による発熱を無視すれば，エネルギ方程式は，
𝜕
𝜕𝑡(𝜌ℎ) + 𝑑𝑖𝑣(𝜌𝑢ℎ) = 𝑑𝑖𝑣 (𝜆
𝐶𝑝𝑔𝑟𝑎𝑑 ℎ) + 𝑆ℎ・・・(1.5)
と書ける．ここでλは熱伝導率，Cp は定圧比熱，Shは単位面積当たりの熱発生率であ る．
エネルギ保存式を考える場合，従属変数として温度そのものを扱った方が便利な場合 も多い．定圧比熱Cpが一定と見なせる場合には，
𝐶𝑝𝑔𝑟𝑎𝑑 𝑇 = 𝑔𝑟𝑎𝑑 ℎ・・・(1.6) が成立するので，エネルギ式は
𝜕
𝜕𝑡(𝜌𝑇) + 𝑑𝑖𝑣(𝜌𝑢𝑇) = 𝑑𝑖𝑣 (𝜆
𝐶𝑝𝑔𝑟𝑎𝑑 𝑇) +𝑆ℎ
𝐶𝑝・・・(1.7)
と書き替えられる．式(1.7)の右辺第1 項はフーリエの法則に基づく流体中の熱伝導を 示している．
- 11 -
個体の熱伝導問題を考える場合，最大の特徴は流れが存在しないことで式(1.7)の左 辺第２項は消滅する．さらに固体は密度一定と考えてよいので，非定常熱伝導の支配方 程式は
𝜌𝑐𝜕𝑇
𝜕𝑡 = 𝑑𝑖𝑣(𝜆 𝑔𝑟𝑎𝑑 𝑇) + 𝑆ℎ・・・(1.8)
と表せる．cは個体の比熱である．流れも含めた一般的なエネルギ式(1.7)を考えれば熱 伝導問題は計算できることになる．
2.2.3 運動方程式
流れの運動方程式はニュートンの運動の第2法則，つまり，
質量×加速度＝力 より得られる．
粘性を無視した場合に１次元流れの場合，運動方程式は，
𝜕
𝜕𝑡(𝜌𝑢) +𝜕𝑢
𝜕𝑥(𝜌𝒖) = −𝜕𝑝
𝜕𝑥＋𝐵𝑥・・・(1.9)
となる．uは速度，pは圧力で左辺第1項は非定常項，第2項は対流項，右辺第1項は 圧力勾配項となり，右辺第２項は体積力(重力項)である．流体が流れるためにはこの勾 配項が必要である．この式は，オイラーの運動方程式を呼ばれる．
式(1.9)に粘性項を追加すると次のようになる．
𝜕
𝜕𝑡(𝜌𝑢) + 𝑑𝑖𝑣(𝜌𝒖𝑢) = 𝑑𝑖𝑣(𝜇 𝑔𝑟𝑎𝑑 𝑢) −𝜕𝑝
𝜕𝑥+ 𝐵𝑥+ 𝑉𝑥・・・(1.9) ここで，μは粘性係数，Vxは左辺第1項以外の粘性項を示している．
2.2.4 質量保存式
流体の密度は，エネルギ方程式(1.5)または，(1.7)によって決定する温度，(1.4)で得 られるガスの質量分率と状態方程式によって関連付づけられる．メッシュについて考え ると，正味流出入質量が，内部の質量の時間変化とつりあわなければならないため，
𝜕𝜌
𝜕𝑑𝑡+ 𝑑𝑖𝑣(𝜌𝑢) = 0・・・(1.10) が成立する．密度一定の場合を考えるならdiv u=0となる．
2.2.5 一般形保存式
これまでに見てきた式を見るとよく似た形をしている．
化学種保存式
- 12 -
𝜕
𝜕𝑡(𝜌𝑚𝑖) + 𝑑𝑖𝑣(𝜌𝑢𝑚𝑖) = 𝑑𝑖𝑣(Γ𝑖𝑔𝑟𝑎𝑑𝑚𝑖) + 𝑅𝑖・・・(1.11𝑎) エネルギー方程式
𝜕
𝜕𝑡(𝜌𝑇) + 𝑑𝑖𝑣(𝜌𝑢𝑇) = 𝑑𝑖𝑣(𝜆
𝐶𝑝𝑔𝑟𝑎𝑑 𝑇) +𝑆ℎ
𝐶𝑝・・・(1.11𝑏) 運動方程式
𝜕
𝜕𝑡(𝜌𝑢) + 𝑑𝑖𝑣(𝜌𝒖𝑢) = 𝑑𝑖𝑣(𝜇 𝑔𝑟𝑎𝑑 𝑢) −𝜕𝑝
𝜕𝑥+ 𝐵𝑥+ 𝑉𝑥・・・(1.11𝑐) 質量保存式
𝜕𝜌
𝜕𝑑𝑡+ 𝑑𝑖𝑣(𝜌𝑢) = 0・・・(1.11𝑑) 一般的な従属変数をφとおくと
𝜕
𝜕𝑡(𝜌𝜙) + 𝑑𝑖𝑣(𝜌𝑢𝜙) = 𝑑𝑖𝑣(Γ𝑔𝑟𝑎𝑑 𝜙) + 𝑆・・・(1.12)
2.2.6 メッシュの離散化
前節において，熱流体の特徴を示す速度，温度，濃度，について数学的表現で満足す べき保存方程式を導いた．しかし，これらの式は非線形項を含むため非常に特殊な条件 下を除き一般的には容易に答えを得ることはできない．そこで，メッシュを使って領域 を分け，そのメッシュ内での諸量の変化は直線的であると考えても大きな誤差を生じな い．対象領域を多数のメッシュで分解し，線形化近似を行えばメッシュ内の真の解を近 似することができる．言い換えると，メッシュの代表点における変数と隣接するメッシ ュの代表点の値を使って表すことが可能になる．
このように微分方程式をメッシュを使って領域を分け，飛び飛びの変数を使って代数方 程式に書き換えることを方程式の離散化と呼び，導かれた代数方程式を離散化方程式と 呼ぶ．
2.2.7 非定常一次元熱伝導の支配方程式と離散化
固体の熱伝導は，密度変化，流れを考慮しないので，一般方程式(1.11a)の左辺第 2 項を省略すると
𝜕
𝜕𝑡(𝜌𝜙) + 𝑑𝑖𝑣(𝜌𝑢𝜙) = 𝑑𝑖𝑣(Γ𝑔𝑟𝑎𝑑 𝜙) + 𝑆・・・(1.12)
𝜕
𝜕𝑡(𝜌𝜙) = 𝑑𝑖𝑣(Γ𝑔𝑟𝑎𝑑 𝜙) + 𝑆 となる．これは結局エネルギ方程式の(1.8)となる．
- 13 - 𝜌𝑐𝜕𝑇
𝜕𝑡 = 𝑑𝑖𝑣(𝜆 𝑔𝑟𝑎𝑑 𝑇) + 𝑆ℎ・・・(1.8) エネルギ方程式(1.8)を1次元方程式で表すと
𝜌𝑐𝜕𝑇
𝜕𝑡 = 𝜕
𝜕𝑥(𝜆𝜕𝑇
𝜕𝑥) + 𝑆ℎ・・・(2.1)
これが，非定常一次元熱伝導方程式の基礎となる．軸対象形の半径方向一次元問題の場 合には
𝜌𝑐𝜕𝑇
𝜕𝑡 =1 𝑟
𝜕
𝜕𝑟(𝑟𝜆𝜕𝑇
𝜕𝑟) + 𝑆ℎ・・・(2.2) となり，管内で考えれば半径方向の温度分布を解くことができる．
これを3 つのメッシュで離散化方程式を導くとρc を一定，時間刻みΔt メッシュの中 心間距離Δxで積分し，Δtで割ると，
𝜌𝑐
∆𝑡∫ ∫ 𝜕𝑇
𝜕𝑡𝑑𝑡𝑑𝑥 = 1
∆𝑡∫ ∫ 𝜕
𝜕𝑥(𝜆𝜕𝑇
𝜕𝑥) 𝑑𝑥𝑑𝑡 + 1
∆𝑡∫ ∫ 𝑆ℎ𝑑𝑥𝑑𝑡・・・(2.3)
𝑒 𝑤 𝑡+∆𝑡 𝑡 𝑒
𝑤 𝑡+∆𝑡 𝑡 𝑡+∆𝑡
𝑡 𝑒 𝑤
(2.3)となり，これを各項について解くと 非定常項
𝜌𝑐
∆𝑡∫ ∫ 𝜕𝑇
𝜕𝑡𝑑𝑡𝑑𝑥 = 𝜌𝑐∆𝑉
∆𝑡 (𝑇𝑃𝑛− 𝑇𝑃0)・・・(2.4)
𝑡+∆𝑡 𝑡 𝑒 𝑤
伝導項 1
∆𝑡∫ ∫ 𝜕
𝜕𝑥(𝜆𝜕𝑇
𝜕𝑥)
𝑒 𝑤 𝑡+∆𝑡 𝑡
𝑑𝑥𝑑𝑡 = 1
∆𝑡∫ {𝜆𝑒(𝑇𝐸− 𝑇𝑃
(𝛿𝑥)𝑒 −𝜆𝑤(𝑇𝑃− 𝑇𝑊) (𝛿𝑥)𝑤 } 𝑑𝑡
𝑡+∆𝑡 𝑡
= 𝑓𝑡{𝜆𝑒(𝑇𝐸− 𝑇𝑃)
(𝛿𝑥)𝑒 −𝜆𝑤(𝑇𝑃− 𝑇𝑤) (𝛿𝑥)𝑤 }
𝑛
+ (1 − 𝑓𝑡) {𝜆𝑒(𝑇𝐸− 𝑇𝑃)
(𝛿𝑥)𝑒 −𝜆𝑤(𝑇𝑃− 𝑇𝑤) (𝛿𝑥)𝑤 }
0
・・・(2.5)
生成項
1
∆𝑡∫ ∫ 𝑆ℎ𝑑𝑥𝑑𝑡 = {𝑓𝑡𝑆ℎ𝑛+ (1 − 𝑓𝑡)𝑆ℎ0}Δ𝑉
𝑒 𝑤 𝑡+Δ𝑡 𝑡
非定常項 伝導項 生成項
P (δx)w
Δx
e E w
W
(δx)e
- 14 -
= [𝑓𝑡{𝑆𝑐+ 𝑆𝑃𝑇𝑃}𝑛+ (1 − 𝑓𝑡){𝑆𝑐+ 𝑆𝑃𝑇𝑃}0]ΔV・・・(2.6) となる．
ft は時間に関する重み係数で，0〜1 の値をとる．ft=0は陽解法，ft=1/2 はクランク-ニ コルソン法，ft=1は完全陰解法と呼ばれる．今回は完全陰解法を使って計算する．
これらの式は，熱流体の特性を示す速度，温度，濃度について，数学的表現で満足すべ き保存方程式を導いているが，特殊な条件下を除くと簡単には解くことができない．そ こで，メッシュを使って領域を分け，そのメッシュ内での諸量の変化は直線的であると 考えても大きな誤差を生じない．従って，対象領域を多数のメッシュで分解し，線形化 近似を行えばメッシュ内の真の解を近似することができる．言い換えると，メッシュの 代表点における変数と隣接するメッシュの代表点の値を使って表すことが可能になる．
このように微分方程式をメッシュを使って領域を分け，飛び飛びの変数を使って代数方 程式に書き換えることを方程式の離散化と呼び，導かれた代数方程式を離散化方程式と 呼ぶ．
離散化方程式を理解した上で(2.3)に(2.4)(2.5)(2.6)を代入し，これを解くと，
{ 𝜆𝑒
(𝛿𝑥)𝑒+ 𝜆𝑤
(𝛿𝑥)𝑤+𝜌𝑐Δ𝑉
Δ𝑡 − 𝑆𝑃Δ𝑇} 𝑇𝑃
= 𝜆𝑒
(𝛿𝑥)𝑒𝑇𝐸+ 𝜆𝑤
(𝛿𝑥)𝑤𝑇𝑤+ (𝑆𝑐Δ𝑇 +𝜌𝑐Δ𝑉
Δ𝑡 𝑇𝑃0)・・・(2.7) となる．ここで，
𝑎𝐸= 𝜆𝑒
(𝛿𝑥)𝑒・・・(2.8𝑎)
𝑎𝑤 = 𝜆𝑤
(𝛿𝑥)𝑤・・・(2.8𝑏)
𝑎𝑃0 =𝜌𝑐Δ𝑉
Δ𝑡 ・・・(2.8𝑐) 𝑏 = 𝑆𝑐Δ𝑇 + 𝑎𝑃0𝑇𝑃0・・・(2.8𝑑) 𝑎𝑃 = 𝑎𝐸+ 𝑎𝑤+ 𝑎𝑃0− 𝑆𝑃Δ𝑉・・・(2.8𝑒) を使って式(2.7)を整理すると
𝑎𝑃𝑇𝑃 = 𝑎𝐸𝑇𝐸+ 𝑎𝑤𝑇𝑤+ 𝑏・・・(2.9) となる．この式を発熱がないと仮定して解くと
𝑎𝑝𝑇𝑃= 𝑎𝐸𝑇𝐸+ 𝑎𝑤𝑇𝑤+ 𝑎𝑃0𝑇𝑃0・・・(2.10)
𝑎𝑃0 =𝜌𝑐Δ𝑥
Δ𝑡 ・・・(2.11) 𝑎𝑃= 𝑎𝐸+ 𝑎𝑤+ 𝑎𝑃0・・・(2.12) となる．
- 15 - これを次のI番目のメッシュに基づいて解くと
AP(I)*T(I)=AE(I)*T(I+1)+AW(I)*T(I-1)+AO(I)*TO(I)・・・(2.13) 各係数は
AE(I)=λ/DG(I)・・・(2.14a) AW(I)= W/DG(I-1) ・・・(2.14b) AO(I)=ρ*C*DX(I)/Δt・・・(2.14c) AP(I)=AE(I)+AW(I)+AO(I) ・・・(2.14d) I＝１の境界条件は
T(1)=TO(1)=TH AE(1)=AW(1)=0.0 AP(1)=AO(1)=1.0 I=NIの境界条件は T(NI)=TO(NI)=TL AE(NI)=AW(NI)=0.0 AP(NI)=AO(NI)=1.0
となる．メッシュの数(NI 個)の連立方程式を解く必要がある．今回は効率的な解法と してトーマス・アルゴリズム TDMA(Tri-Diagonal Matrix Algorithm)と呼ばれる方法を利 用する．
Ti=PiTi+1+Qi と仮におく・・・Z Ti-1=Pi-1Ti+Qi+1
apiTi=aeiTi+1+awiTi-1+bi に Ti-1=Pi-1Ti+Qi+1 を代入する apiTi=aeiTi+1+awi(Pi-1Ti+Qi+1)+bi
Ti(api-awiPi-1)=aeiTi+1+awiQi+1+bi
Ti=((aei)/(api-awiPi-1))Ti+1+(awiQi+1+bi)/(api-awiPi-1) Zの式と同じ形式になる したがって
Pi=aei/(api-awiPi-1)
Qi= (apiT0i+awiQi-1)/ (api-awiPi-1)
となり，この式を使ってエクセルで解くことが可能になる．
3
1 2 I-2 I-1 I N1-1
DG(I-1) DG(I) X(I)
DX(I)
N1 T(NI)=TL
DX(1) DX(NI)
I+1 T(1)=TH
- 16 -
I=0 ap0T0=ae0T1+aw0T-1+b0 （T-1：-1は温度もないメッシュのため意味なし）
1 = 0 + 0 + ap0Tp00
I=1 ap1T1=ae1T2+aw1T0+b1
・・
I=N-2 apN-2TN-2=aeN-2TN-1+awN-2TN-3+bN-2
I=N-1 apN-1TN-1=aeN-1TN +awN-1TN-2+bN-1
apN-1TN-1= 0 + 0 +apN-1TpN-10
1 TN-1= 0 + 0 + 1 TpN-10
2.2.8 流れおよび熱伝導の数値解析法
流動解析では新たに圧力と各速度成分が未知数として加わるが，数値解析手法として は境界条件および初期条件の設定が容易でかつ三次元流れへの拡張も容易な，SIMPLE 法(Semi Implicit Method for Pressure-Linked Equations)を採用した．
2.2.9 時間平均方程式の誘導
支配方程式としては質量保存式（連続の式），運動方程式（ナビエ・ストークスの式）
およびエネルギ方程式を考える必要がある．これらの式を非圧縮性流体と考え，流体密 度ρ一定として，テンソル形式で書き直すと次になる．
𝜕𝑢𝑗
𝜕𝑥𝑗= 0・・・(3.1)
𝜕𝑢𝑖
𝜕𝑡 + 𝜕
𝜕𝑥𝑗(𝑢𝑗𝑢𝑖) = 𝜕
𝜕𝑥𝑗(𝜈𝜕𝑢𝑖
𝜕𝑥𝑗) −1 𝜌
𝜕𝑝
𝜕𝑥𝑖+ 𝑔𝑖𝛽(𝑇 − 𝑇𝑟𝑒𝑓) + 𝜕
𝜕𝑥𝑗(𝜈𝜕𝑢𝑗
𝜕𝑥𝑖)・・・(3.2)
𝜕𝑇
𝜕𝑡+ 𝜕
𝜕𝑥𝑗(𝑢𝑗𝑇) = 𝜕
𝜕𝑥𝑗(𝜈 𝑃𝑟
𝜕𝑇
𝜕𝑥𝑗)・・・(3.3) ここで ν は流体の動粘性係数(ν=μ/ρ)，βは体積膨張率(𝛽 = − (𝜕(𝑙𝑛𝜌)𝜕𝑇 )
𝑝)，Pr はプラン
トル数(Pr=μCp/λ)という無次元数である．また，温度さによる浮力は，
𝑔𝑖𝛽(𝑇 − 𝑇𝑟𝑒𝑓)
より近似（ジプネ近似）しており，運動方程式(3.2)の右辺第3項，第4項は(1.9)式の 体積力項βxおよび粘性項Vxに対応している．
さて，ここで実際に時間平均の操作をするに当たり各変数を 𝑢𝑖 = 𝑢̅ + 𝑢′𝑖 𝑖・・・(3.4𝑎)
𝑇 = 𝑇̅ + 𝑇′・・・(3.4𝑏) 𝑝 = 𝑝̅ + 𝑝′・・・(3.4𝑐)
時間平均化された成分𝑢̅ , 𝑇̅ , 𝑝̅𝑖 とそのまわりの不規則成分u’i ,T’ ,p’との和に置き換える．
重力ベクトル 体積膨張率 参照温度のズレ ref は参照の意味
