第４章　層流予混合火炎の振動燃焼の発生に関するシミュレー…

(1)

第４章層流予混合火炎の振動燃焼の発生に関するシミュレーション

４．１試作燃焼器で発生した振動燃焼の特徴と抑制のための課題ガス給湯器は集合住宅から戸建まで様々な形態の住居に設置されるため、設置場所によっては排気筒を設置する必要がある。予混合燃焼を用いた試作給湯器の排気筒設置時と排 気筒を設置しないときの騒音スペクトルの比較を図４．１（ａ）に示す。105 mm の排気筒設置時には175 Hz とその奇数倍の振動モードが観測されている。排気筒が無い場合は特定の周 波数ピークは観測されないが排気口からの低周波領域の燃焼騒音が観測される。また図４． １（ｂ）は排気筒が950 mm に延長した場合の騒音スペクトルを示している。排気筒長さが長ければ振動燃焼が必ず発生するということではなく、この場合は着火時である低負荷でファン回転数の低い場合に75 Hz の振動燃焼が観測されている。 0 10 20 30 40 50 60 70 80 90 100 0 200 400 600 800 1000 Frequency Hz

Sound pressure level dB(A)

Combustion rate = 25 kW, Fan rotational speed = 3,550 rpm, Over all = 52.04 dB(A) Combustion rate = 25 kW, Fan rotational speed = 2,940 rpm, Over all = 67.98 dB(A)

0 10 20 30 40 50 60 70 80 90 100 0 200 400 600 800 1000 Frequency Hz

Sound pressure level dB(A)

Combustion rate = 9.07 kW, Fan rotational speed = 2,000 rpm, Over all = 44.39dB(A) Combustion rate = 25 kW, Fan rotational speed = 3,600 rpm, Over all = 68.17dB(A)

Exhaust pipe length = 105 mm

Exhaust pipe length = 950 mm Exhaust pipe length = 0 mm

(a) (b)

Fig. 4.1 Noise spectrums of domestic water heater with various exhaust pipes

上記振動燃焼が発生しているときの燃焼排気ガス中NOx 濃度は 30 ppm 程度であり CO 濃度も低い良好な燃焼条件である。したがって燃焼条件の設定を変化させて排気筒設置時の振動燃焼を抑制することは難しくファン給気口を絞るなどによる対策を施している。ただし給気圧損が増大するためファン回転数を増加させることによる騒音の増加は避けられない問題が生じる。したがってこれらの課題を解決するためには振動燃焼の発生条件を把握する手法を開発するとともに発生に影響するパラメータを抽出することでバーナ形状の設計変更に反映させる必要がある。４．２振動燃焼の発生に関するシミュレーション手法の開発４．２．１音響モード解析を用いた振動燃焼の発生に関するシミュレーション複雑な構造を持つ実機燃焼器形状で決定される音響固有モードの把握が振動燃焼の発生

(2)

をシミュレーションする上でまず必要であると考えられる。（１）Rayleigh 判定式とローカルインピーダンスの位相の関係神野ら(1)_{が検討した振動燃焼の発生機構を参考に試作燃焼器を用いたシミュレーション} を行った。火炎による発熱量変動と流速変動の位相差を

φ

、ローカルインピーダンスの位相（圧力変動と流速変動との位相差）を

ϕ

としたときに燃焼器媒体への熱力学的仕事は式（４．１）のようになる。

(

)

cos

W

π

φ ϕ

ω

∝

−

（４．１）このとき振動燃焼が発生する場合の位相条件を検討すると、ローカルインピーダンスの位相によって次のように場合分けができる。

2 ϕ

= −

π

の場合（４．２）

0 < <

φ π

,

W

<

0 0,

W

0 π φ

− < <

>

2 ϕ π

=

の場合（４．３）

0 < <

φ π

,

W

>

0 0,

W

0 π φ

− < <

<

またローカルインピーダンスの極大値や極小値に火炎が存在する場合、圧力変動や流速変動がゼロとなるため振動燃焼は発生しないと考えられる。（２）試作燃焼器を用いた振動燃焼の発生シミュレーション以上の検討から前述した予混合燃焼型給湯器の音響モード解析と火炎位置の関係を検討した。このとき給湯器の形状は複雑であるため各部を単純化し各断面の代表半径と長さを 求めて音響計算を行った。この結果を図４．２に示す。

Exhaust pipe length = 105 mm, f=175Hz

0 20 40 60 80 100 120 140 160 0 0.2 0.4 0.6 0.8 1 1.2 Position m

Sound pressure level dB

0 20 40 60 80 100 120 140

Radius of water heater model mm

Sound pressure Local impedance Radius

(3)

火炎が存在する位置のローカルインピーダンスの位相は

ϕ=−π

2

であり、発熱量変動と流速変動の位相差が

− π < φ <

0

のときに振動燃焼が発生することが予想できる。上記検討では振動発生時に火炎の位置を変化させたり燃焼器形状を変更することによって振動発生の位相条件を避ける対策が可能となるが現実問題として火炎位置や形状を大幅に変更することは困難である。したがって発熱量変動と流速変動との位相差

φ

を変更させる対策が必要となると考えられる。４．２．２火炎の簡易動的特性モデルを用いた振動燃焼のシミュレーション単純な燃焼器を使用して上記音響モード解析を用いた発振予測を発展させた。バーナ火炎の発熱量変動モデルを用いて位相と増幅率を考慮してRayleigh 判定式と組み合わせることによりミキシングチャンバの長さが変化したときの振動燃焼発生の有無を検討した。（１）音響管内の音圧と粒子速度定在波が存在するときの圧力変動（音圧）と流速変動（粒子速度）の関係を導出する。閉端壁面をx=0 として、平面波が壁に向かって x が負の方向から進行しているとするとき、圧力変動および流速変動は式（４．４）で表される(2)_。

(

)

(

)

1 1

sin

m m

p

t

kx

p

u

t

kx

c

ω

ρ

=

−

=

−

（４．４）

k

c

f

c

=

ω

=

2 π

（４．５）反射波は入射波とは進行方向が反対で、x 方向に対して負の音速をもつので、上式括弧内の符号は＋となる。また、x=0 では p1=p2, u1=u2=0 となるため、式（４．６）のようになる。

(

)

(

)

2 2

sin

m m

p

t

kx

p

u

t

kx

c

ω

ρ

=

+

= −

+

（４．６）入射波と反射波が同時に存在することから、合成された圧力変動と流速変動は式（４．７）のようになる。

( ) ( )

( )

1 2 1 2

2 cos

sin

2 sin

sin

2

m m

p

kx

t

p

u

kx

t

c

ω

π

ω

ρ

=

+

=





= + =

_

−

_





（４．７）以上から定在波が存在するとき、圧力変動と流速変動は常に90 deg.の位相差を持つことがわかる。（２）火炎の簡易動的特性モデルの導入

(4)

は、振幅振動モード、音響管内の位置により符号が異なる。これをδp, δuで示す。

( )

2 cos

2 sin

p m m u

p

kx

p

kx

c

δ

ρ

=

（４．８）またこの開端付近でのδp, δuの符号と、圧力変動−発熱量変動の時間遅れτp-q’は、それぞ れ１∼３次モードの場合、表４．１のように表される。

Table 4.1 Delay time between p and q’ (τp-q’) at each mode 1st_mode（

δ

p

: ,

+

u

:

−

） 2nd mode（

δ

p

: ,

−

δ

u

:

+

） 3rd mode（

δ

p

: ,

+

δ

u

:

−

）

(

)

1 mod

,

4

p q

T

u q

T

τ

₋ _′

=

−

τ

_{′ ′}₋

3 mod

(

,

)

4

p q

T

u q

T

τ

₋ _′

=

−

τ

_{′ ′}₋

1 mod

(

,

)

4

p q

T

u q

T

τ

₋ _′

=

−

τ

_{′ ′}₋ このときτu’-q’は、第２章での検討結果より周波数によらず一定と考えられる。したがって火炎が正弦波的に式（４．９）のように変動すると仮定すると、１周期あたりの熱力学的仕事量W は、p, q’の積分（Rayleigh 判定式）で式（４．１０）のように表される。

(

)

sin

amp p q

q

′ =

A

ω φ

t

+

₋ _′ （４．９）

( )

cos

_{p q} T

W

pqdT

π

φ

ω

− ′

′

=

_Ñ

_∫

∝

（４．１０）式（４．１０）よりW が負の時は混合気の振動エネルギが発熱体に吸収され、振動は減衰する。一方、W が正の時は発熱体により仕事がなされ振動エネルギが増加し、振動が増幅され、振動燃焼の発生に至る。また式（４．１０）より熱力学的仕事量の大きさは固有振動数に反比例し、伝達関数の増幅率も周波数が高くなるにつれて低くなることがわかる。これは高周波の振動は励起されにくいことを示している。したがって以上のことから振動燃焼の発生には位相差と増幅率の双方が寄与していることがわかる。（３）解析結果と熱力学的仕事量に関する考察上記解析を基に本実験結果で検証を行った。炎口径 2 mm の多炎口金属バーナを使用し てミキシングチャンバ長を変化させていったところ、図４．３に示すように、①１次モード （140 mm~230 mm）と２次モード（350 mm~470 mm）で発生、②１次と２次モードの移行が起こり、250 mm~250 mm の範囲で振動が停止することがわかる。 図４．４に上記モデルを用いた解析結果を示す。このとき火炎の発熱量変動振幅は第２章 の解析結果より周波数に反比例すると仮定した。その結果、W>0 となる長さは１次モードで150 mm~240 mm、２次モードで 360 mm~490 mm となり、実験結果と良く一致する事がわかった。また各モード間におけるW の大きさの比較を行うと、どのモードで振動燃焼が発生するか予測可能になることがわかる。ところが発生予測と振動燃焼のレベルを精度 W とともに損失エネルギも考慮しなければならない。

(5)

7 5 φ 4-φ8 P.C.D. 23-φ1.35 ピッチ20 65 140 5 500 40 20 φ 35 77 97 110 φ φ P.C.D.

Length of mixing chamber

200 300 400 500 600 700 800 900 1000 100 150 200 250 300 350 400 450 500

Length of mixing chamber mm

Oscillation Frequency Hz Oscillation First mode Second mode Third mode

Fig. 4.3 Combustion oscillation observed on model combustor

-0.005 0 0.005 100 150 200 250 300 350 400 450 500 W1(1st mode) W2(2nd mode) W3(3rd mode)

Length of mixing chamber mm

Oscillation Oscillation

W

(6)

４．２．３火炎による圧力変動増幅を考慮した振動燃焼のシミュレーション前項で検討した熱力学的仕事量と損失エネルギの比較により振動燃焼の発生が予測できると考えた。ここではバーナを通過する圧力変動の減衰量のみを損失エネルギとして考慮し、系外への散逸や燃焼器壁面での粘性散逸は考慮していない。振動燃焼の発生経路の概 念を図４．５に示す。供給管内における混合気流量のわずかな乱れは火炎を変動させ、これ が発熱量変動となる。発熱量の変動は、周囲流体の体積を膨張・収縮させ圧力変動を引き起こす。ここで、発熱量変動が圧力変動と同期するときに正の振動エネルギを生成し、周囲流体へ仕事を行う。この仕事は、圧力変動の形で供給系にフィードバックして混合気流量をさらに変動させ、さらに発熱量の変動を増幅させる(3)_{。このように燃焼と供給系の間で} サイクルを形成し、振動発生初期におけるわずかな擾乱を増大させていく。 Attenuation Expansion f q′→ p′ Rayleigh’s criterion

( )

0 f W p q d t W ω ′ ′ = ⋅ ⋅ >

Ñ

∫

Dispersion Combustion chamber Short length f p′ Resonator p′ Supply system p′→m&′ m&′→q′ Flames Burner mass fraction Supply tube . m const& =

Fig. 4.5 Conceptual mechanism of feed-back loop

このようにRayleigh の発振条件が正となることが振動発生の必要条件となるが、この仕事が圧力変動の形で供給系に伝わるかどうかは、バーナでの圧力変動の減衰との大小で決 まると考えられる。すなわち図４．６に示すように火炎で発生する圧力変動がバーナを通過 する際に完全に減衰しないで供給系に伝わるときに供給系の圧力変動が増幅されると考えた。したがって燃焼器にバーナ火炎を組み込む前にあらかじめRayleigh の判定式による検討に加えて火炎による圧力変動の増幅率とバーナでの圧力変動の減衰率との比較を行うことによって、そのバーナが燃焼器に組み込まれたときの振動燃焼が発生する可能性のある燃焼条件の予測を行えると考えた。そこで、以下に示す実験により、Rayleigh 判定式による評価と火炎による圧力変動の増幅率の実験的な検討を行った。

(7)

Short length Resonator Combustion chamber Burner Supply tube

( )

0 f W p q d t W ω ′ ′ = ⋅ ⋅ >

∫

Ñ

Rayleigh’s criterion f p′ p′ f Oscillation p′ > pressure attenuation f Without oscillation p′ < pressure attenuation f p′ : sound pressure generated by flames f p′ : sound pressure generated by flames Attenuation

Fig. 4.6 Attenuation of sound pressure at the burner

（１）火炎による圧力変動振幅の増幅 図４．７に示すような燃焼器を用いて燃焼時および非燃焼時のバーナ前後の圧力変動と発 熱量変動を測定した。バーナ上流に設置したスピーカを用いて燃料ガス・空気混合気へ強制的に流速変動を与え、その加振周波数を変化させた。このとき混合気を加振するスピーカは150∼1500 Hz までほぼ平坦な周波数特性を持つため、印加電圧を一定にした正弦波信号を与えた。燃料には都市ガス13A (CH4: 88.5 %, C2H4: 4.6 %, C3H8: 5.4 %, C4H10: 1.5 %）を用いた。燃料ガスと空気はあらかじめベンチュリミキサによって十分に混合され、供給管壁面の小孔(直径 4 mm)から供給される。本実験で使用したバーナは厚さ 3 mm のステンレス板に、直径1 mm の炎口を 3 mm のピッチで配置した多炎口バーナである。また保炎のために各炎口にステップを設けてある。圧力変動の測定にはプローブ型マイクロフォン（B&K 製）を用いて供給管内のバーナから上流2 mm の壁面付近と、バーナ下流 7 mm のバーナ端近傍で行った。マイクロフォンのプローブ長は100 mm である。また発熱量および発熱量変動と強い相関を持つ OH ラジカルの発光強度(4)_{は中心波長}_{306.4 nm のバンドパスフィルタを受光部に設置した光電子増} 倍管（浜松ホトニクス製）で測定した。また本システムは図中のトリガによりスピーカの駆動開始、圧力変動と発光強度の測定を同期させることができる。出力信号はＦＦＴアナライザ（小野測器製CF-6400）に取り込まれ時間波形の位相差の解析が可能である。本測定において、スピーカで加振したときの圧力振幅は6 Pa 程度である。振動燃焼が生じたときの燃焼器内圧力振幅はこれの10 倍以上あるため、微小な圧力変動に対する火炎の周波数特性を測定していることになる。

(8)

220

mm

Photo-multiplier with a Band-pass filter (OH: 306.4 nm) Probe type microphone mixture mixture loud speaker φ65 Probe type microphone 3 mm φ1 φ2 Function Generator with Sinusoidal Wave Audio amp. A/D Converter FFT Analyzer Amp. Trigger Amp. Amp. PC p’, q’ p’₁ p’₂ p’1 p’2

Fig. 4.7Experimental burner set-up

スピーカにより混合気に強制的に正弦波の変動を与えたときのバーナ上流、下流の圧力変動p’1、p’2とOH ラジカルの発光強度の時間波形を図４．８に示す。このときのスピーカによ る加振周波数は550 Hz である。非燃焼時では圧力波形に変化が現れず、ほぼ同位相であったが、燃焼時のp’2の波形は、スピーカと火炎から発生する圧力変動が重なることで波形が歪み、また位相に遅れが生じる。圧力変動の振幅はスピーカにより与えたp’1がバーナ通過により減衰した結果の圧力変動と火炎変動に伴い発生する圧力変動の和がp’2として観測されると考えられる。発光強度変動は光電子増倍管の出力電圧波形を示している。この波形はノイズを多く含むためＦＦＴにてローパスフィルタをかけた後に位相差を解析した。OH ラジカルの発光強度波形は正弦波に近く圧力変動と同様にp’1に対して位相に遅れを生じていることがわかる。 本実験では図４．９に示すように燃焼時のスピーカによる圧力波と火炎で発生した圧力波 との合成波をp’fと定義し、このときの発熱量変動を示すと考えられるOH ラジカルの発光強度変動をq’とした。

(9)

Combustion rate = 1.74 kW Excess air ratio = 1.1

Excitation Frequency = 550 Hz OH radical (Luminescence with PM) Sound pressure Sound pressure 0 2 4 6 8 Time msec Phase difference Phase difference 1 p′ 2 p′

Fig. 4.8 Measured pressure oscillations and luminescence of OH radical

f

p

′

p

2

′

1

p

′

(a) With flames (b) Without flames

1

p

′

Speaker Speaker

Fig. 4.9 Definition of the sound pressure generated by the flame

p’fはp’1に対して位相差φを持つとすると式（４．１２）のように表せる。またq’も p’1 に対して位相差ψを持つと考えられるため式（４．１３）のように表せる。

(

)

cos

f f

p

′ =

A

ω φ

t

+

（４．１２）

(

)

cos

q

′ =

A

ω ϕ

t

+

（４．１３）上式にRayleigh の判定式を適用すると式（４．１４）のように変形され、|φ-ψ|< 90 deg. のときに熱力学的仕事が正となる。

(10)

( )

cos

(

)

f

W

p q d

ω

t

π

φ ϕ

ω

′ ′

=

_Ñ

∫

⋅

∝

−

（４．１４） このとき一例として図４．１０に燃焼量1.74 kW、空気比 1.1 の火炎近傍の圧力変動 p’fと発光強度変動q’の p’1に対する位相差を示す。0∼800 Hz までの周波数帯では両者の位相差は 90 deg.以下であり上式（４．１４）の W は正となる。図中にハッチングした 800∼1,050 Hz の周波数帯で位相差が90 deg.以上となり、それ以上の周波数帯で再び 90 deg.以内になった。したがって式（４．１４）のRayleigh 判定式からは 800∼1,050 Hz 以外の周波数帯で熱力学的仕事量W が正となり振動燃焼が発生する可能性を示している。 Ph as e de la y de g. -700 -600 -500 -400 -300 -200 -100 0 0 500 1000 1500 2000 Frequency Hz Combustion rate : 1.74 kW Air ratio : 1.1 p'1-p'f p'1-q' |φ-ψ| ≒ 75 deg. |φ-ψ| ≒ 100 deg. |φ-ψ| ≒ 40 deg.

Combustion rate = 1.74 kW, Excess air ratio = 1.1

Excitation frequency Hz

Phase shift deg.

Fig. 4.10 Phase difference of sound pressure and the heat release

またスピーカにより与えた圧力変動p’1に対して非燃焼時にバーナ通過による圧力振幅の増幅率（減衰）と火炎による圧力振幅の増幅率をそれぞれ下式のように定義した。このとき非燃焼時の増幅率は加振周波数50∼2,000 Hz の範囲で約-8 dB と一定であった。燃焼時におけるバーナの圧力損失は燃焼ガスの粘度や温度上昇により非燃焼時よりも増加するが、バーナ炎口を圧力変動が通過する際の振幅の減衰は音響インピーダンスZ（=ρc）に依存すると考えた。このとき燃焼によるバーナ温度の上昇により密度ρが減少、音速c は増加するため音響インピーダンスZ は非燃焼時との変化が小さいと考えられる。

( )

2

( )

1

20log

p

8 gain dB

dB

p



′



=

_

_′

_

= −





（４．１５）

( )

1

20log

p

f

gain dB

p

′





=

_

_′

_





（４．１６）

(11)

式（４．１６）を用いた圧力変動の増幅率の周波数特性を図４．１１に示す。燃焼条件は図４． ９と同じく燃焼量1.74 kW、空気比 1.0 と 1.1 である。ハッチングされた 800∼1,050 Hz 以外の周波数領域でRayleigh 判定式による熱力学的仕事量 W は正となり、かつ空気比 1.0 では100 Hz と 500 Hz 付近で増幅率が 0 dB 以上、空気比 1.1 では 400∼750 Hz にわたって増幅率が0 dB 以上となっていることがわかる。したがって本燃焼条件では上記周波数領域でRayleigh 判定式による W が正となっているものの、バーナでの増幅率が-8 dB であることからバーナ上流へ火炎により発生した圧力変動が伝播する際に減衰され振動燃焼が抑 制されると考えられる。また同時に図４．１１に示す空気比1.0 での特性は空気比 1.1 と異なり100 Hz と 500 Hz 付近で増幅率が 0 dB 以上となっている。しかし増幅率が低いため、この燃焼条件においてもRayleigh 判定式が正であったとしても振動燃焼は発生しないと考えられる。

0

500 1000

1500

2000

10

5

0 -5

-10

-15

-20

-25

-30

Excitation frequency Hz

Gain dB

Combustion rate = 1.74 kW, Excess air ratio = 1.0 and 1.1

λ= 1.0

λ= 1.1

λ= 1.0

λ= 1.1

Fig. 4.11 Gains of sound pressure amplification at various excess air ratios

前述した検討を表４．２に示す燃焼条件と加振周波数で圧力変動の振幅と位相の測定によ り行った。

Table 4.2 Experimental conditions

Excess air ratio 0.8∼1.7

Combustion rate kW

(kcal/h) (800∼2,300) 0.93∼2.67 Excitation frequency Hz 50∼2,000

(12)

このバーナを設置したときの振動燃焼の発生を予測、検証するための実験を行った。図４． １２（ａ）に示す軸対称燃焼器は前述の火炎の周波数特性を測定する燃焼器とは、供給管長さ、および形状が異なり、バーナプレートのみ同じものを使用している。供給管長さは460、560、 660 mm であり燃焼室として 40 mm の円筒を設置している。燃料には都市ガス 13A を使用し、燃料ガスと空気の混合気は供給管壁面の小孔(直径 4 mm)から供給される。さらに供給管壁面および燃焼室壁面には圧力測定用の小孔を約30 mm 間隔に設けており、燃焼器内部の圧力変動分布を測定できるようになっている。振動発生領域の測定は燃焼量を固定して空気流量を連続的に変化させて行った。供給管長さ460 mm のときの上記燃焼器を用いて 測定した振動領域の一例を図４．１２（ｂ）に示す。550 Hz と 900 Hz 付近の２つのモードの振 動燃焼が観測され、さらに両モードが重なる領域がある。これは例えば2.91 kW において高空気比側から徐々に空気流量を絞っていくと空気比1.1 まで 550 Hz のモードで振動が持続し、その後900Hz のモードに瞬時に移る現象を示している。また低空気比側から流量を上げると1.3 まで 900 Hz の振動が持続する。 40 mm φ65 460, 560, 660 mm mixture mixture

(a) Experimental combustor (b) Areas of observed combustion oscillation 0.6 0.8 1.0 1.2 1.4 1.6 1.8

Excess air ratio

550 Hz mode 900 Hz mode 550 Hz mode Without oscillation Oscillation 3.5 3.0 2.5 2.0 1.5 1.0 0.5 Combustion rate kW

Fig. 4.12 Combustion oscillation observed on empirical combustor

このとき図４．１３に示すように２つの周波数での振動燃焼はバーナ下流が開端、供給管下端 が閉端となっている音響モードが励起された550 Hz が２次、900 Hz が３次モードであることが圧力変動分布よりわかる。また空気比1.0 付近で高周波の振動燃焼が観測され希薄側と過濃側で低周波となっている理由として式（４．１４）における圧力変動と発熱量変動の位相差|φ-ψ|の大きさが異なることと空気比による若干のルイス数変化によって起こる火炎面の擾乱による発熱量変動強度の違いが考えられる。

(13)

0 20 40 60 80 100 0 100 200 300 400 500 600 552.5 Hz 1105 Hz 1657.5 Hz Position mm Burner So un d pr es su re P a Burner 895 Hz 1785 Hz Position mm Burner So un d pr es su re P a

air ratio = 1.3, combustion rate = 1.74 kW

air ratio = 1.0, combustion rate = 2.33 kW

0 20 40 60 80 100 0 100 200 300 400 500 600 Burner

Excess air ratio = 1.0, Combustion rate = 2.33 kW Excess air ratio = 1.3, Combustion rate = 1.74 kW

Fig. 4.13 Sound pressure distribution in the supply tube at 1st and 2nd mode

スピーカによる加振実験による測定値を式（４．１６）で整理した550 Hz の振動数にお ける圧力変動の増幅率を図４．１４に示す。このとき Rayleigh の判定式（４．１４）が正となる空気比および燃焼量におけるデータを用いている。この結果をバーナを通過する圧力変動の増幅率（減衰）である8 dB 以上の領域で区切ると２つの領域があらわれることがわ かる。これは振動燃焼の発生領域の測定結果と定性的に一致する。図４．１５に示すもう１つ のモードである 900 Hz でも同様の処理を行うと振動が増幅率の高い領域で発生していることがわかる。同様に供給管長さを560 mm（図４．１６、図４．１７）、660 mm（図４．１８、 図４．１９）としたときの振動燃焼の測定を行った。両者とも２つの振動モードが観測され、 増幅率による評価により振動燃焼の発生領域の予測が可能なことがわかる。以上の実験よってRayleigh の判定式が正となる場合でも火炎による圧力変動振幅がバーナによる減衰量を上回らないと振動燃焼が発生しないことを定性的に示した。このことは実機においてバーナの圧損が高いと振動が発生しにくくなる現象を説明できると考えられる。

(14)

10 0 8 6 0.6 0.8 1 1.2 1.4 1.6 1.8 0.6 0.8 1 1.2 1.4 1.6 1.8 0.5 1.0 1.5 2.0 2.5 3.0 3.5

550Hz

no oscillation

oscillation

Co

mb

us

ti

on

r

at

e

kW

Gain : 8 dB Gain : 0 dB

550Hz mode

900Hz mode

550Hz mode

no oscillation

Excess air ratio Excess air ratio

Combustion rate kW

(a) Distribution of the sound pressure gain (b) Areas of observed combustion oscillation

Oscillation

Without oscillation

Fig. 4.14 Evaluation of occurrence of combustion oscillations with the gain of sound pressure (Supply tube length = 460 mm, 550 Hz mode)

0.6 0.8 1 1.2 1.4 1.6 1.8 -5 0 -10 15

900Hz

no oscillation

oscillation

Gain : 8 dB Gain : 0 dB 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -10

Co

mb

us

ti

on

r

at

e

kW

0.6 0.8 1 1.2 1.4 1.6 1.8

550Hz mode

900Hz mode

550Hz mode

no oscillation

8 -10 Combustion rate kW

Oscillation

Without oscillation Without oscillation

(15)

0.6 0.8 1 1.2 1.4 1.6 1.8 0.6 0.8 1 1.2 1.4 1.6 1.8 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Co

mb

us

ti

on

r

at

e

kW

Air ratio

0 0 2 4 4 4 4 6 6 8 1012 10

450 Hz

no oscillation

oscillation

8

no oscillation

450 Hz mode

750 Hz mode

450 Hz mode

Combustion rate kW

Oscillation

Without oscillation

0.6 0.8 1 1.2 1.4 1.6 1.8

Co

mb

us

ti

on

r

at

e

kW

Air ratio

0.5 1.0 1.5 2.0 2.5 3.0 3.5 -10 -4 -2 0 2 6 8 12

750 Hz

oscillation

no oscillation

0.6 0.8 1 1.2 1.4 1.6 1.8

no oscillation

450 Hz mode

750 Hz mode

450 Hz mode

Combustion rate kW

Oscillation

Without oscillation Without oscillation

(16)

0.6 0.8 1 _1.2 _1.4 1.6 1.8

Co

mb

us

ti

on

r

at

e

kW

385 Hz mode

635 Hz mode

no oscillation

0 ₂ 2 ₄ 6 8 10 0.5 1 1.5 2 2.5 3 3.5 0.6 0.8 1 1.2 1.4 1.6 1.8

385 Hz

Gain : 8 dB_{Gain : 0 dB} 8

no oscillation

Combustion rate kW

Excess air ratio _{Excess air ratio}

Oscillation Without oscillation

Without oscillation

0.6 0.8 1 _1.2 _1.4 1.6 1.8

Co

mb

us

ti

on

r

at

e

kW

385 Hz mode

635 Hz mode

no oscillation

0.5 1 1.5 2 2.5 3 3.5 0.6 0.8 1 1.2 1.4 1.6 1.8 0 2 4 10 14 16 ₁₈ 6 8

635 Hz

Gain : 8 dB_{Gain : 0 dB}

no oscillation

Combustion rate kW

Oscillation

Without oscillation

(17)

４．２．４有限要素法を用いた振動燃焼のシミュレーション上記の予測手法は計算が簡易である長所があるが一方で複雑な燃焼器を対象とする場合例えば音響モードの解析が困難である。また火炎の発熱量変動モデルを導入した計算が可能になったが振動抑制対策へ結びつけることが難しいと考えられる。そこで複雑な燃焼器形状に対応させるため、有限要素法を用いた振動燃焼解析用プログラムを開発した。（１）発熱量変動を考慮した波動方程式と簡易動的特性モデルの導入熱流体を支配する以下の保存式（４．１７）∼（４．１９）と状態方程式（４．２０）から波動方程式(5)_{（４．２１）が導出できる。このとき①振動時の変動量の大きさは平均量} に比べて小さい、②平均流は音速に比べて小さい、③平均流による乱れや熱の輸送現象を無視する、という仮定を用いている。

0 D

Dt

ρ ρ

+ ∇ ⋅ =

u

（４．１７）

0 D

Dt

ρ

u

+ ∇ =

P

（４．１８）

DT

Cv

p

Q

Dt

ρ

+ ∇ ⋅ =

u

&

（４．１９）

P

=

ρ

RT

（４．２０）

(

D Dt

=

∂ ∂

t

+ ⋅∇

u

)

2 2 2 2

p

R

Q

c

p

t

Cv

t

∂

− ∇ =

∂

&

（４．２１）式（４．２１）において生成項中の

&q

は単位体積あたりの発熱量変動であり、本研究では圧力変動と周波数の関係を式（４．２２）で表される関数と仮定した。ここで、

A

( )

ω

,

p

は発熱量変動の振幅、θは圧力変動と発熱量変動との位相差、

q

は平均発熱量（燃焼量）である。

(

)

{

}

( , ) exp

,

Q

&

=

A

ω

p

⋅

_



i

ω θ ω

t

−

p



_

+

q

（４．２２）さらに、本研究では、周波数領域での数値計算を行うため、圧力変動も単純な周期関数（例えば正弦波関数）と仮定すると、以下のようにヘルムホルツ型の波動方程式に変換できる。 2 2 2 2

i

R

p

Q

c

Cv

ω

∇ +

= −

⋅

&

_{（４．２３）} このとき、右辺の

&q

は、式（４．１８）のように、圧力変動に対する発熱量変動の伝達関数

( )

TF

ω

,

p

によって表されると仮定する。すなわち式（４．２４）は燃焼領域に圧力変動

p

( )

ω

が作用したときの圧力変動と発熱量変動との関係を示している。

(18)

本研究では圧力変動

p

( )

ω

をバーナ炎口上での値と定義した。このとき工業用ボイラなどに使用される乱流火炎などを例にとると、圧力変動を定義した場所から火炎の局所までの距離がそれぞれ異なるため、圧力変動と発熱量との関係はそれぞれ異なる特性を示すと考えられる。すなわち異なる伝達関数

TF

( )

ω

,

p

が火炎面全体に分布することになる。しかし局所の伝達関数を測定することは非常に困難であるため本研究では発熱量は分布を持たず火炎が存在する領域に均一に分布し、さらにその特性が１つの伝達関数で表されると仮定した。したがって式（４．２４）の伝達関数

TF

( )

ω

,

p

は、単位体積の火炎領域に作用する圧力変動と、位相差を含む発熱量変動との関係を示していると考えることができる。本研究では以上の仮定のもとに式（４．２３）を解析するために圧力変動に対する発熱量変動の伝達関数

TF

( )

ω

,

p

を第２章に述べたの火炎可視化システムと画像処理手法により 測定した。図４．２０に圧力変動に対する発熱量変動の位相差と発熱量変動の振幅と強制的に 燃料ガス・空気混合気に与えた圧力変動振幅との比の周波数特性を示す。図４．２０（ａ）の破線が示すように位相差は300 Hz 程度まで線形に大きくなることがわかる。また圧力変動振幅に対する発熱量変動の振幅は加振周波数の-1.85 乗に比例することがわかる。以上の測定結果を式（４．２４）の型式に変換し数値計算における火炎の伝達関数として与えた。 1 10 100 1000 10 100 1000 Excitation frequency Hz

Amplitude of the transfer function W/Pa

measured Predicted -200 -150 -100 -50 0 50 100 150 200 0 100 200 300 400 Excitation frequency Hz

Phase shift deg.

measured predicted

(a) Phase difference of the heat release to sound pressure (b) Amplitude of the transfer function of heat release Combustion rate = 1.163 kW, Excess air ratio = 1.0

Fig. 4.20 Phase difference and amplitude for the transfer function to sound pressure

（２）有限要素法を用いるための線形化と境界条件振動燃焼においては、火炎による発熱量変動は音響場に作用し、また逆に音響場によって、発熱量変動が大きな影響を受ける。さらに発熱量変動は圧力変動に対して非線形な関係にある。したがって式（４．２３）の生成項に含まれる発熱量変動が圧力変動に対して任意の伝達関数で表される場合にも解析可能な方法が必要となる。一般の音響解析において広く利用されている有限要素法(6)_{は分割された要素に関して、そ} の節点間の圧力を比較的単純な内挿関数（たとえば１次関数）を用いることで要素内の圧力変動を近似し、系全体のポテンシャルエネルギ

L

を求め、これが最小または停留するよ

(19)

うに圧力変動が分布したときに、その系が平衡状態に達するという変分原理を利用したものである。このとき、式（４．２５）に示すように、系全体のインピーダンスマトリクスを

Z

( )

ω

とすると、圧力変動の分布

p

と

Z

( )

ω

の積がポテンシャルエネルギに相当する。また、

g p,

( )

ω

は、火炎における発熱量変動に伴い生成するエネルギを示しており、

h

は壁面などに与えられる境界条件のエネルギを示す。したがって、式（４．２５）の左辺が停留するように圧力変動

p

が分布したときに、式（４．１７）を満足する結果が得られる。

( )

ω

= =

(

p

,

ω

)

+

Z

p

f

g

h

（４．２５）通常の音響場のみの解析では火炎が存在しないので、生成項

g p,

( )

ω

を含まず

( )

Z

ω = =

p

f

h

となり、逆行列

Z

( )

ω

−1を求めれば

p

を解析できるが、本解析では圧力の関数である非線形の生成項

g p,

( )

ω

を含む式（４．２５）を求めるために繰り返し計算を行う 必要がある。このときの解析フローを図４．２１に示す。 Read Data Assemble Decompose Z( )ω

Set up Initial Solution Form Load Vector

(p,ω)

f = g + h

Solve the System_{( )} ₍ ₎

, p ω ω Z p = g + h Convergence ? Interpret Results Next Frequency ? End No Yes No Yes ( )ω Z

Fig. 4.21 Analytical flow chart of numerical simulation

（３）モデル燃焼器を用いた振動燃焼の発生に関するシミュレーション

図４．２２に示す２種類の長さ（208, 408 mm）の燃焼室を持つ軸対称燃焼器での振動燃焼

現象を上記の伝達関数を使用して数値解析を行い、実験との比較を行った。このとき、予備実験において、燃焼室長さを徐々に長くしていくと、完全に振動燃焼が停止している状

(20)

態から、ある燃焼室長さを境に振動燃焼が発生し、さらに燃焼室を長くすると振動燃焼が激しくなり火炎が吹き飛ぶ現象が観察された。そこで燃焼室長さは、安定に振動燃焼が発生する状態と、完全に停止する状態の代表として上記の２種類の長さを選んだ。 Combustion room L mm L = 208∼408 mm φ100 φ65 φ12 450 Water cooling

Natural gas and air mixture

Infinite boundary condition

Initial condition

p=0.1 mPa Flame region

Temperature distribution Temp. = 288 K

Fig. 4.22 Axisymmetric combustor and analytical model of combustor

式（４．２３）における生成項には式（４．２４）の伝達関数を実験において火炎が存在する領域とほぼ同じ領域においてのみ与えた。燃焼器内の温度分布は、燃焼室中心軸上での熱電対による測定値をもとに与えたが、火炎近傍では輻射による測定誤差が大きかったため、断熱火炎温度になるように補正して与えた。また、実験ではバーナリムを水冷しているため、バーナリムより上流への熱伝導は少ないと仮定して、ミキシングチャンバ内の温度は288 K に設定した。また、燃焼器を囲むように、燃焼器端面から上下に200 mm、半径方向に 100 mm の外部計算領域を設け、その境界に無限境界条件である放射インピーダンスZ（=ρc）を定義することで、燃焼器からの音響エネルギの散逸を考慮した。また、ミキシングチャンバ下部に微少な圧力変動（0.1 mPa）を与えることにより、混合気がミキシングチャンバに流入する際の圧力変動を模擬している。このとき、ミキシングチャンバ全域に圧力変動を与えた場合の計算も行ったが、計算結果への影響はないことを確認している。本解析モデルは、要素数3181、節点数 3497、解析のための CPU 時間は 3.5 min であった。 図４．２３に燃焼器出口付近の騒音測定結果と計算結果を示す。実験において振動燃焼が発 生した燃焼室長さが408 mm の場合、周波数に関しては、良い一致を示したものの、音圧レベルの一致までは得られず、音圧のピークレベルを比較すると計算結果の方が大きくなった。計算では、管壁で音は全反射すると仮定しており、管壁での減衰や管壁を振動させることによる管壁からの音の放射を考慮していない。また、開口部のエッジでは流速変動

(21)

の減衰が起こり、このことにより音圧が減少することも考えられる。これらの影響が計算結果と実験結果のピークレベルの違いの原因と考えられる。したがって、より精度の高い数値解析を行うには、管壁における音圧の減衰や透過および燃焼器開口部での音圧の減衰 を考慮する必要がある。図４．２４に燃焼器内の音圧分布の計算結果と実測値を示す。燃焼室 の音圧レベルがミキシングチャンバよりも大きいことから燃焼室を共鳴器とする振動燃焼が発生していることがわかる。 L = 208 mm _{(L = 208 mm)}Experiment L = 408 mm Calculation (L=208, 408 mm) Experiment 0 100 200 300 400 500 Frequency Hz 3.0 2.5 2.0 1.5 1.0 0.5 0 Sound pressure Pa

(22)

Fig. 4.24 Distribution of sound pressure at the center of the combustor 一方、実験で振動燃焼が発生しない燃焼室長さが208 mm の場合では、計算においても L=408 mm の場合のような音圧のピークは現れなかった。今回の計算では、伝達関数を 400 Hz までしか測定していないため、400 Hz 以上についての計算を行っていない。このため、さらに高い周波数に音圧のピークが現れている可能性がある。そこで、実験により求めた伝達関数を1000 Hz まで外挿して与え、計算を行った。また同時に燃焼室長さを変化させ た計算も行ったので、その結果をあわせて図４．２５に示す。燃焼室長さを短くしていくと、 ピークとなる周波数は高周波側に移動するが、ピークレベルもそれにつれて小さくなり、 L=208 mm では L=408 mm と比較するとそのレベルは極端に小さく、事実上振動燃焼は発生していないと考えてよい。このように、本計算では、ピークレベルを比較することで振動燃焼の有無を予測できることが明らかとなった。

(23)

200 300 400 600 800 1000 3.0 2.5 2.0 1.5 1.0 0.5 0 Sound pressure Pa Frequency Hz

Fig. 4.25 Effect of the length of combustion chamber

（４）実機燃焼機器を対象としたシミュレーション 図４．２６に示す低NOx 型燃焼器を対象とした振動燃焼の発生に関するシミュレーションを行った。バーナは予混合燃焼方式を用いており、その下流に熱交換器が配置されている。前列の水管群により火炎冷却を行い NOx 排出濃度を抑えるとともに、その際に発生した CO を後列水管群との間にもうけた断熱空間で完全燃焼させる構造になっている。本燃焼器で空気比 1.2、燃焼量 20.5 kW にて強い振動燃焼が発生し、そのピーク周波数は 960 Hz であった。数値計算に用いる火炎の動的特性をモデル化するために炎口面積、燃焼量および空気比の設定条件から簡易動的特性モデルを適用し、発熱量変動振幅と遅れ時間を算出した。これを式（４．２４）の伝達関数型式に変換し、設計図面から作成した計算用メッシュの火炎領域に設定した。このとき音速および密度は実測した燃焼室内温度を用いて算出した値を用いている。この結果、音圧に関するスペクトルは実測値と比較的良い一致を示してい る。また図４．２７に示すように燃焼器内の音圧分布の結果よりバーナ上流のミキシングチャ ンバが共鳴器となる振動燃焼が発生していることがわかった。以上の解析によりミキシングチャンバ形状の変更やバーナ構造の改善対策を施し、振動燃焼を抑制することができた。

(24)

0.01 0.1 0 500 1000 1500 calculation 0.01 0.1 measured

MFB(Metal Fiber Burner) Air Natural gas Mixing chamber 10 1.0 0.1 0.01 Sound pressure Pa 10 1.0 0.1 0.01 Sound pressure Pa 0 500 1000 1500 Frequency Hz

Fig. 4.26 Experimental and simulation results of combustor with heat exchanger

60 50 40 30 20 10 0 Sound pressure Pa Length m 0 200 400 600 800 1000 1200 1400 1600 1800 Pressure (Pa) 5.493 4.813 4.132 2.771 2.091 1.410 0.729 0.0488 3.452

(25)

４．３まとめ混合気流速変動に対する発熱量変動の振幅と位相関係を用いた簡易動的特性モデル用いた振動燃焼の発生シミュレーション手法の開発とその検証を行った。Rayleigh 判定式に動的特性モデルを適用してその熱力学的仕事量による振動発生モードが予測できること、そして複雑な形状を持つ燃焼器に対応できる有限要素法による振動燃焼シミュレーションプログラムの開発を行い、実機燃焼器でも振動燃焼の解析が可能であることを示した。以下に得られた知見を示す。（１）混合室長さを変化できる燃焼器において動的特性モデルと音響的なローカルインピーダンスの位相ψを Rayleigh 判定式に導入して振動燃焼が発生する混合室長さとその振動モードを熱力学的仕事量の大きさで予測できることを示した。（２）火炎による圧力変動の増幅とバーナによる減衰との関係を定量的に評価することにより Rayleigh 判定式による熱力学的仕事量が正であっても振動燃焼が発生しない場合があることを実験的に示し、経験的にバーナ圧損を増加させると振動燃焼が抑制されるという現象を説明した。（３）発熱量変動の動的特性モデルを波動方程式の生成項の適用した振動燃焼の発生シミュレーションプログラムを開発した。モデル燃焼器による実験で燃焼室長さを変化させたときの振動燃焼発生の有無を予測できることを示した。（４）上記シミュレーションプログラムは有限要素法を用いるため複雑な形状の燃焼器に対応でき、実機燃焼器にて発生した振動燃焼の発生を解析し抑制対策を施すことができた。４．４参考文献 (1) 神野, 石川, 第 28 回燃焼シンポジウム前刷集, p.533-535, 1990 (2) 福田, 奥田, 騒音対策と消音設計, p.20-23, 1967, 共立出版 (3) 岸本, 燃焼研究, vol. 101, p.15-25, 1995

(4) M.Katsuki, 21st_{International Symposium on Combustion, p.1543-1550, 1986} (5) T. Takeno, Institute of Space and Aero. Science, University of Tokyo Report, No.430,

p.309-347, 1968

第４章 層流予混合火炎の振動燃焼の発生に関するシミュレー…