予混合乱流伝ぱ火炎構造の解明と直接数値シミュ レーション

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

予混合乱流伝ぱ火炎構造の解明と直接数値シミュ レーション

城戸, 秀樹

九州大学工学機械科学

https://doi.org/10.11501/3150820

出版情報：Kyushu University, 1998, 博士（工学）, 課程博士 バージョン：

権利関係：

予混合乱流伝ば火炎構造の解明と 直接数値シミュレーション

城戸秀樹

目次

L i s t  o f  T a b l e s  and F i g u r e s  

m 

記 号 V 

1 1 1 3 4 6   1.  緒 論

1.1  燃 焼 と 乱 れ

1.1.1 

乱れの構造

1.1.2 

乱流火炎モデ、ル

1.2 

数値解析手法

1.3 

本論文の概要

づ / 寸 / づ / O o n JyハU今ムA守

噌EEA唱EEA唱︐EA

2. 

計 算 手 法

2.1 

スベクトル法

2.1.1 

重み付き残差法

2.1.2 

スベクトノレ法による離散化

2.2 

フーリエ・スベクトノレ選点法

2.3 

ギブス現象

2.4  非線形項の取り扱いとエイリアシング誤差 2.5  時間前進法

/OfOQV

ハ

UAV‑1‑i

今ノ

h A

時E

F3

1111

今 ム 今 ム 今 ' 'u 勺 ん 弓 ノ ム ウ

&

今 ん 今 ん

3. 

計算条件

3.1 

支 配 方 程 式 3.2  物 性 値

3.2.1 

当量比

3.2.2 

粘 性 係 数

3.2.3 

熱伝導率 3.2.4  拡散係数

3.2.5  定圧比熱及びエンタノレピー 3.2.6  その他の物性値

3.3 

初期乱れ場の発生方法

QノQ

ノ 今 ム 今/ 何 今 ん 今 コ

4. 

乱れ強さと乱れのスケールの影響

4.1 

計算方法及び初期条件

4.2  燃 焼 反 応

4.3 

物 性 値

3 2  

4.4 

計算結果

3 3  

4.4.1 

乱れ強さの影響

3 3  

4.4.2 

乱れのスケーノレの影響

3 3  

4.5 

結論

3 8  

5. 

分子の選択拡散の影響

39 

5.1 

計算方法及び初期条件

39 

5.2 

燃 焼 反 応

40 

5.3 

物 性 値

40 

5.4 

計算結果

42 

5.4.1 

当量比1.

0

の場合

42 

5.4.2 

当量比

0 . 7

及び1.

2

の場合

48 

5.5 

結 論

5 1  

6. 

局所的な当量比変化を考慮した反応機構

5 2   6.1 

計算方法及び初期条件

5 2  

6.2 

燃焼反応

5 2  

6.3 

物 性 値

54 

6.4 

計算結果

54 

6.5 

結 論

6 3  

7. 

結論

64 

参考文献

66 

謝 辞

74 

付 録

R e d u c e d  k i n e t i c  m e c h a n i s m   ^{7 5}  

11 

L i s t  o f  T a b l e s  and F i g u r e s  

T a b l e s  

Table 3.1  Coefficients to calculate properties  Table 4.1  Properties of mixture used by calculation 

Table 4.2  The relation between turbulence scales and turbulent burning velocities  Table 5.1  Properties of mixture used by calculation (か0.7)

Table 5.2  Properties of mixture used by calculation (o= 1.0)  Table 5.3  Properties of mixture used by calculation (砂=1.2) Table 6.1  Properties of mixture used by calculation (砂=1.0) Table 6.2  Properties of mixture used by calculation (か1.2)

Fig.2.1  Fig.2.2  Fig.3.1  Fig.4.1  Fig.4.2  Fig.4.3 

Fig.4.4  Fig.4.5  Fig. 5.1  Fig.5.2  Fig. 5.3  Fig. 5.4  Fig.5.5  Fig.5.6  Fig.5.7  Fig.6.1  Fig.6.2  Fig.6.3  Fig.6.4 

F i g u r e s  

Dirichlet kernel for N=8， N= 16 

Several smoothings for the square wave  Comparison ofCalculated data with JANAF  Initial condition of vorticity distributions 

Initial condition of temperature and mass 企actionof reactant  The effect of turbulence intensities on the propagating flame  (k_fllax=2mm‑^l， t=30μs) 

Comparison between measured data and simulated results 

The effect ofturbulence scales on the propagating flame (u'=20m/s， t=50μs)  Initial condition oftemperature and mass fraction (o=1.0) 

Distribution of mass企actionat t=20μs (o=1.0)  Variation of mass fraction at t=20μs(o=1.0)  Distribution of equivalence ratio at t=1 0μs (o=1.0)  Distribution of equivalence ratio at t=20μs (砂=1.0) Distribution of equivalence ratio at t=20μs 砂(=0.7) Distribution of equivalence ratio at t=20μs (砂=1.2) Initial condition of temperature and mass fraction (o= 1.0)  Variation of laminar burning velocity with equivalence ratio  Distribution of temperature and reaction rate (ケ1.2，DJ..  D₂) 

Distribution of temperature and turbulence (ケ1.2)

Fig. 6.5  Distribution of mass fraction at t=80μsc

o = 1 .

0)  Fig. 6.6  Distribution of mass fraction at t=80μsC

o = 1 .

2) 

Fig.6.7  Distribution oftemperature and mass f1ux at t=40μs C

か1.

2) Fig. 6.8  Distribution of equivalence ratio at t=80μsc

o = 1 .

0) 

Fig. 6.9  Distribution of equivalence ratio at t=80μsC

o = 1 .

2) 

lV 

号

司い コ 一

一口

本研究で使用した主な記号と定義を以下に記す.

B 

:反応の頻度因子[lI

s ]

C₁ 化 学 種iの質量分率[ー]

Cp，l 化 学 種iの定圧比熱

[ J / m o l .K ]  

DI 化 学 種iの拡散係数[m²/s]

Dυ  :化学種 i，化 学 種j聞の相互拡散係数[m²/s] e

，

全エネルギー

[ J / k g ]

E  :活 性化エネルギー

[ J / m o l ] E ( k )  

:乱れのエネルギースペクトノレ

ム

:化 学 種iのエンタノレピー

[ J / m o l ] k 

:ボノレツマン定数[J/K]

k  :波 数

[mm ‑

^I]

k_max 乱れのエネルギースペク トノレが最大となる波数[mm‑^I]

L  :乱れのマイクロスケール[m] Lf  :積分長さスケーノレ[m]

~ 化学種 i の分子量[kg/kmol]

p  :圧力[Pa] Pο  :初期圧力[Pa] p，^附 :理論燃焼圧力[Pa]

Q  :発 熱 量

[ J / k g ]

ん :燃焼室実体積の等価球半径[m] R^O  • 一般ガス定数[J九nol.

K ]  

Sしo 層流燃焼速度[m/s] S1 ' 乱 流 燃 焼 速 度[m/s] t  :時間[s]

企t :時間刻み幅[s] T  :温 度[K]

T*  .無 次 元 温 度[ー]

u 速度ベクトノレ[m/s]

u' 乱れ強さ[m/s]

U  :内部エネノレギー[J/Kg]

v

，

拡散速度[m/s]

~ 化学種 i のモノレ分率[-]

6/^，0 層流火炎帯厚さ[m]

~x) :デイラクのデ、ノレタ関数 ε  :エネルギー消散率[m²/s^3]

ε  :レナード ・ジョーンズのポテンシヤノレの極小値[1]

ゆ :当量比[ー]

ηk コルモゴロフスケール[m]

K 比熱比[ー]

λ  :熱伝導 率[W/m.K]

必

:テイラーのミクロスケーノレ[m] μ  :粘性 係数[kg/m.sec] 

v 動粘度[m²/s] θ  :無 次 元 温 度[‑]

ρ  :密 度[kg/m³]

σI 分子 iの特性直径

[ A ]

ω  :反応速度[kg/m^3・s]

0/) 衝突積分[ー] Qμ  :衝突積分[ー]

Vl 

1  .緒論

近年，化石燃料の枯渇等のエネノレギー供給問題や二酸化炭素増加による地球 温暖化現象等の環境汚染問題が大きく取り上げられている.このような背景の もと，現在のエンジン設計は従来までの高性能化，高効率化に加えて，低公害 化も要求されている.このような複雑な要求を満たすためにはエンジン稼働時 の全現象を把握する必要がある.このため，内燃機関において，重要な現象で ある予混合乱流火炎の研究はその領域，内容ともに増加している.火炎構造の 定量的な領域区分については，まだ多くの課題が残されているが，次第に明ら かになりつつある.

一方，流体分野で発達した数値シミュレーションは，計算機の高速化に伴い 乱流燃焼場についての解析ができるようになってきた.乱流燃焼場についての 数値シミュレーションは乱流火炎に影響を及ぼすと考えられる支配方程式中に 含まれる因子の影響を，それぞれ独立に調べることができるため，実験では困 難な点の検討が行える.

本章では，これまでの乱れの構造や予混合乱流火炎に関する従来の研究及び 数値解析的手法の発展を要約し，本研究の背景を述べ，本研究の目的を述べる.

1.1  燃焼と乱れ 1.1.1  乱れの構造

乱流場においては，粘性力に対して↑貫性力が支配的となってくる.Osborne  Reynoldsは乱流場においては，乱れの変動速度 (fluc国atingvelocity)により，分 子運動で生じる粘性力に匹敵する慣性力を生じることを示した.これは， レイ ノルズ応力 (Reynoldsstress) と呼ばれる.乱流場で、は流れに沿って次々に渦粒 子を発生し，平均流れの運動エネルギーを渦エネルギーに変える.大きい渦は 分裂して小さい渦にエネノレギーを伝え，小さい渦の渦エネルギーは粘性により 摩擦熱に変えられる.大きい渦から小さい渦への乱れのエネルギー輸送過程を

カスケード・ダウン (cascadedown)と呼ぶ.したがって，流れには大小さまざ まな渦が含まれ，大渦の中に小渦が，さらにその中に微細渦が含まれている. このように，乱流場においては，渦の大きさや強さがエネルギー伝達に重要で、 ある[114].

乱れは，広い周波数スベク トノレを持ち，構造が非常に複雑である.その性質 を表すために統計的手法が用いられる. Taylorの等方性乱流の理論[115]による

と，相関関数g(y)は，

ふ

( y ) =

1

ポ

^{E'︑/︑噌EEA E ‑A  ︐ ︐︐︑︐ ︐ ︑}

で表せる放物線に接する.ん は等方性乱れのテイラーのミクロスケール (Taylor microscale) である.乱流運動のエネルギー

k=

炉 + 手

+w² )/2  (1.2)  は流速勾配により作られ，小渦?微渦を経て分子運動へと消散して行く.ここ で

，

U

， 

V， Wはそれぞれ

，

X

， 

y

， 

Z方向流速の変動成分である.テイラーのミ クロスケーノレと等方性乱れのエネノレギー消散率e=‑dk/dt(energy dissipation rate)  との聞には次のような関係がある.

ε= 15v竺マ λL g 

(1.3) 

ここで， iは乱れ強さ (turbulenceintensity)， vは動粘度 (kinematicviscosity)  である.

テ イ ラ ー の ミ ク ロ ス ケ ー ル は エ ネ ル ギ ー 消 散 率 と 密 接 な 関 係 を 持 つ が ， Kolmogorovは一歩進め，エ不ルギー消散率と動粘性係数だけを用いて，コルモ

ゴロフスケール (Kolmogorovscale，ηk )を定義した[115].すなわち，

ηk = 

I

¥ &  ) と

r

^(1.4) 

である.このスケーノレは，分子運動に変わる寸前の微細渦の代表寸法と見られ，

粘性の影響でこれより小さな渦は存在し得ないと考えられている.

等方性乱流場において，乱れのマクロスケーノレL(macroscale of turbulence)が 次式で定義される[119].

ε=At~')(与)

^(1.5) 

ここで， Alは約 11である.

Lは大きな渦の代表半径， λgはこの渦に含まれる小渦の代表半径と見なすこ

とができる.大渦が小渦とエネノレギー交換する割合は，乱れのエネルギーに比 例し，渦運動の特性時間に逆比例する.したがって，Lと令の関係は，

2 

L  A^{，0 4 n}  m  一一=ーづ=ごK.e

λR

、

^{/15 (1.6)} 

である.ここで，Reはレイノルズ数 (Reynoldsnumber) である.

これらのことから，乱流中には仇 ~L のスケールの渦が含まれていると考え られる[115].

1.1.2  乱流火炎モデル

燃焼に及ぼす乱れの影響は以前から定性的には知られていた.

1 8 8 3

年に MallardとLeChatelierは，乱れはエネノレギーの伝達を促進し，火炎表面積を増加

させ，新しい点火点を作ると述べている[39].

このような乱れの影響は 1947年にDamkohlerにより，より定量的に考察され た.Danu(ohlerの仮説[2]によれば，流れ場に存在する渦のスケーノレが層流火炎厚 さと比較して小さい場合，乱れにより反応域での伝達率(仕組sportrates)が増加 し，渦のスケールが層流火炎厚さより大きかった場合，乱れによるしわにより 火炎面積が増加する.このように，予混合乱流燃焼は乱れのスケールと層流火 炎厚さを比較することにより，おおまかに2つの燃焼モデルに分けることがで きる.それらは，分散反応火炎モデル (dis廿ibutedreaction model)としわ状層流 火炎モデル (wrinkledlaminar f1ame model) と呼ばれている. しかしながらこれ

も1つの仮説に過ぎず，多くの議論がある

[ 9

之

1

，

1 1 2 ] .

Kovasznayはより包括的に火炎構造における乱れの影響について述べている [40].乱れのスケーノレに加えて，重要な乱れのパラメータである乱れ強さj，重 要な混合気の物性値である層流燃焼速度 SLO'火炎構造を判断する基準である無 次元パラメータのコバツニー数 (Kovasznaynumber)  [40]を考慮に入れる必要が あると述べている.コパツニー数Kzは化学反応特性時間と乱流運動特性時間と を比較した無次元パラメータである.

Kz‑u

y λ g 

δLO/S_ω 

ここで， δLOは層流火炎帯厚さである.

(1.7) 

KlimovとWilliamsは，乱流においてしわ状層流火炎 (wrinkledlaminar f1ame)  が存在する条件として， η仇k>δ

九

LO (恨Klinlo仰V

Baι11凶alとLefe肘bvr印eらは，乱流火炎モデ、ノレにおけるコルモゴロフスケールの重要 性に着目し，未燃焼混合気中の乱れの諸特性の詳細な計測と乱流火炎の瞬間シ ュリーレン写真に基づいて， U'/SLO'η

' k /

δLOの大小関係、により火炎構造が異なる ことを主張し， 3領域モ、ノレを提案したデ [2.3，4].

乱流燃焼をいくつかの基準をもと に異なった領域に分類する試みは，さらに

発展して，無次元パラメータを座標とした平面上で燃焼領域に分ける燃焼モデ ルの線図表示に至った.無次元パラメータとしては，レイノルズ、数 (Reynolds number) ，ダムケラー数 (Damkohlernumber) ，コパツニー数 (Kovasznaynumber) ， 

仇/δLO' U'jSLO' Lf / ð

LO

などがある.ここで，んは縦方向積分尺度である• Williams  は，U'jSLOと仇

δ /

^LOを座標軸に適用し，平面上にstirredreaction， weak turbulence，  reaction sheetの領域を示した[112].同様な試みは続いて Bray，Peters， Abrahan1，  Williams and Bracco， Kidoらによってなされた[9，22，43，44]. これらの線図表示法 は領域間の基準値の特性など実験結果との比較に関してさらなる研究が望まれ る.

1.2  数値解析的手法

乱流燃焼現象は，流動，化学種の拡散，反応，反応による熱発生，熱の移動，

などであり，これらの現象が並行して生じる.従って，舌し流燃焼現象は燃焼反 応を伴う流れ場とみなし，数値流体力学的手法が適用できる.すなわち，圧縮 性ナビ、エ・ストークス方程式を数値解析する際に，化学反応を考慮に入れてシ

ミュレーションを試みる.

数値流体力学の計算手法には，差分法，有限要素法，境界要素法，スペクト ル法などが挙げられる.それぞれの特徴を簡単に述べると，差分法は，微分方 程式に現れる微分項を，テイラー級数展開を利用して表現する手法である.有 限要素法は，基礎方程式をいったん重み付き残差方程式に変換し，これを近似 するという間接的手法を用いる.また，差分法とは異なり，分割に用いられる 有限要素の形や大きさや向きは任意に選択することができる.境界要素法は，

境界を有限個のセグメントに分割し，各セグメント上で未知数を離散化する手 法である.スペクトル法は，関数を直交関数により級数近似をする[116]. 

近年，スベクトノレ法が，少ない格子点で良い精度の計算を行えるため，様々 な分野で使われてきている.スベクトル法は，差分法，有限要素法などの数値 解法と同レベルの重要な方法になっている. しかし，スペクトル法が流体力学 問題へ多く応用されるようになったのは， 1970年以後，特に 1980年代のことで ある. 1977年に，スペクトル法の理論は， Gottlieb， Orszag[116]により展開さ れた.それとは対照的に， Canutoらは，理論から応用に至るまで広い範囲にわ たり，スペクトノレ法のほとんどについて触れている[116].その後，スペクトル 法による様々な流体計算が行われてきた.層流に関しては， Orszagによるテー ラー ・グリーン渦の計算， Moin， Kim，や Zang，Hussainiによるチャンネル流 れの計算， Moserらや Marcusによるテーラー・クエット流れの計算， Quere，  Roquefortによる自然対流の計算， Hatziavran1idis， Kuの浮力により誘発された流 れの計算， Kuらによるキャビティ流れの計算， Tanによる層流馬蹄渦流れの計

4 

算， Goldhirschらによる熱対流の計算などがある[116].また安定性に関しては，

Patera， Orszagによる軸対称パイプ流れの計算， Bramley， Dennis， Bridges， Morris，  Gardnerらによるポアゼーユ流れの計算， Orszag， Pateraや Bayley，Orszagによ

る努断流れの計算， Ha11， Malikによる付着線境界層の計算などがある[116].さ らに，遷移に関しては， Orszag， Kellsによる 2次元ポアゼーユ流れと 2次元ク エット流れの計算， Laurien， Kleiserによる境界層の計算などが[116]，また，乱 流に関しては， Orszag， PattersonやBasdevantによる等方性乱流の計算， Fornberg  による 2次元乱流の計算， Brachetらによるテーラー・グリーン渦の計算， Eidson 

らによる自然対流の計算などがある[116]. 

スベクトル法は乱流燃焼場における数値シミュレーションにも適用されてし る[97]"‑' [108].乱流燃焼に適用する場合には，燃焼反応を記述する化学反応に関 しでも，解析手法に適した取り扱いをしなければならない.燃焼現象に関与す る化学反応機構は非常に複雑である.水素やメタン等に関してはほぼ正確な化 学反応機構(完全化学反応機構白11kinetic mechanism，詳細化学反応機構 detailed kinetic mechanism)が明らかにされている.これらの詳細化学反応機構は数十 数百の化学種とそれ以上の素反応から構成されている.水素・酸素予混合火炎 及び拡散火炎における詳細化学反応機構は，窒素が反応に関与しない場合にお いても最小で 8種類の化学種と 21組の素反応から構成されている.従って，水 素・酸素予混合乱流火炎あるいは乱流拡散火炎の数値解析を行う場合，連続の 式，運動量保存式，エネノレギ一保存式に加えて 7化学種に対する保存方程式を 同時に解析しなければならない.さらに，反応機構を構成する多くの素反応の 時間スケールは，流体運動の時間スケーノレよりはるかに短いため，乱流燃焼の 解析を困難にしている.

これらのことから，現在，比較的重要度の高くない化学種が関与する素反応 については定常状態を仮定することにより，化学反応機構を簡単にする試みが 行われている.これらは 簡略化学反応機構 (reducedkinetic  mechanisln) と呼

ばれており，乱流燃焼の実用的な数値解析には必要不可欠である.

偏微分方程式の解法は，重み付き残差法により表すことができ，この重み関 数をどのように表すかで様々な方法が考えられる.スベクトル法で、は，直交多 項式を用い，無限回微分可能である.主に，スベクトノレ法は以下の 3つの手法 に分けられる[116].スベク トノレ・ガラーキン法 (SpectralGalerkin Method)は級 数展開によって係数問の関係式を導き，それを解く方法である.スベクトル・

タウ法 (SpectralTau Method)は境界条件についても級数展開し，係数聞の関係、

式を同時に解く方法である.スペクトノレ選点法 (SpectralCollocation Method)は

級数展開された関数を離散点であわせる方法である.

本論文ではフーリエ ・スペクトノレ選点法を用いて，予混合乱流伝は火炎の直 接数値シミュレーションを試みる.

1.3  本論文の概要

以上のように予混合乱流燃焼に関する研究は，多くの研究者によりなされて おり，かなりの成果を上げている.本論文では，予混合乱流燃焼の直接数値シ ミュレーションを行い，乱流燃焼現象の解明を試みている.本論文の 2章以下 は次のような構成である.

2章では，直接数値シミュレーションに用いた計算手法，スベクトル法(主 にフーリエ ・スペクトノレ選点法)及びルンゲ、・クッタ ・ギノレ法について述べる.

3章では，直接数値シミュレーションに用いた支配方程式，物性値及び初期 乱れの発生方法について述べる.

4章では， 2次元予混合乱流伝ぱ火炎において，フーリエ ・スペクトル選点 法による直接数値シミュレーションを行う .乱れの火炎に及ぼす影響を調べる.

5章では，燃料と酸素の分子拡散を考慮、した予混合乱流伝ぱ火炎の直接数値 シミュレーションを行い，分子拡散の火炎に及ぼす影響を調べる.

6章では，局所的な当量比変化を考慮、できる反応速度の式を提案し，予混合 乱流燃焼において直接数値シミュレーションを行う .

7章では，本研究で得られた計算結果をまとめて提示する.

2. 計算手法

2.1  スベクトル法

直接数値計算を行うためには，精度の高い離散化手法が要求される.スベク トノレ法は，差分法や有限要素法などに比べ，高い精度を有するので多くの直接 数値計算に利用されている[110]. スベクトノレ法にはガラーキン法，選点法，タ

ウ法などがある.スベクトル法で、は，非線形項の取り扱いに注意を要し，擬ス ペクトル法と呼ばれる手法を用いることにより計算時間を短縮することができ

る.その場合，エイリアシング誤差を取り除く必要がある.

2.1.1  重み付き残差法

スベクトル法は重み付き残差法の一種である.重み付き残差法においては，

試行関数(展開関数，近似関数)及び試験関数(重み関数)が重要である.試 行関数は微分方程式の解を有限級数に展開する場合の近似関数で，試験関数は 微分方程式が有限級数展開によって満足されることを保証するために用いられ

る.

ある関数u(x)を (xの範囲はα β)以下のように有限級数関数utlx)

u(x)~

U N ( X ) =   L

^{αnOn (x)  (2.1)} 

n=l 

で近似する場合を考える.ここで α11 は未知係数，ゆ~(x) は基底関数で，多項式，

三角関数，チェビ、シェフ関数，ノレジャンドル関数などの収束性の高い関数列か ら適当なものを採用する.

重み付き残差法とは，残差関数

R ( x )

と試験関数w^/l

l x )

とのα β問で、の内積が O になるように係数列引を決定する方法である.

(

収 R ( 仇

ベ¹

ω

ふ

^ゆ

刷^x

^寸

(

^ゆ

か

^{) )   =} : f ^{R (}

^{かx中以仇か仇九}

^パ

^m

^川

^ベ1(

ここでで、，試験関数WIl1(x)の選択により，重み付き残差法は以下のように分類で きる.

( 1 ) 選点法

W m ( x )  

=δ(x ‑X/l

J 

(2.3) 

e ¥

^x)は，デイラクのデルタ関数で，xの範囲α βの適切に選択された N点でu(x) とutlx)が等しい値をとるという条件を与えることに相当する.

(  2 )  

最小自乗法

W_I1I

( x )   =出 _。

_α ^(2.4) 

最小自乗法で，残差関数

R ( x )

の自乗のα β問で、の積分値が最小となる条件を 与えることに相当する.

( 3 ) ガラーキン法

W m ( x )  

^=ν111 

( x ) o m  ( x )  

(2.5)  試験関数として基底関数を用いる.νぷx)は，仇I令)に正規直交性をもたせるた めの重みで，関数列によって以下のように定義される.

一角関数

チェビ、シェフ多項式 ノレジャンド、ル多項式

( 4  ) 

モーメント法

W m ( x ) = x m ‑ l  

九

j ( x ) =   π

^川

L ^{( 2 . 6 )}   m ( x ) =   。 ^{‑ X 2 } 2 ( 2 . 7 )}  

ν m ( x )   = 

1 

( 2 . 8 )  

(2.9) 

いずれの場合も，式

( 2 . 2 )

は

N

元の連立方程式に帰着する.この連立方程式を 解くことにより，係数列 αnが決定し，近似関数が得られることになる.この場 合，基底関数列と試験関数列が直交関係，

(

仇 ( x )

^，

W m ( x ) ) = 

0 

( m

^;t: 

n )  

(2.10)  にあると，連立方程式を表す係数行列は，対角成分のみが零でなくなり，連立 方程式を解く作業がきわめて容易になる.

基底関数列は 近似すべき関数の性質に応じて選択すべきであるが，一般に，

計算が容易であること(直交関数列，高速変換法 [FFT等]が利用できる)， N 

→∞で近似誤差が 0に収束すること，小さい Nで精度の高い近似が得られるこ と，などの性質を持っていることが望ましい.

2.1.2  スベクトル法による離散化

離散化においては，関数の近似手法と同様の操作を微分方程式に施すことに なるが，試験関数の選択，境界条件の取り扱い方によって以下の3種類に分類 できる.

( 1 )ガラーキン法

試行関数は境界条件を満足する.たとえば，周期境界条件の場合，三角関数

を試行関数として選択することで u

l A

X)の周期性は自動的に満たされる.

( 2 )選点法

計算領域中の N個の選点上で，近似誤差を 0とすることに相当する.通常，

境界も選点とする.

( 3 )タウ法

試行関数は境界条件を満足しないので，境界条件の数だけ展開の次数を増し，

境界条件を満足させるための式を付加する，

2.2  フーリエ・スベクトル選点法

本研究では，試行関数に三角関数で展開するフーリエ級数を使用するフーリ エ・スベクトル選点法を用いた.フーリエ・スベクトノレ選点法の計算手法は以 下のとおりである.0'""2π間において，関数u(X)は，

u

か ) = : L ^{u 〆 伽 ( 0}

^壬X三

^{2 π )}

^(2.11 ) 

k=ー∞

のように展開できる.一般に，数値計算では N 個の有限項で打ち切り近似計算 を行う. したがって，式(2.11)は次のように近似できる.

N /2‑1 

u(X)= 

: L ι e

^伽

k=‑N 12 

ただし，_Ukは級数係数で，

u_k =

会 ^{f o}

^27r

^{U ( X )}  

^e‑I 

である.次に，

X. = _‑2rcj 

J  N 

として式(2.12)の離散化を行う.

フーリエ変換 (FFT) を用いて，

似値

u J .

毛)及び級数係数らは，

山 j ) = 三 ん

^I灯j

k=‑N'.2 

Uk =

岩山 ^{j h ‑ I h}

となる.式(2.15)のn回微分は，

ぞ:'In N/2‑1 

(2.12) 

(-N/2~k 三 N/2-1) (2.13) 

j = O，

t " '

，N ‑1  (2.14) 

この場合，離散点巧が選点となる.これを高速 フーリエ変換を行うと，離散点における級数近

(2.15)  (‑N /2三k_~ N /2‑1)  (2.16) 

士山 j ) = 乞

(ik

Y 

^{ukelkXj  (2.17)} 

V A   k=‑NI2 

となり，空間微分を解析的に行うことができる.この解析的微分が行える点が スベクトル法の大きな特徴の一つである.

2.3  ギブス現象

ギブス現象は 関数の不連続点付近に現れるフーリエ級数・離散フーリエ級 数に特徴的な振動現象である.含まれる波数が大きいほど，不連続点に近い最 大の振動強さはある限界値に近づき，オーバーシュートは不連続点に接近する.

これは次のように説明できる.関数U(X)のフーリエ級数展開は，

九

U(X)

= 工 ら

^e山

Ikl~N 2 

x u

ノd

JG   ρしU︑hE

︐ /

bJ

/v

  vu 吋寸￨川￨￨﹂

レ

Fザ

y

い

/'

t︑ ︑

k

u e

2 π︑d

J

‑ O

一 ︑

/

‑N

{lJ

宮 ︑ ︐

A

︿ 一

l一π﹁￨￨九州‑

pi一

今 ム

π

2 ρ lJ

o 

1嗣

y ]

一 勾

一 一 一 一

(2.18) 

である.よって

P N U

は

P N U ( X ) = 去 f

^Jr

^{D N ( X} ー ル ( y ) d y

のように表すことができる.ここで，

D N ( ω c )

はDi凶d巾ωchωh吐制ll

N 2 

D_N(c)=1+2Lω kc 

(2.19) 

k=l 

r 

sin((N + 1);/2) 

= ~ sin(c/2)  IN+1 

c 

^=1‑2~j

jEZ  (2.20) 

c = 2

^勿.

で表される.波数N=8，16の場合における Dirichlet核をFig.2.1に示す.

24 

N =   8 

N=16 

1 6 

D N  8 

‑8 ‑1 .0  ‑0.5  O  XIπ 

0.5  1 .0 

Fig. 2.1 Dirichlet kernel for N=8， N=16 

10 

離散フーリエ級数の場合は，

仰い)=万ーさんか -Xj~(Xj)

^(2.21 ) 

のようになる.

理論的にも実用的にも，物理空間における不連続関数の忠実な再現のために，

不連続点付近におけるギブス現象を小さくすることが望まれる.ギブス現象は 不連続関数におけるフーリエ係数の減衰の遅さに関係があるため，高次の係数 にスムージングを施す. しかし，あまりに強いスムージングは結果として関数 の近似の精度を下げる.そのため，適切なスムージングの方法が考慮、されてき た.

高次のフーリエ係数にσkをかけることによりスムージンク令する方法が，最も 簡単な方法である.これにより，フーリエ級数PNUは，

SNU(

χ ) 二 ヱ

^{σAGIK (2.22)} 

k=‑N /2 

のように置き換えられる.通常， σkは負でない実数値をとり， σ。二 0，σk二 σ‑k

である.

Cesaroはフーリエ級数の古典的なスムージング法で，

σ

ι 二 1 ̲ ̲1 k l k=‑21..N l 

凡 N/2 + 1  2' . 2  (2.23) 

である.その他にも， Lanczosによる方法，

σ sin(2kπ/N) 

k ‑ 2kπ/N 

N一2

N一2

v k   (2.24) 

Raised cosineによる方法，

σk=l + c o s ( 2 K M N ) k二 ‑ N N  l  (2.25) 

2  2  2 

がある.Figure 2.2に矩形波の離散フーリエ級数展開における， (a)スムージング なし， (b)Cesaroによるスムージング， (c)Lanczosによるスムージング， (d)Raised  cOSlneによるスムージングを示す.

1.5  1.5 

Unsmoothed  Cesaro 

1.0  1

1 Jい

t h h u t

︑ ︑ ︑

一 一 f13  l62  一 一 一 一 一 ‑

M川MNM川νリr

1.0 

0.5  0.5 

0.0  0.0 

‑{).5 

O.

。

0.5  ‑{).5 

1.0 ・1.5 2.0  0.0  0.5  1.0  1.5 

xlπ XIπ 

.(a)__~ .  . . " .   (b) 

Fig."2.2 Several smoothings for the square wave 

2.0 

1.5  1.5 

Lanczos  Raised cosine 

1.0  0.5 

‑0.5  0.0 

1.0 

0.5 

~ N=8 ^{一ーーー ¥¥':''‑}

~ N= 8^{一 一‑} N=16一一一

N=32‑‑ N=16^一一一

N=32 ‑‑

1 0 1 5 2 0 4 5 0   0 5 1 0 1 5   0.5 

XIπ XIπ 

(c)  (d)  Fig. 2.2 Several smoothings for the squ紅ewave

2.0 

本研究では，ギブ、ス現象を小さくするために， Lanczosによるスムージングを 用いた.

2.4  非線形項の取り扱いとエイリアシング誤差

非線形項の計算の際に生じる高波数成分が，低波数成分として読み換えられ ることによって生じる誤差をエイリアシング誤差という.非線形項の計算を物 理空間で行う差分法などでは，一般にこの誤差を取り除くことができない.以 下で，スペクトノレ法における非線形項の取り扱い方とエイリアシング誤差の除 去法について説明する.

次のような二つの関数

N /2‑1 

u(x)=zdnPInt 

m=‑N /2 

ν

( x ) =   L v n e

^ln 

n=‑N/2 

(2.26) 

(2.27) 

の積を等間隔におかれた N 個の格子点上で評価することを考える. N個の格子 点上で扱えるモードは ，N/2次までであり，式(2.26)，式(2.27)で表される関数は 厳密に評価できる.ところが，この二つの関数の積

w ( x )  

= U

か ) v ( x )

(2.28)  は，最高 N 次までのモードを有するため，N/2次より大きいモードは N個の格 子では分離できず，N/2次以下の低波数のモードとして読み換えられる.これが エイリアシング誤差である ここで，波数 k=‑N/2""'‑'N/2‑1での

U (

ゆか)の正確な 級数係数九は，

12 

九=計九(ゆ ( x ) e

^‑Ikx 

x 

yd  

x 

'k

e  

E  

¥1 11 1_ノ

d

u x   l e 4  

1ハ

ν yh

2州2一︑_F

‑N叩I

N‑

a一 一

' e

n R  

/

¥ 2 i  

¥V

HH HH fa it

‑A  

川

A ν

i I  

明H l A U 2

1 u

一日ι

す人]

イ

一 ︑

︑

/E

N K

・N︑

4

= F 月 刊

4

fl i‑

‑¥ Le

﹁dJπC︑_/

‑寸j

ρl Jo N︑

a p

l一Ml一

μ

一 一 一

一 (2.29) 

= 

^{L UIIIV}n 

lII+n=k 

‑N 2，5m<N 2 

‑N 2，5n<N 2 

となる.ガラーキン法を用いると，式(2.29)から正確な級数係数_Wkを求めること ができる.しかしながら，式(2.29)の最終項の計算にはN²オーダーの演算回数を 要し，計算効率の点からあまり用いられない.これに対し，スペクトル選点法 では ，FFTを利用して U，Vのそれぞれの物理空間での値を計算し，これらの積 を計算した後に，再度 FFTを用いてフーリエ級数を求める.この場合，演算回 数はM凶fに比例する程度ですむ.しかし，この方法で、はエイリアシング誤差を 生じることになる.この際，二つの関数

山 ) ) ト = 工 ; 九 ν ん ν ル 〆

^m/

〆 ^/

^/e

^〆

^/

^e

^{1fJF川川川附f1川11I111ばHX}

山 ^j ト ⁾ = ^{N 2}

^三￡:う}

〆

^/n/

^〆

^{/ffJ附l川nばeH7門円mX勺j} 

(2.30) 

(2.31 ) 

n=‑N 2 

x ̲ 2Trj 

I  ‑N  ^j = 0，1，' "，N‑1  (2.32) 

の積

︑ ︑

EEE

︐ /

︐ ︐

J X 

︐ ︐

SEE‑︑N ︑_{t W}

J

︐J X ''SEE‑¥ N U 

‑ ‑

︑EEE︐J

︐ ︐

d x 

s a a ‑ ‑ ‑ ︐ ︐ _{・ ︑}

︑N W  (2.33) 

の級数係数W_kは，

玖=三山)

₎₌₀ 

~N

^(Xj 

~

^‑/h; 

ρしV以

L J

︹

ハν

斗

A U 2 2‑ p‑

N 

N.

4F 

ト

﹁ 'd

NN可4戸1一N (2.34) 

L U/)/Vn +  L UIIIVn 

lII+n=k  lII+tl=ki:.N 

‑N 2^，5III<N  2  ‑N 2^，5III<N  2 

‑N 2，5n<N 2  ‑N 2，5n<N 2 

となる.式(2.34)の右辺第2項がエイリアシング誤差である.このエイリアシン グ誤差は精度低下や不安定性の原因になるため，何らかの除去法が必要になる. 本研究では，この誤差を除去するため， 2/3則及び 3/2則というこつのエイリア

シング誤差の除去法を用いた.

2/3員Ijは分割数Nで，U， Vの積を計 算 した後，Ikl>N/3の波数に対する係数を 0 に換えることにより ，高周波成分を除去する方法である.これに対し， 3/2則 で は 物 理 空間で U，ν を求 め る 前に，M?̲1.5Nと し て ‑ M/2三k<‑N/2， N/2三k<M/2の波数に対する係 数 に 0を補充してスベク トノレモード数をNか ら M に増やし，エイリアシング誤差を除去する方法である.

2.5  時間前進法

スペクトノレ法を編微分方程式へ適用する場合，空間的な離散化はスベクトル 法 に よ り 行 わ れ る が ， 時 間 的 な 離 散 化 は 通 常 の 有 限 差 分 法 に よ り 行 わ れ る . 本 研究では，時間積分に 4次精度ノレンゲ・クッタ・ギノレ法を用いた.

ノレンゲ・クッタ・ギル法[113] 一般的な連立微分方程式は，

4 ( ; j y t )   =  F ;   (~J h  J  .  .  .  J  ， t J 

i = 

t 2

J・・・ (2.35) 

のように示すこ とができる.

時 間 ち に お け る 値 ん か ら 時 間 らlにおける値/;J+Iを求める方法は，次のように 表せる.

1⁵¹ Stage 

k

l o =M( ん ん )

^(2.36) 

︑ . ︐ ノ

ハU

0 4  

ペノ﹄ AU   LK

ル/d・︑︑l一2

r I   一 一 (2.37) 

( 1

メ = ん 十 円 )  l (2.38) 

qil =仏o

+3~1 ーか0

2^nd Stage 

(2.39) 

ktl 

=~吋λ勺(1)

^{J. • . )}  

' i

，

 

=  1 ( 甘い ^{i l )}

(2.40)  (2.41 ) 

メ

(2)

̲ メ

(1)十円2 (2.42) 

q，2 = qd +什尚叶3玖

h

片ペら

叶

V^2一 k'2 =色t

爪 f ? ， ¥ 〕 f ; r 1

吟

¥ 1 ， ケ

ら

= ( 1 引い

(2.43) 

(2.44) 

(2.45) 

メ

(3)

̲ メ

(2)+ら

q

  3 l

= q'2 + 

3

^ら

‑ ( 1 す ) ん

Final Stage 

kl3=M(f)pλ(3) ^J' ..) 

(2.46)  (2.47) 

(2.48) 

ペ4

=~(九-

^2q

^{l 3 )}  

(4)  ./'(3) 

メ = メ + ら

(2.49) 

(2.50) 

今︑J

v κ  

l一2

4U守r. 1 

﹁3

+ 

今J

Q4  

鍋A1 一 一

0 A   (2.51 ) 

qiOの初期値はすべて 0とする.その後のq，oは一つ前のステップのqi4を置き 換えて使用する.

計算条件 3. 

燃焼現象には発 熱による大きな温度及び密度変化が伴うために，燃 焼 場 を 支 配する基礎方程式が非常に複 雑となる.ここでは， 厳密な基礎方程式から，燃 焼場においての直接 数値 計算に適用する基礎方程 式の導出を行う .

( 3 . 1 )  

( 3 . 2 )  

(3.4) 

( 3 . 3 )  

3.1 

支配方程式

燃焼場を支配する厳密な支配方程式は次のように与えられる

[ 1 1 2 ] .

連 続の式

与 +

^V . 

^{[ p u ]   =} 

⁰ 

運動量保存式

弓

^{+ρu.Vu=‑V.P+ρZCl}N 

^メ

エネノレギ一保存式(内部エネノレギー) θU  ~~T ~ ~ /~ ^¥ 

. . f ‑

ケ

+ρu.V U  = ‑V .q‑P: (Vu)+ρZCI

メ

κ 化学種保存式

。 ^c

て!̲+ u . V 

c

， = ‑~ V . 

( 〆

Iκ)+^斗 L

O[ ρ ρ   ⁱ 

= 

，12，・・.，N

(3.5) 

TB ︑︐ ︐

ju 

v 

/at︑+ ︑l

v 

u  ︐ /

rL HH

いμ 

Ft  

寸Ill1111﹂︑

︑ . . .  

v 

u  ︐ ︐ ︐

〆' a

・ ・ ・ ︑

︑¥・1 11 11 1J

ノ

' κ

/ 32一 μ   

FI l‑

‑¥  

+ P ﹁Illi‑‑L

P  一 一

ここで，

(3.6)  q=‑λVT+

心 ^{c ， V ;}

^+R

吃土空と 肱

^‑Vj)+qR

同 ，=1

= } へ

M，D^Ij) 

VX， 

=  ~(三子J(vj ^{‑ V ; ) +}  

(3.7)  i = 1.2 ・・・ N 

+ 

~[(コ1そ一号JJ芋

(3.8)  p=

バ 叫

16 

である.

u= む C 1 ‑ 1

hI=h f+ l J c p l T 

X ー_一 C，/M，

之 や

^j/M.J

j=1 

(3.9) 

i = 1.2 ・・・ N  (3.10) 

i = 1，2，・"，N (3.11) 

上記の支配方程式の直接数値計算を行うには，多大な計算時間と記憶容量が 必要となるため，乱流燃焼の直接数値計算を目的とする場合は，し¥くつかの仮 定をおく必要がある.本研究では，以下のような仮定をおく.

( 1 ) 外力なし

メ = 

0  i 

= 

1，2，・・，N  (3.12) 

(  2 )  

ソレ

( S o r e t )

効果無視

会 η ( [ ^{J ( 守 一 号 ] J 子 。}

^(3.13) 

(  3  ) 

デ、ュフォ

( D u f o u r )

効果無視

咋会(世 ^{} V ; ‑ κ ) =}  

⁰ 

(4 )圧力勾配拡散無視

( C

， ‑

x， ) 豆 =

0  (3.15)  P 

(  5 

)体積粘性無視

k = 0  (3.16) 

(  6 

)轄射熱流束無視

q R 

=  0 

(3.17) 

これらの仮定により ，上記の支配方程式は以下のように変形できる.

主 +

^V. [pu] 

⁼ 

^{0 川}

ρ

学

^{+ρu.Vu=‑V.p (3.19)} 

dt 

ρ一θU +ρ

u .

'1

u=‑

'1

. q ‑ p : (

'1

u )

(3.20)  θf 

8

:

^， + 

^{u  .} 

^{'1 C， = ‑}

^{̲ ! ̲}  

^{'1 • (ρ引 )+ι i = 1.2 ・・・ N  (3.21 )} 

θt ρ ρ  

ここで

P=[ ぺ

^μ(¥7.

^{u ) }}

^μ[('1

^{u )}  

^+ ('1

^{u y ]}  

^(3.22) 

q =‑

λ'1

T+

ρ

I

N 

h

^，C，

κ 

^(3.23) 

'1X 

= _0l  割 引

_D JI

^{( V}

^j‑

^{V ; )}  

^{i = 1，2，・・'，N (3.24)} 

である.

次に，内部エネノレギ一保存式を全エネノレギ一保存式に変形する. 先ず，運動量保存則の両辺に uを乗じると，

ρ

守

^+ρ

^{u u .}

^'1u = 

^{‑ u (}

^'1. 

^{P )}  

^(3.25) 

。(1  i  ̲( 1  i 

ρ ~I~u.ul+ ρu' '11

。

^{t¥2 )'  ¥2} ~u.u ⁾ 

1 = ‑ u (

'1

. p )  

^(3.26) 

となり，運動エネルギ一保存則を得る.この式と内部エネノレギ一保存則の和を とると，

ρ

三 (U+:u'ul+

ρ

u .

'1( 

U  +~u' ^u  1 = 

‑'1 

. q ‑ P :  

('1

u ) ‑ u (

'1 

. p )  

(3.27)  θt¥  2  )'  ¥  2  ) 

ρ

号 + 仰 向 ^，

^=‑'1

^{. q ‑}

^'1

^{. ( p . u )}  

^{問 )}

となる.連続の式より，

生!̲

+ '1 pue

，  = 

‑'1. q ‑'1. 

( p .  

u)  (3.29)  8t 

となる.ここで，

1  N  ~ 1 

ef=U+fu‑u=2JhICl‑L+

士

u.u  (3.30) 

L  ，⁼¹ ρ L  

は，全エネルギーを表している.また，運動量保存式は連続の式より，

手 + ' 1 レ

u u ] = 

‑'1 • P  (3.31)  となる.

18