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

3.3 初期乱れ場の発生方法

初期の乱れ場は，以下の手順で発生させた • x‑yの 2次元実空間において， xシ.

方向 (XI=X，X2=y)の一辺の長さを lとする場合，この実空間をフーリエ変換によって波数空間に変換したときの波数仁，

k

_mは，

21m

k J 1 = 7 ( n = 0 3 1 p p N ) ( 3 6 5 )

k^nJ

=

^一

~竺

とおける.ここで N. Mはフーリエ級数の打ち切り値で，十分大きい値とする.

この場合，ある座標X，Yでの毛方向の速度Uj(X，y)は以下の複素フーリエ級数に展開できる.

Uj(X，y)=

すエエ

^c^j⁽^kⁿ^lⁿ^k^l⁽^{いん}^y⁾ ⁽³^.⁶⁷⁾

， _{n=‑N m=‑M}

ここで，Cj(kn，k_{m )}は波数_kn，んでの複素フーリエ係数で、ある.

一方，舌しれのエネルギーは以下の式で、表される.

(m=O

，

L‑vM)

⁽³^.⁶⁶⁾

j F = H i l u j ( x p Y 1 2

^的 ^間 ⁾

式(3.68)に式(3.67)を代入すると以下のように変形できる.

;平安川之主

^C^J⁽^kⁿ^，^k^m

⁾ ^e

ⁱ⁽^川 ^，⁾

_I

』「二一

^{ー寸} ⁽³^.⁶⁹⁾

x

l 土芝 ^ξ(kA ^》

^‑^t⁽^k^p^X⁺^k^ρ)

^I ^d ^x

^砂

ここで，CJはCjの共役複素数である.式(3.69)の右辺を変形すると以下のようになる.

﹁t Il l1 11

ー 1

﹂

﹃l

i‑

‑ J

Y4 ︐d

e M

刈m

vJ 7AM

e H

刈

7 κ

p・

ムは

(3.71 )は以下のように展開できる.

ここで，

1 N M r 、

式(3.71)

=占 I

^IlCj(kn^，^k^m

} C

ん k " J J

Lι n=‑N III=‑M

Cj(kn，km

) =

^α^j⁽^k_n^，^k_m)‑ⁱ^bj (kn^，km)

ε ; ^仇

^km

) =

αj(k

λ 小

ⁱ^b^j⁽^kⁿ^，_km)

とおく. ここで， α久a_j(パ(仇丸んP比

ι

ん

k

ふ

(3.72)

(3.73)

でで、あり，どちらも実数である.式(3.72)に式(3.73)を代入すると，式(3.67)は以下のようになる.

1 ‑ τ 1 ../!‑.. ~r ^?I. ^{¥ .}^?I. J

‑u.‑=

寸エ L :

^l^a^j²⁽^kⁿ^，^k^m⁾⁺^b^j²⁽^kⁿ^，^k^m⁾^J ⁽³^.⁷⁴⁾

2J Li n=‑N HI=‑J 」

ここで，波数空間での等方性を考慮して式(3.74)を変形すると，以下のように表される.

j i 2 = j F = 会

^[^aJ

2

⁽^k^o^'^k^O⁾⁺^b^j

2

₍_k_o_'_k_O₎_]

今 M N r ， ⁽³^.⁷⁵⁾

+ 戸エエドノ ( ん

^k^m)+

ザ ( 仁

^，km)J

一方，時刻 0では非圧縮での連続の式を満足するように乱れを発生させる.

即ち，式(3.67)から，

竺、 ‑ : 1 . . : N ルf

竺L+立主主=*

L : L :

^[^kⁿ^C_t⁽

仁川 +

kmC2 (kn^，km

) ] e

l(川 y)= 0 (3.76) dx dx JL n~N

が全てのあ yで成り立つためには，

knCt(k

λJ+kλ

⁽³^.⁷⁷⁾

となり，これは波数ベクトルと速度振幅ベクトノレが直交していることを示す.

実数部，虚数部ともにこの関係、を満たすので，

kna₁(kn^，km)+ kma2 (kn^，km)

=

仁川

kn，km)+k^，

九 ( 久

^I，km

) =

0 (3.78) これから，初期の乱れ場で連続の式が成り立つには，複素フーリエ係数に関して以下の関係式を満たしていなければならないことが分かる.

はい/11)=十

^t⁽^kn， k

/ l l )

九 ( 久

I，km

) = 十 ! 仇

⁽³^.⁷⁹⁾

一方，乱れのエネノレギースペクトル分布

E ( k )

を，

時

)=uF2

会 e m ( ‑ ￡ )

⁽³^.⁸⁰⁾

とおくと，乱れのエネルギーは以下の式で、表される.

1，2 1 r^<^X^>‑^r^/^，^̲^¥^，^，^̲ ^U，2

r^∞ ^k ^{( k}

i

‑u=‑h(KM=‑ I ‑仰

‑ ‑ ; ‑ I d k

2 2^ぬ 2内丸間 _~ k_ma_x₎ ⁽³^.⁸¹⁾

ここで，kは波数い2

̲ ピ

+kfP12)，ルはエネルギースペクトルが最大となる波数，

jは乱れ強さである.また，式(3.65)，式(3.66)から

dkn ⁼dkm

与

⁽³^.⁸²⁾

であるので，k2 =kn

2 +km

2より以下の関係が成り立つ.

dk=

竺 (

^kⁿ

^~k 斗 (3.83)

[¥ k )

波数空間の等方性を考慮し，式(3.83)を用いると，式(3.81)は以下のようになる.

U'2 ( 川 ì~~1 ⁽^k^.^.+ ^k^.^.^.

i (

...lk^H2

+ k..^，2

1 1

~

^U^'²⁼

^U~-ll~

_l₁_氏似)~凶^1C

I ! ! I I

_Il

^Kⁿ^L_k⁺_{似 )}^K^m

¹ ^α ¹

⁷ ¹ ^̲

^~^K^‑ⁿ¹^̲^T K^‑^m

¹ ¹

⁽³^.⁸⁴⁾

1 k_max _J_I

式(3.75)と式(3.84)は等しいので，係数を比較して次の関係が成り立つ.

j ( んん)=ムい

^O^J^k^O⁾⁼^U ⁽³^.⁸⁵⁾

よって，式(3.75)は式(3.79)，(3.85)を用いて次式のように変形できる.

juF2=tZS

レ

^[²⁽^k^l^l^J^k^m

⁾ ⁺ ^{内丸}

⁾⁺^α2²⁽^kⁿ

^九 ⁾ ⁺ ^b

²²⁽^kⁿ

^ム

^，

J ]

= す培 2 2 [ ^旧 ^長 ^が ^い ^仏 ^い ^f ^l ^ペ ^川 ^川 ² ^刊 ^勺

⁽ ^い ^川

⁽³^.⁸⁶⁾

廿 t [ F 民 ^手 ^が ^い ^仏 ^い ² ^パ ^J ^川 ^凶 ² ^刊 ^勺 ⁽ ^い

^k^.^.^.^k

^山 ^九

式(3.84)と式(3.86)を比較すると， K2=kn2+

九

²^の関係から以下の式が成り立つ.

U'2[3₁_l_士2(k.+k..I ( k I α12(kn3km)+b12(kjfYI)= _k2 _L_kml n mlmp￨

ーァ￨

max ~ k"肌 )^A _~ ^kmax) ヲ U'2[3₁_l_士2(k..+k...I ( ^kI 022(kjnI)+b22(kjnI)= 2 nl n "~

l e x p l ーァ￨

k市f附~ k，附 )^A ^¥ ^kmax)

(3.87)

これは，連続の条件を満たす.ここで，

θ =tan - ^I (- 叫ん kll/ ~ I

⁼^t^aⁿ^‑¹⁽^̲

b ² (k/) ， kll/ ~I

一

l

a_l

( k

₁₇，

k

ll/

) ) ‑ L

^L⁴^"

l 叫ん k

ll/) (3.88) とすると，

e

は位相遅れであり O壬O壬2πの値をとる.。の値と，乱れ強さur，乱れのエネノレギースペクトノレが最大となる波数仁川を任意に定めることで式 (3.67)より乱れ場を発生させることができる.また，渦度は以下の式で定義される.

ω=

I I θ ゆ

I，X2)一伽

I (

い

1 2 1

21 みl 批2 1 ⁽³^.⁸⁹⁾

4. 乱れ強さと乱れのスケーノレの影響

これまで，予混合乱流火炎の構造についてイオンプロープによる火炎面の検知やレーザシート法による可視化などの実験，あるいは，乱れ強さや層流火炎厚さなどの混合気の特性値に基づくモデル化などの研究が盛んに行われてきた . 火炎構造の定量的な領域区分については，まだ多くの課題が残されているが，上述の手法などにより次第に明らかになりつつある.

一方，乱流燃焼についての数値シミュレーションでは，舌

L

流火炎に影響を及ぼすと考えられる支配方程式中に含まれる因子の影響を，それぞれ独立に調べることができるため，実験では困難な点、の検討を行えるものと考えられ

る.

本章では，予混合乱流伝ぱ火炎について，フーリエ・スペクトル選点法

[ 1 1 6

，

1 1 7 J

による直接数値シミュレーションを行い，乱れの火炎に及ぼす影響を調べた.

4.1 計算方法及び初期条件

数値解析には，フーリエ・スペクトル選点法を用いた.フーリエ・スペクトル選点法により非線形項を計算する場合に生じるエイリアシング誤差を除くために，

2 / 3

則

[ 1 1 0 J

を適用した.また，時間積分は

4

次精度ノレンゲ・クッタ・ギノレ法により行った.

2

次元計算領域は一辺の長さが

2mm

の正方形で，各方向の分割格子数は

1 2 8

とした.境界は，周期境界条件とした.時間刻み幅は，クーラン数が lを超えないように定めた.

初期速度は非圧縮性の連続の式を満たし，かっ，乱れのエネノレギースペクトル分布E(k)が式(4.1)となるように与えた.

( k )

E(k)=u'2 7仰￨一一

I

⁽⁴^.¹⁾

一~

^k^m^山 ^j

平均流速が零になるように乱数を用い，一様で等方的な乱れ場を与えた.

F i g u r e 4 . 1

に，丸山が

1

，

2

及び

3mm ‑

¹の場合の初期乱れの渦度を示す.

混合気の初期圧力は

O . l

:NfP

a

とし，初期温度分布には

300K

の領域の中心に

F i g

4.2に示すような円形の断熱燃焼温度の高温部を与えた.未燃混合気の質量分率分布は，温度が高くなるにつれて小さくなるように温度分布と逆の分布形状で与えた.高温部と低温部の境界の厚さは，計算による温度伝導率と実測による層流燃焼速度から求めた予熱域厚さの

4

倍とした.図には計算領域のy=l.

Omm

でのx方向の各値の変化も示している.

2.0 2.0

〉、〉、

8 1 .

o ^2.0

⁸ ¹^.0

^2.0

x mm ^X ^打¹^打¹

( a ) k

_max

= l m m ‑

¹ (b

) k m

ax:

=2mm‑

〉、

8 1 .

o 2.0

× 町1m

( c

)k_max

= 3 mm‑

Fig. 4.1 lnitial condition ofvorticity distributions

2.0

、

2.0

〉、

1.0

× 町1m

ドキュメント内予混合乱流伝ぱ火炎構造の解明と直接数値シミュレーション (ページ 34-40)

ト ー

3.3 初期乱れ場の発生方法

k

=

~竺

す エ エ

(m=O

L‑vM)

j F = H i l u j ( x p Y 1 2

;平安川之主

) e

I

』 「 二 一

x

l 土芝 ξ(kA 》

I d x

7 κ

=占 I

} C

ん k " J J

) =

ε ; 仇

) =

λ 小

ι

k

ふ

寸エ L :

j i 2 = j F = 会

2

2

+ 戸 エ エ ド ノ ( ん

ザ ( 仁

L : L :

仁川 +

) ] e

λJ+kλ

=

仁 川

九 ( 久

) =

はい/11)=十

/ l l )

九 ( 久

) = 十 ! 仇

E ( k )

時

会 e m ( ‑ ￡ )

i

‑ ‑ ; ‑ I d k

̲ ピ

与

竺 (

~k 斗 (3.83)

i (

1 1

~

U~-ll~

I ! ! I I

1 α 1

7 1 ̲

1 1

j ( んん)=ムい

レ

) + 内 丸

九 ) + b

ム

J ]

= す 培 2 2 [ 旧 長 が い 仏 い f l ペ 川 川 2 刊 勺

( い 川

廿 t [ F 民 手 が い 仏 い 2 パ J 川 凶 2 刊 勺 ( い

山 九

九

ーァ￨

l e x p l ー ァ ￨

θ =tan - I (- 叫ん kll/ ~ I

b 2 (k/) ， kll/ ~I

l

( k

k

トー

すエエ

⁾ ^e

_I

』「二一

l 土芝 ^ξ(kA ^》

^I ^d ^x

ε ; ^仇

+ 戸エエドノ ( ん

仁川

^~k 斗 (3.83)

^U~-ll~

¹ ^α ¹

⁷ ¹ ^̲

¹ ¹

⁾ ⁺ ^{内丸}

^九 ⁾ ⁺ ^b

^ム

= す培 2 2 [ ^旧 ^長 ^が ^い ^仏 ^い ^f ^l ^ペ ^川 ^川 ² ^刊 ^勺

⁽ ^い ^川

廿 t [ F 民 ^手 ^が ^い ^仏 ^い ² ^パ ^J ^川 ^凶 ² ^刊 ^勺 ⁽ ^い

^山 ^九

l e x p l ーァ￨

θ =tan - ^I (- 叫ん kll/ ~ I

b ² (k/) ， kll/ ~I

o ^2.0

^2.0