(57)【要約】【課題】内燃機関の熱発生率波形を作成するための工数を低減すると共に、例えば燃焼速度について、所要の精度を確保しながら従来よりも簡易に推定し、評価できるようにする。【解決手段】混合気の着火後に熱発生率が増加する熱発生率増加期間（例えば前半燃焼期間ａ）における、この熱発生率のクランク角度変化に対する増加率を、熱発生率波形の特性値の一つである熱発生率傾きｂ／ａとして規定する。この熱発生率傾きを、熱発生率増加期間において予め設定された所定時期の燃料密度（例えば熱発生率最大時燃料密度ρfuel@dQpeak）に基づいて推定し、この推定した熱発生率傾きを用いて熱発生率波形を作成する。【選択図】図２
2015162983.tif 000002
【特許請求の範囲】
【請求項１】
火花点火式の内燃機関における熱発生率波形を算出するための装置であって、
混合気の着火後に熱発生率が増加する熱発生率増加期間における、この熱発生率のクランク角度変化に対する増加率を、前記熱発生率波形の特性値の一つである熱発生率傾きとして規定し、
前記熱発生率傾きを、前記熱発生率増加期間において予め設定された所定時期の気筒内の燃料密度に基づいて推定し、この推定した熱発生率傾きを用いて熱発生率波形を算出する構成となっていることを特徴とする内燃機関の熱発生率波形算出装置。
【請求項２】
請求項１に記載の内燃機関の熱発生率波形算出装置において、
混合気の着火時期から熱発生率が最大になる熱発生率最大時期までの期間における熱発生率の平均的な増加率を前記熱発生率傾きとし、この熱発生率傾きを前記熱発生率最大時期の燃料密度の一次関数として算出する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項３】
請求項１または２のいずれかに記載の内燃機関の熱発生率波形算出装置において、
前記熱発生率傾きを、燃焼効率に基づく補正係数を乗算して算出する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項４】
請求項１〜３のいずれか一つに記載の内燃機関の熱発生率波形算出装置において、
前記熱発生率傾きを、機関回転速度に基づく補正係数を乗算して算出する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項５】
請求項１〜４のいずれか一つに記載の内燃機関の熱発生率波形算出装置において、
前記熱発生率傾きを、熱発生率最大時期の筒内体積に基づく補正係数を乗算して算出する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項６】
請求項１〜５のいずれか一つに記載の内燃機関の熱発生率波形算出装置において、
前記熱発生率傾きを、空燃比に基づく補正係数を乗算して算出する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項７】
請求項１〜６のいずれか一つに記載の内燃機関の熱発生率波形算出装置において、
前記熱発生率傾きを、ＥＧＲ率に基づく補正係数を乗算して算出する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項８】
請求項１〜７のいずれか一つに記載の内燃機関の熱発生率波形算出装置において、
混合気の着火から燃焼終了までのクランク角度期間を底辺とし、熱発生率最大時期における熱発生率を頂点とする三角波形によって熱発生率波形を近似し、
前記三角波形において着火時期から熱発生率最大時期までの熱発生率を表す斜辺の傾きを、前記熱発生率傾きとして規定する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項９】
請求項８に記載の内燃機関の熱発生率波形算出装置において、
前記三角波形における着火時期から熱発生率最大時期までの期間が、機関負荷率、空燃比、ＥＧＲ率および油水温のいずれにも依らないものとして、前記三角波形を作成する構成となっている、内燃機関の熱発生率波形算出装置。
【請求項１０】
火花点火式の内燃機関における熱発生率波形を算出する方法であって、
混合気の着火後に熱発生率が増加する熱発生率増加期間における、この熱発生率のクランク角度変化に対する増加率を、前記熱発生率波形の特性値の一つである熱発生率傾きとして規定し、
前記熱発生率傾きを、前記熱発生率増加期間において予め設定された所定時期の気筒内の燃料密度に基づいて推定し、この推定した熱発生率傾きを用いて熱発生率波形を算出することを特徴とする内燃機関の熱発生率波形算出方法。
【請求項１１】
請求項１０に記載の内燃機関の熱発生率波形算出方法において、
混合気の着火時期から熱発生率が最大になる熱発生率最大時期までの期間における熱発生率の平均的な増加率を前記熱発生率傾きとし、この熱発生率傾きを前記熱発生率最大時期の燃料密度の一次関数として算出する、内燃機関の熱発生率波形算出方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、火花点火式の内燃機関における熱発生率波形を算出するための装置、およびその方法に関連し、特に、混合気の着火後に熱発生率が増加する際の燃焼の速度に着目して、熱発生率波形を得るための技術に係る。
【背景技術】
【０００２】
従来より、内燃機関の燃焼状態を表現するために、気筒内の熱発生率をウィーベ（Wiebe）関数によって近似することが行われている。ウィーベ関数は、複数のパラメータを特定することによって、熱発生率波形を好適に表現することができるものであって、内燃機関の燃焼による熱発生率や質量燃焼割合などを推定するために利用されている。【０００３】
例えば、特許文献１に記載されているウィーベ関数パラメータの決定方法では、熱発生率が最大となるクランク角度での燃焼割合に基づいて、所定の式によりウィーベ関数の形状パラメータｍを同定するようにしている。そして、ｋ、ａ／θpm+1、θbといった他のパラメータについてもそれぞれ所定の式により同定して、実際の熱発生パターンに高い精度で適合するように、ウィーベ関数を決定することができる。【０００４】
前記特許文献１には、そのように複数のパラメータｍ，ｋ，ａ／θpm+1，θbを同定し、ウィーベ関数を決定する作業を種々の運転条件に対して行うことにより、これらのパラメータと内燃機関の運転パラメータ（負荷率、回転速度、空燃比、点火時期など）との関係を把握することができる、と記載されている。また、そうして把握した関係を利用すれば、内燃機関のあらゆる運転条件に対してウィーベ関数を決定し、当該内燃機関の燃焼状態を精度良く表すことができると記載されている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００７−１７７６５４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら前記特許文献１には、ウィーベ関数のパラメータｍ、ｋ、ａ／θpm+1、θbと内燃機関の運転パラメータとの関係を特定する具体的な方法について開示されていない。このため、実際には殆ど全ての運転条件についてパラメータｍ、ｋ、ａ／θpm+1、θbを同定し、それぞれの運転条件についてウィーベ関数を決定しなくてはならない。つまり、前記従来の方法には、燃焼率波形を作成するための工数をさらに低減し、コストを削減する余地が残されている。【０００７】
また、前記の方法では、パラメータｍ、ｋ、ａ／θpm+1、θbをそれぞれ同定し、ウィーベ関数を決定して初めて熱発生率波形全体を表現することができ、これに基づいて燃焼状態を評価できるものである。よって、熱発生率波形全体を表現することなく、例えば初期の燃焼速度だけを簡易な方法で推定し、評価することはできない。【０００８】
本発明はかかる諸点に鑑みてなされたものであり、その目的とするところは、燃焼状態を表す指標の一つである燃焼速度に着目し、熱発生率波形を作成（算出）するための工数を低減すると共に、例えば燃焼速度について所要の精度を確保しながら従来よりも簡易に推定し、評価できるようにすることにある。
【課題を解決するための手段】
【０００９】
−発明の解決原理−
本発明の発明者は、火花点火によって着火された混合気が燃焼し、熱発生率が増加していく際に、その熱発生率の増加率（熱発生率傾き）が燃料密度と高い相関を有しており、その熱発生率傾きに対する機関負荷率および点火時期の影響が燃料密度によってまとめて表現できる、という新規な知見を得た。【００１０】
かかる新規な知見に基づいて、本発明の解決原理は、熱発生率波形の特性値の一つとして前記熱発生率傾きを用い、この熱発生率傾きを燃料密度に基づいて推定するようにしたことにある。【００１１】
−解決手段−
具体的に、本発明は、火花点火式の内燃機関における熱発生率波形を算出するための装置を対象とするもので、混合気の着火後に熱発生率が増加する熱発生率増加期間における、この熱発生率のクランク角度変化に対する増加率を、熱発生率波形の特性値の一つである熱発生率傾きとして規定する。そして、この熱発生率傾きを、前記熱発生率増加期間において予め設定された所定時期の気筒内の燃料密度に基づいて推定し、この推定した熱発生率傾きを用いて熱発生率波形を算出する構成とする。【００１２】
前記の特定事項によると、内燃機関における気筒内の混合気の燃焼をモデル化し、熱発生率波形を算出するに際して、その熱発生率波形の特性値の一つとして、混合気の着火後に熱発生率が増加する間の燃焼速度を表すパラメータ、即ち、クランク角度変化に対する熱発生率の増加率（熱発生率傾き）を用いている。この熱発生率傾きは、内燃機関の負荷率や回転速度、空燃比、点火時期、ＥＧＲ率、油水温等、種々の運転条件によって変化するものであるが、機関負荷率および点火時期の影響については、熱発生率増加期間における燃料密度という一つのパラメータによって、まとめて表現することができる。【００１３】
そこで、熱発生率増加期間において予め設定した所定時期（例えば、最大熱発生率の８０％というように予め設定した熱発生率に達する時期）の気筒内の燃料密度に基づいて、熱発生率傾きを推定するようにすれば、機関負荷率および点火時期の両方に基づいて推定するのに比べて、熱発生率傾きの推定にかかる工数を低減できる。こうして推定した熱発生率傾きを用いることで、熱発生率波形を算出する工数も低減できる。【００１４】
しかも、熱発生率波形の全体を作成することなく、前記のように燃料密度に基づいて熱発生率傾きだけを推定できるので、着火後に熱発生率が増加する間の燃焼速度について、所要の精度を確保しながら従来よりも簡易に推定し、評価することができる。【００１５】
前記の所定時期は熱発生率増加期間において設定すればよいが、熱発生率の変化と混合気の燃料密度との相関を調べたところ、特に混合気の着火時期から熱発生率最大時期までの期間における熱発生率の平均的な増加率（熱発生率傾き）と、その熱発生率最大時期における気筒内の燃料密度との相関が高く、機関運転条件の変化などによる影響を受け難い（つまり、熱発生率傾きの推定のロバスト性が高い）ことが解った。【００１６】
そこで、前記所定時期を熱発生率最大時期とし、混合気の着火時期から熱発生率最大時期までの熱発生率の平均的な増加率（熱発生率傾き）を、熱発生率最大時期における燃料密度に基づいて算出するのが好ましい。【００１７】
具体的に熱発生率傾きは、熱発生率最大時期の燃料密度の一次関数として算出することができる。また、熱発生率最大時期の燃料密度は、熱発生率が最大となるときに気筒内に存在する燃料の量を、そのときの気筒内の容積（筒内体積）によって除算することによって算出できる。好ましいのは、未燃分を除いて気筒内で実際に燃焼する燃料の量を筒内体積によって除算することであり、そのためには、前記のように熱発生率傾きを推定する際に、燃焼効率に基づく補正係数を乗算して算出すればよい。【００１８】
また、熱発生率傾きを推定する際に、機関回転速度に基づく補正係数（例えば機関回転速度の指数関数）を乗算するようにしてもよい。すなわち、機関回転速度が変化すると一般的に気筒内の流動の強さが変化することから、乱れの影響によって燃焼速度が変化することになる。そこで、機関回転速度に基づく補正を行うことにより、燃焼速度の変化を反映させて、熱発生率傾きをより精度良く推定することができる。【００１９】
このように気筒内の乱れの影響を考慮するのであれば、前記のように熱発生率傾きを推定する際に、熱発生率最大時期の筒内体積に基づく補正係数（例えば熱発生率最大時期の筒内体積の指数関数）を乗算するようにしてもよい。気筒の燃焼（膨張）行程においてはピストンの下降によって筒内体積が増大するに連れて乱れが弱くなり、燃焼速度が低下するからである。【００２０】
さらに、そうして熱発生率傾きを推定する際に、空燃比やＥＧＲ率に基づく補正係数（例えば空燃比の指数関数など）を乗算するようにしてもよく、こうすれば、空燃比やＥＧＲ率の変化による燃焼速度への影響を好適に反映させて、熱発生率傾きをより精度良く推定することができる。具体的には、例えば空燃比が理論空燃比よりもリーンになったり、ＥＧＲ率が高くなったりして、混合気の燃焼速度が低下することを熱発生率傾きに好適に反映させることができる。【００２１】
そして、前記のようにして算出した熱発生率傾きを用いて算出する熱発生率波形としては、例えば、混合気の着火から燃焼終了までのクランク角度期間を底辺とし、熱発生率最大時期における熱発生率を頂点とする三角波形が挙げられる。この三角波形によって熱発生率波形を近似すれば、着火時期から熱発生率最大時期までの熱発生率を表す三角形の斜辺の傾きが熱発生率傾きとなる。【００２２】
その場合に好ましいのは、混合気の着火時期から熱発生率最大時期までの期間が、機関負荷率、空燃比、ＥＧＲ率および油水温のいずれにも依らず、主に機関回転速度および点火時期に依って決まるものとして、前記三角波形を作成することである。すなわち、機関負荷率や空燃比、ＥＧＲ率などの変化によって混合気の燃焼速度が変化しても、着火時期から熱発生率最大時期までの期間は変化せず、熱発生率傾きの変化の分は、最大熱発生率が変化するものとして、前記三角波形を作成する。こうすることで、熱発生率波形を作成する工数のさらなる低減が図られる。【００２３】
見方を変えると本発明は、火花点火式の内燃機関における熱発生率波形の算出方法に係る。すなわち、まず、混合気の着火後に熱発生率が増加する熱発生率増加期間における、この熱発生率のクランク角度変化に対する増加率を、前記熱発生率波形の特性値の一つである熱発生率傾きとして規定する。そして、この熱発生率傾きを、前記熱発生率増加期間において予め設定された所定時期の気筒内の燃料密度に基づいて推定し、この推定した熱発生率傾きを用いて熱発生率波形を算出するものである。【００２４】
この場合も、前記所定の時期は熱発生率最大時期とし、混合気の着火時期から熱発生率最大時期までの期間における熱発生率の平均的な増加率（熱発生率傾き）を、熱発生率最大時期の燃料密度の一次関数として算出すればよい。
【発明の効果】
【００２５】
本発明によると、熱発生率波形の特性値の一つとして、混合気の着火後に熱発生率が増加する間の熱発生率傾きを用い、この熱発生率傾きを燃料密度に基づいて推定するようにしたので、熱発生率波形の作成（算出）にかかる工数を低減できると共に、熱発生率波形の全体を作成することなく、初期の燃焼速度について所要の精度を確保しながら、従来よりも簡易に燃焼状態を推定し、評価できるようになる。
【図面の簡単な説明】
【００２６】
【図１】実施形態に係る熱発生率波形算出装置の構成、および、この熱発生率波形算出装置の入出力情報を表す図である。
【図２】熱発生率波形算出装置から出力される熱発生率波形の一例を示す図である。
【図３】熱発生率波形算出装置において行われる熱発生率波形の作成手順を示すフローチャート図である。
【図４】ＢＴＤＣ着火の場合における、点火時期ＳＡでの筒内の燃料密度ρfuel@SAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示す図である。
【図５】式（１）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示す図である。
【図６】ＡＴＤＣ着火の場合における、着火時期ＦＡでの筒内の燃料密度ρfuel@FAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示す図である。
【図７】式（２）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示す図である。
【図８】ＢＴＤＣ着火する場合の点火時期ＳＡと熱発生率波形とを示す図である。
【図９】ＡＴＤＣ着火する場合の点火時期ＳＡと熱発生率波形とを示し、図９（ａ）は点火時期ＳＡがＢＴＤＣである場合を、図９（ｂ）は点火時期ＳＡがＡＴＤＣである場合をそれぞれ示す図である。
【図１０】負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。
【図１１】ＥＧＲ率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。
【図１２】空燃比のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。
【図１３】油水温のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。
【図１４】点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示した図である。
【図１５】エンジン回転速度Ｎｅのみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。
【図１６】あるエンジンに対して、式（３）で算出された予測前半燃焼期間と、実機において計測された実測前半燃焼期間との関係を検証した結果を示す図である。
【図１７】他のエンジンに対して、式（３）で算出された予測前半燃焼期間と、実機において計測された実測前半燃焼期間との関係を検証した結果を示す図である。
【図１８】負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。
【図１９】点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示した図である。
【図２０】熱発生率最大時燃料密度ρfuel@dQpeakと熱発生率傾きｂ／ａとの関係を、互いに異なるエンジン回転速度Ｎｅそれぞれに対して調べた実験の結果を示す図である。
【図２１】気筒内の燃焼場の状態を模式的に示す図であって、図２１（ａ）は気筒の軸線方向に見た平面図であり、図２１（ｂ）は側面図である。
【図２２】空燃比の変化による熱発生率傾きｂ／ａへの影響を示す図である。
【図２３】空燃比に基づく補正係数を乗算した熱発生率最大時燃料密度ρfuel@dQpeakと熱発生率傾きｂ／ａとの関係を検証した結果を示す図である。
【図２４】ＥＧＲ率の変化による熱発生率傾きへの影響を示す図である。
【図２５】ＥＧＲ率に基づく補正係数を乗算した熱発生率最大時燃料密度ρfuel@dQpeakと熱発生率傾きｂ／ａとの関係を検証した結果を示す図である。
【図２６】油水温の変化による熱発生率傾きへの影響を示す図である。
【図２７】油水温の影響を含む総熱発生量Ｑallを熱発生率最大時筒内体積Ｖ@dQpeakで除算して、熱発生率傾きｂ／ａとの関係を調べた実験の結果を示す図である。
【発明を実施するための形態】
【００２７】
以下、本発明の実施の形態を図面に基づいて説明する。本実施形態では、自動車用のガソリンエンジン（火花点火機関）を対象とした熱発生率波形を算出（作成）する熱発生率波形算出装置に本発明を適用した場合について説明する。
【００２８】
図１は、本実施形態に係る熱発生率波形算出装置１の構成、および、この熱発生率波形算出装置１の入出力情報を表す図である。この熱発生率波形装置１は、エンジンの状態量、制御パラメータの制御量および物理量の各種情報が入力される。これら入力情報としては、エンジン回転速度、負荷率、点火時期、ＥＧＲ率、空燃比、油水温、吸排気の各バルブの開閉タイミング（バルブタイミング）等が挙げられる。また、熱発生率波形算出装置１は、各入力情報に基づいて以下の各推定モデルを格納した推定部２〜５により熱発生率波形の各種特性値を推定して、この各種特性値を利用して作成した熱発生率波形を出力する。
【００２９】
−熱発生率波形の各特性値の推定部−
熱発生率波形算出装置１は、熱発生率波形の特性値として着火遅れ、前半燃焼期間、熱発生率傾きおよび熱発生量をそれぞれ推定するために、着火遅れ推定モデルを格納した着火遅れ推定部２、前半燃焼期間推定モデルを格納した前半燃焼期間推定部３、熱発生率傾き推定モデルを格納した熱発生率傾き推定部４、および、熱発生量推定モデルを格納した熱発生量推定部５を備えている。
【００３０】
着火遅れ推定部２は、エンジンの点火プラグによって混合気に点火された時期（以下、点火時期という）、即ち点火プラグの電極間に火花放電が行われてから、この火花によって混合気に着火し初期火炎核が形成される時期（以下、着火時期という）までの期間（以下、着火遅れ期間という）を、着火遅れ推定モデルを使用して推定する部分である。この着火遅れ期間はクランク角度［ＣＡ］で表される。なお、前記着火時期の定義として、本実施形態では、点火時期以降に熱発生率（クランクシャフトの回転の単位クランク角度当たりの熱発生量）が１［Ｊ／ＣＡ］に達した時期とする。この値はこれに限らず適宜設定可能であり、例えば、点火時期以降の熱発生量が総熱発生量に対して所定割合（例えば５％）に達した時期を着火時期としてもよい。さらに、総熱発生量に対する熱発生量の割合が所定値に達した時期（例えば１０％に達した時期でのクランク角度位置）と前記熱発生量の割合が他の所定値に達した時期（例えば５０％に達した時期でのクランク角度位置）とに基づいて着火時期を定義するようにしてもよい。つまり、これらのクランク角度位置と熱発生量の割合とによって、熱発生率が上昇していく期間での熱発生率波形に近似した三角形（三角波形）を作成し、この三角波形に基づいて着火時期を定義するものである。また、熱発生率が上昇していく期間での一般的な熱発生率波形の形状を、前記クランク角度位置と熱発生量の割合との関係が成立するように当て嵌めて熱発生率波形を作成し、この熱発生率波形に基づいて着火時期を定義するようにしてもよい。前記各値はこれに限定されることなく適宜設定可能である。
【００３１】
前半燃焼期間推定部３は、混合気の燃焼期間のうち、前記着火時期から、火炎核の成長に伴って熱発生率が最大となるタイミング（点火時期から燃焼終了時期までの期間中で熱発生率が最大となるタイミング）までの期間である前半燃焼期間を、前半燃焼期間推定モデルを使用して推定する部分である。以下、この熱発生率が最大となるタイミングを熱発生率最大時期という。また、この熱発生率最大時期および前半燃焼期間はそれぞれクランク角度［ＣＡ］で表される。
【００３２】
熱発生率傾き推定部４は、前記前半燃焼期間、即ち着火時期から熱発生率最大時期までの期間における、クランク角度変化に対する平均的な熱発生率の増加率（熱発生率の傾き）を、熱発生率傾き推定モデルを使用して推定する部分である。すなわち、本実施形態では、図２を参照して以下に述べるように、熱発生率波形に近似した三角波形を作成するようにしており、熱発生率傾き推定部４は、その三角波形における着火時期から熱発生率最大時期までの熱発生率を表す斜辺の傾きを推定するものである。この熱発生率の傾きの単位としては［Ｊ／ＣＡ2］で表される。
【００３３】
熱発生量推定部５は、混合気の燃焼によって発生した熱発生量（燃焼期間の全期間において発生した熱発生量であって、点火時期から燃焼終了時期までの期間中における熱発生率の積算値）を、熱発生量推定モデルを使用して推定する部分である。この熱発生量の単位としては［Ｊ］で表される。
【００３４】
各推定部２〜５での推定動作によってそれぞれ、着火遅れ、前半燃焼期間、熱発生率傾き、および熱発生量という熱発生率波形の特性値が求められ、これらの特性値を利用して熱発生率波形が作成される。そして、この作成された熱発生率波形が熱発生率波形算出装置１の出力となる。
【００３５】
このため、本実施形態に係る熱発生率波形算出装置１にあっては、図３に示すフローチャートの如く、着火遅れ推定部２における着火遅れ期間の推定動作（ステップＳＴ１）、前半燃焼期間推定部３における前半燃焼期間の推定動作（ステップＳＴ２）、熱発生率傾き推定部４における熱発生率傾きの推定動作（ステップＳＴ３）、熱発生量推定部５における熱発生量の推定動作（ステップＳＴ４）が順に行われた後、これら推定された特性値を利用した熱発生率波形の作成動作（ステップＳＴ５）が行われることになる。
【００３６】
図２は、各推定部２〜５で推定された特性値を利用して作成されて、熱発生率波形算出装置１から出力される熱発生率波形の一例を示している。この図２では、図中のタイミングＳＡが点火時期であり、図中のタイミングＦＡが着火時期である。このため、図中のτが着火遅れ期間となる。また、図中においてｄＱｐｅａｋＡが熱発生率最大時期であって、この熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率が図中のｂである。つまり、この熱発生率ｂが燃焼期間中における最大熱発生率となっている。また、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの期間である図中のａが前半燃焼期間となっている。このため、前半燃焼期間ａにおける熱発生率の傾きはｂ／ａとして表される。さらに、熱発生率最大時期ｄＱｐｅａｋＡから燃焼終了時期ＥＡまでの期間である図中のｃが後半燃焼期間となっている。また、図中のＱ１は前半燃焼期間ａにおける熱発生量であり、Ｑ２は後半燃焼期間ｃにおける熱発生量である。そして、燃焼期間の全期間において発生する熱発生量（総熱発生量Ｑall）は、これら熱発生量Ｑ１と熱発生量Ｑ２との和として表される。
【００３７】
言い換えると、本実施形態の熱発生率波形算出装置１は、混合気の着火から燃焼終了までのクランク角度期間（図中のＦＡからＥＡ）を底辺とし、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂを頂点とする三角波形によって熱発生率波形を近似するようにしている。本実施形態では、この熱発生率波形算出装置１の出力である熱発生率波形を利用し、エンジン設計に際してのシステムの検討、制御の検討、適合値の検討が行われる。
【００３８】
以下、各推定部２〜５での推定処理について具体的に説明する。
【００３９】
−着火遅れ推定部−
着火遅れ推定部２は、前述した如く、点火時期ＳＡから着火時期ＦＡまでの期間である着火遅れ期間τを推定する部分である。
【００４０】
この着火遅れ推定部２において行われる着火遅れ期間τの推定処理は以下のとおりである。
【００４１】
この着火遅れ期間τは、以下の式（１）および式（２）の何れか一方を利用して推定される（これらの式が着火遅れ推定モデルに相当する）。
【００４２】
【数１】
2015162983.tif 000003
【００４３】
【数２】
2015162983.tif 000004
【００４４】
ρfuel@SAは点火時期ＳＡにおける筒内の燃料密度（筒内燃料量［ｍｏｌ］／点火時期の筒内体積［Ｌ］）である。ρfuel@FAは着火時期ＦＡにおける筒内の燃料密度（筒内燃料量［ｍｏｌ］／着火時期の筒内体積［Ｌ］）である。Ｎｅはエンジン回転速度である。Ｃ1，Ｃ2，χ，δ，φ，ψはそれぞれ実験等に基づいて同定される係数である。
【００４５】
これら式（１）および式（２）は、空燃比が理論空燃比であり、ＥＧＲ率が「０」であり、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸排気の各バルブの開閉タイミングが固定されていることを条件に成立する式となっている。
【００４６】
式（１）は、ピストンが圧縮上死点に達するタイミング（ＴＤＣ）よりも進角側（ＢＴＤＣ）で混合気が着火した（以下、ＢＴＤＣ着火という）場合の着火遅れ期間τの算出式である。また、式（２）は、ピストンが圧縮上死点に達したタイミング（ＴＤＣ）よりも遅角側（ＡＴＤＣ）で混合気が着火した（以下、ＡＴＤＣ着火という）場合の着火遅れ期間τの算出式である。
【００４７】
これらの式に示すように、着火遅れ期間τは、所定タイミングでの筒内の燃料密度ρfuelおよびエンジン回転速度Ｎｅを変数とする演算式によって算出される。
【００４８】
これら演算式によって着火遅れ期間τが算出可能である根拠について以下に説明する。
【００４９】
図４は、ＢＴＤＣ着火の場合における、点火時期ＳＡでの筒内の燃料密度ρfuel@SAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示すグラフである。この実験は、空燃比を理論空燃比とし、ＥＧＲ率を「０」とし、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸排気の各バルブの開閉タイミングを固定して行われたものである。また、この図４では、「○」「△」「□」「◇」「×」「＋」「▽」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「△」は１０００ｒｐｍ、「□」は１２００ｒｐｍ、「◇」は１６００ｒｐｍ、「×」は２４００ｒｐｍ、「＋」は３２００ｒｐｍ、「▽」は３６００ｒｐｍである。
【００５０】
この図４に示すように、ＢＴＤＣ着火した場合、点火時期ＳＡにおける筒内の燃料密度ρfuel@SAと着火遅れ期間τとの間にはエンジン回転速度Ｎｅ毎に相関がある。つまり、これらの相関を概ね１本の曲線で表すことができる。図４では、エンジン回転速度Ｎｅが１０００ｒｐｍの場合および２４００ｒｐｍの場合のそれぞれについて、点火時期ＳＡにおける筒内の燃料密度ρfuel@SAと着火遅れ期間τとの相関を１本の曲線で表している。
【００５１】
図４に示すように、点火時期ＳＡにおける筒内の燃料密度ρfuel@SAが高いほど着火遅れ期間τは短くなっている。これは、燃料密度ρfuel@SAが高いほど、点火プラグ周辺の燃料分子の数が多くなり、点火プラグの点火後における火炎核の成長が急速に行われるためであると考えられる。また、エンジン回転速度Ｎｅは着火遅れ期間τに影響を与えている。つまり、エンジン回転速度Ｎｅが高いほど着火遅れ期間τは短くなっている。これは、エンジン回転速度Ｎｅが高いほど気筒内での混合気流れの乱れ（以下、単に乱れという）が強くなって前記火炎核の成長が急速に行われるためであると考えられる。このように、点火時期ＳＡにおける筒内の燃料密度ρfuel@SAおよびエンジン回転速度Ｎｅは着火遅れ期間τに影響を与えるパラメータとなっている。
【００５２】
図５は、式（１）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示すグラフである。この予測着火遅れ期間を求めるに当たっては、式（１）におけるＣ1，χ，δの各係数をエンジン運転条件に応じて同定することにより得られた予測式を使用している。この図５では、「○」「△」「□」「◇」「×」「＋」「▽」「☆」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「△」は１０００ｒｐｍ、「□」は１２００ｒｐｍ、「◇」は１６００ｒｐｍ、「×」は２０００ｒｐｍ、「＋」は２４００ｒｐｍ、「▽」は３２００ｒｐｍ、「☆」は３６００ｒｐｍである。
【００５３】
この図５から明らかなように、予測着火遅れ期間は実測着火遅れ期間に略一致しており、式（１）によって、ＢＴＤＣ着火した場合の着火遅れ期間が高い精度で算出されていることが解る。
【００５４】
図６は、ＡＴＤＣ着火の場合における、着火時期ＦＡでの筒内の燃料密度ρfuel@FAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示すグラフである。この実験は、エンジン回転速度を固定し、空燃比を理論空燃比とし、ＥＧＲ率を「０」とし、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸排気の各バルブの開閉タイミングを固定して行われたものである。また、この図６では、「○」「×」「＋」「△」の順で負荷率が高くなっている。例えば「○」は負荷率２０％、「×」は負荷率３０％、「＋」は負荷率４０％、「△」は負荷率５０％である。
【００５５】
この図６に示すように、ＡＴＤＣ着火した場合、着火時期ＦＡにおける筒内の燃料密度ρfuel@FAと着火遅れ期間τとの間には負荷率に依らず（負荷率に関わりなく）相関がある。つまり、これらの相関を概ね１本の曲線で表すことができる。
【００５６】
図６に示すように、着火時期ＦＡにおける筒内の燃料密度ρfuel@FAが高いほど着火遅れ期間τは短くなっている。これは、前述したように、燃料密度ρfuel@FAが高いほど、点火プラグ周辺の燃料分子の数が多くなり、点火プラグの点火後における火炎核の成長が急速に行われるためであると考えられる。このように、着火時期ＦＡにおける筒内の燃料密度ρfuel@FAは着火遅れ期間τに影響を与えるパラメータとなっている。また、前述した場合と同様に、エンジン回転速度Ｎｅも着火遅れ期間τに影響を与えるパラメータになっていると想定される。
【００５７】
図７は、式（２）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示すグラフである。この予測着火遅れ期間を求めるに当たっては、式（２）におけるＣ2，φ，ψの各係数をエンジン運転条件に応じて同定することにより得られた予測式を使用している。この図７では、「○」「×」「＋」「△」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「×」は１２００ｒｐｍ、「＋」は３６００ｒｐｍ、「△」は４８００ｒｐｍである。
【００５８】
この図７から明らかなように、予測着火遅れ期間は実測着火遅れ期間に略一致しており、式（２）によって、ＡＴＤＣ着火した場合の着火遅れ期間が高い精度で算出されていることが解る。
【００５９】
本発明の発明者は、これらの新たな知見に基づいて前記式（１）および式（２）を導き出した。
【００６０】
以下、着火遅れ期間τの算出に当たり、着火時期で層別する理由について説明する。つまり、ＢＴＤＣ着火した場合と、ＡＴＤＣ着火した場合とを層別し、それぞれに対して異なる演算式（前記式（１）および式（２））を用いて着火遅れ期間τを算出する理由について説明する。
【００６１】
まず、ＢＴＤＣ着火する場合、図８（点火時期ＳＡと熱発生率波形とを示す図）に示すように、点火時期ＳＡも、ピストンが圧縮上死点に達するタイミングよりも進角側（ＢＴＤＣ）である。この場合、点火時期ＳＡを迎えた後、ピストンは圧縮上死点に向かって移動することになる。つまり、筒内体積が小さくなっていき、これに伴って燃料密度ρfuelは大きくなっていく。このため、燃料密度ρfuelとしては、着火時期ＦＡでの燃料密度ρfuel@FAよりも点火時期ＳＡでの燃料密度ρfuel@SAの方が小さくなっている。そして、着火遅れ期間の最大値（想定される着火遅れ期間のうち最も長い着火遅れ期間）に相関のある値である点火時期ＳＡでの燃料密度ρfuel@SAに対して、予め同定された各種係数を乗算することにより、着火遅れ期間τを高い精度で求めることができる。
【００６２】
一方、ＡＴＤＣ着火する場合、図９（点火時期ＳＡと熱発生率波形とを示す図）に示すように、点火時期ＳＡは、ピストンが圧縮上死点に達するタイミングよりも進角側（ＢＴＤＣ）である場合（図９（ａ）を参照）と、遅角側（ＡＴＤＣ）である場合（図９（ｂ）を参照）とがある。これらの場合、着火時期ＦＡを迎えた後、ピストンは下死点に向かって移動することになる。つまり、筒内体積が大きくなっていき、これに伴って燃料密度ρfuelは小さくなっていく。このため、燃料密度ρfuelとしては、点火時期ＳＡでの燃料密度ρfuel@SAよりも着火時期ＦＡでの燃料密度ρfuel@FAの方が小さくなっている可能性が高い。そして、着火遅れ期間の最大値（想定される着火遅れ期間のうち最も長い着火遅れ期間）に相関のある値である着火時期ＦＡでの燃料密度ρfuel@FAに対して、予め同定された各種係数を乗算することにより、着火遅れ期間τを高い精度で求めることができる。
【００６３】
また、これら式（１）および式（２）の何れを使用するかを判断する手順（着火時期がＢＴＤＣ着火であるかＡＴＤＣ着火であるかを判断する手順）および着火遅れ期間（後述する真の着火遅れ期間）を算出する手順は以下のとおりである。まず、仮想の着火時期を設定し、この仮想の着火時期における筒内体積を求める。この筒内体積は仮想の着火時期に対応するクランク角度位置（ピストンの位置）によって幾何学的に求めることが可能であるので、仮想の着火時期から筒内体積は一義的に決まる。そして、この筒内体積と燃料噴射量とから燃料密度を求める。そして、仮想の着火時期をＢＴＤＣ着火として設定した場合には、この仮想の着火時期における燃料密度およびエンジン回転速度を式（１）に代入して推定着火遅れ期間を算出する。一方、仮想の着火時期をＡＴＤＣ着火として設定した場合には、この仮想の着火時期における燃料密度およびエンジン回転速度を式（２）に代入して推定着火遅れ期間を算出する。そして、仮想の着火時期に対し、前記算出された推定着火遅れ期間だけ進角した時期を仮想の点火時期とする。ここで、この仮想の点火時期と実際の点火時期（入力情報である点火時期）とを比較する。仮想の点火時期と実際の点火時期とが一致しない場合、前記仮想の着火時期を変更する。例えば仮想の着火時期を遅角側に変更する。そして、再度、この仮想の着火時期における燃料密度およびエンジン回転速度を式（１）または式（２）に代入して（仮想の着火時期をＢＴＤＣ着火として設定した場合には式（１）に代入し、仮想の着火時期をＡＴＤＣ着火として設定した場合には式（２）に代入する）推定着火遅れ期間を算出して、仮想の点火時期を求め、これと実際の点火時期（入力情報である点火時期）とを比較する。この動作を繰り返し、仮想の点火時期と実際の点火時期とが一致した場合における仮想の着火時期が真の着火時期として得られることになる。そして、この際（真の着火時期が得られた際）、式（１）または式（２）において算出された推定着火遅れ期間が真の着火遅れ期間として得られることになる。また、前記真の着火時期がＢＴＤＣであった場合（ＢＴＤＣ着火であった場合）に前記求められた着火時期を改めて式（１）に代入して着火遅れ期間τの算出を行い、前記真の着火時期がＡＴＤＣであった場合（ＡＴＤＣ着火であった場合）に前記求められた着火時期を改めて式（２）に代入して着火遅れ期間τの算出を行うようにしてもよい。
【００６４】
以上の動作を言い換えると以下のとおりとなる。前記実際の点火時期と仮想の着火時期との間の期間（仮想の着火時期で着火したとした場合の仮想の着火遅れ期間）と、式（１）または式（２）で算出された（推定された）推定着火遅れ期間とを比較し、これらが一致しない場合には、前記仮想の着火時期を変更する。そして、再度、式（１）または式（２）によって推定着火遅れ期間を算出した後、前記実際の点火時期と仮想の着火時期との間の期間（仮想の着火遅れ期間）と、式（１）または式（２）で算出された推定着火遅れ期間とを比較する。この動作を繰り返し、これらが一致した場合（仮想の着火遅れ期間と推定着火遅れ期間とが一致した場合）における推定着火遅れ期間が真の着火遅れ期間として得られることになる。
【００６５】
以上のような着火遅れ推定部２による着火遅れ期間τの推定により、エンジン運転領域の全領域に対して着火遅れ期間τの推定が可能になる。
【００６６】
以上のようにして着火遅れ期間τが算出されると、前記点火時期ＳＡに着火遅れ期間τを加算することで着火時期ＦＡを求めることができる。
【００６７】
−前半燃焼期間推定部−
前半燃焼期間推定部３は、前述した如く、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの期間である前半燃焼期間ａを推定する部分である。
【００６８】
この前半燃焼期間推定部３において行われる前半燃焼期間ａの推定処理は以下のとおりである。
【００６９】
この前半燃焼期間ａ［ＣＡ］は、以下の式（３）を利用して推定される（この式が前半燃焼期間推定モデルに相当する）。
【００７０】
【数３】
2015162983.tif 000005
【００７１】
Ｖ@dQpeakは前記熱発生率最大時期ｄＱｐｅａｋＡにおける筒内体積［Ｌ］であり、以下では、熱発生率最大時筒内体積ともいう。Ｎｅはエンジン回転速度である。Ｃ，α，βはそれぞれ実験等に基づいて同定される係数である。
【００７２】
この式（３）は、吸気バルブの開閉タイミングが固定されていることを条件に成立する式となっている。また、この式（３）は、負荷率、ＥＧＲ率、空燃比、油水温の影響を受けることなく成立するものとなっている。つまり、式（３）は、前半燃焼期間ａが、負荷率、ＥＧＲ率、空燃比、油水温の影響を受けないことに基づいて成立している。
【００７３】
この式（３）によって前半燃焼期間ａが算出可能である根拠について以下に説明する。
【００７４】
図１０〜図１３それぞれは、互いに異なるエンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示したものである。図１０は、負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。図１１は、ＥＧＲ率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。図１２は、空燃比のみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。また、図１３は、エンジンの暖機運転の途中のように油水温のみが互いに異なる場合に得られる熱発生率波形を重ねて表示したものである。
【００７５】
これら図１０〜図１３に示すように、負荷率、ＥＧＲ率、空燃比、油水温の何れが変化しても前半燃焼期間ａは一定に維持されている。つまり、前半燃焼期間ａは、負荷率、ＥＧＲ率、空燃比、油水温の影響を受けないものであることが解る。
【００７６】
一方、図１４は、点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。この図１４から解るように、点火時期ＳＡが遅角されるほど前半燃焼期間ａは長くなっている。
【００７７】
また、図１５は、エンジン回転速度Ｎｅのみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示したものである。エンジン回転速度Ｎｅが高いほど単位時間［ｍｓ］当たりにおけるクランクの回転角度［ＣＡ］は大きくなるので、その分、前半燃焼期間ａは長くなる（クランク角度軸上で長くなる）はずであるが、図１５に示すものでは、エンジン回転速度Ｎｅが異なっても前半燃焼期間ａは殆ど変化していない。これは、エンジン回転速度Ｎｅが高いほど前半燃焼期間ａが短くなる要因が存在しているものと考えられる。つまり、エンジン回転速度Ｎｅが高いほど単位時間当たりにおけるクランクの回転角度が大きくなることに起因して前半燃焼期間ａが長くなることとは別に、「他の要因」によって前半燃焼期間ａが短くなっていると想定できる。
【００７８】
このように前半燃焼期間ａは、点火時期ＳＡおよびエンジン回転速度Ｎｅの影響を受けるものであることが解る。
【００７９】
前半燃焼期間ａが、点火時期ＳＡおよびエンジン回転速度Ｎｅの影響を受ける要因としては、点火時期ＳＡおよびエンジン回転速度Ｎｅが、気筒内の乱れに影響を与えることが考えられる。
【００８０】
つまり、点火時期ＳＡが遅角側に移行するほど着火時期ＦＡおよび熱発生率最大時期ｄＱｐｅａｋＡも遅角側に移行し、この熱発生率最大時期ｄＱｐｅａｋＡにおける筒内体積（熱発生率最大時筒内体積Ｖ@dQpeak）が大きくなると共に、気筒内の乱れは弱くなる。そして、気筒内の乱れが弱くなると、火炎伝播が緩慢になって前半燃焼期間ａは長くなる。逆に、点火時期ＳＡが進角側に移行するほど着火時期ＦＡおよび熱発生率最大時期ｄＱｐｅａｋＡも進角側に移行し、熱発生率最大時筒内体積Ｖ@dQpeakが小さくなると共に、気筒内の乱れは強くなる。これにより火炎伝播が急速になって、前半燃焼期間ａは短くなる。
【００８１】
また、エンジン回転速度Ｎｅが低いほど吸気系から気筒内に流入する空気の流速が低くなって気筒内の乱れは弱くなる。そして、気筒内の乱れが弱くなると、火炎伝播が緩慢になって前半燃焼期間ａは長くなる。逆に、エンジン回転速度Ｎｅが高いほど吸気系から気筒内に流入する空気の流速が高くなって気筒内の乱れは強くなる。そして、気筒内の乱れが強くなると、火炎伝播が急速になって前半燃焼期間ａは短くなる。前述した「他の要因（前半燃焼期間ａを短くする要因）」は、このエンジン回転速度Ｎｅが高いほど気筒内の乱れが強くなることに起因して火炎伝播が急速になっていることである。
【００８２】
本発明の発明者は、この新たな知見に基づいて前記式（３）を導き出した。そして、この式（３）においては、制御量である点火時期ＳＡに相関のある物理量として筒内体積、特に、熱発生率最大時筒内体積Ｖ@dQpeakを変数として用いている。つまり、前述した如く、点火時期ＳＡが遅角側に移行するほど熱発生率最大時期ｄＱｐｅａｋＡも遅角側に移行し、筒内体積Ｖ@dQpeakが大きくなることから、点火時期ＳＡに相関のある物理量として、熱発生率最大時筒内体積Ｖ@dQpeakを変数として用いている。
【００８３】
前記式（３）の変数である熱発生率最大時筒内体積Ｖ@dQpeakを求める手順および前半燃焼期間ａを算出する手順は以下のとおりである。まず、仮想の熱発生率最大時期を設定し、この仮想の熱発生率最大時期における筒内体積を求める。この筒内体積は仮想の熱発生率最大時期に対応するクランク角度位置（ピストンの位置）によって幾何学的に求めることが可能であるので、仮想の熱発生率最大時期から筒内体積は一義的に決まる。そして、この仮想の熱発生率最大時期における筒内体積およびエンジン回転速度を式（３）に代入して推定前半燃焼期間を算出する。そして、前記仮想の熱発生率最大時期に対し、この算出した推定前半燃焼期間だけ進角した時期を仮想の着火時期とする。前述した着火遅れ推定部２では着火遅れ期間τが算出されているため、前記点火時期ＳＡに着火遅れ期間τを加算することで着火時期ＦＡが算出できる。ここで、前記仮想の着火時期と算出された着火時期ＦＡとを比較する。仮想の着火時期と算出された着火時期ＦＡとが一致しない場合、前記仮想の熱発生率最大時期を変更する。例えば仮想の熱発生率最大時期を遅角側に変更する。そして、再度、この仮想の熱発生率最大時期における筒内体積およびエンジン回転速度を式（３）に代入して推定前半燃焼期間を算出して、仮想の着火時期を求め、これと算出された着火時期ＦＡ（前記点火時期ＳＡに、着火遅れ推定部２で算出された着火遅れ期間τを加算することで得られる着火時期ＦＡ）とを比較する。この動作を繰り返し、仮想の着火時期と算出された着火時期ＦＡとが一致した場合における仮想の熱発生率最大時期が真の熱発生率最大時期ｄＱｐｅａｋＡとして得られることになる。そして、この際（真の熱発生率最大時期ｄＱｐｅａｋＡが得られた際）、式（３）において算出された推定前半燃焼期間が真の前半燃焼期間として得られることになる。また、前記真の熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率最大時筒内体積Ｖ@dQpeakを幾何学的に求め、これを改めて式（３）に代入して前半燃焼期間ａの算出を行うようにしてもよい。
【００８４】
以上の動作を言い換えると以下のとおりとなる。前記着火時期ＦＡ（実際の点火時期に従って求められた着火時期）と仮想の熱発生率最大時期との間の期間（仮想の前半燃焼期間）と、式（３）で算出された（推定された）推定前半燃焼期間（仮想の熱発生率最大時期における物理量に基づいて推定された前半燃焼期間）とを比較し、これらが一致しない場合には、前記仮想の熱発生率最大時期を変更する。そして、再度、式（３）によって推定前半燃焼期間を算出した後、前記着火時期ＦＡと仮想の熱発生率最大時期との間の期間（仮想の前半燃焼期間）と、式（３）で算出された推定前半燃焼期間とを比較する。この動作を繰り返し、これらが一致した場合（仮想の前半燃焼期間と推定前半燃焼期間とが一致した場合）における推定前半燃焼期間が真の前半燃焼期間ａとして得られることになる。
【００８５】
式（３）における各係数について具体的に説明すると、Ｃおよびαは実験等に基づいた同定によって得られる。また、βは気筒内でのタンブル比に応じた値となっており、タンブル比が大きいほど大きな値として与えられる。なお、各係数それぞれを実験等に基づいた同定によって設定するようにしてもよい。また、これら係数は、吸気バルブの開閉タイミングの変化に対して同定することも可能である。
【００８６】
図１６および図１７は、互いに異なるエンジンに対して、式（３）で算出された予測前半燃焼期間と、実機において計測された実測前半燃焼期間との関係を検証した結果を示すグラフである。この予測前半燃焼期間を求めるに当たっては、式（３）における係数Ｃをエンジン運転条件に応じて同定することにより得られた予測式を使用している。図１６では、「○」「△」「□」「◇」「×」「＋」「▽」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「△」は１０００ｒｐｍ、「□」は１２００ｒｐｍ、「◇」は１６００ｒｐｍ、「×」は２４００ｒｐｍ、「＋」は３２００ｒｐｍ、「▽」は３６００ｒｐｍである。また、図１７では、「○」「×」「＋」「△」「□」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「×」は１２００ｒｐｍ、「＋」は２４００ｒｐｍ、「△」は３６００ｒｐｍ、「□」は４８００ｒｐｍである。
【００８７】
これら図１６および図１７から明らかなように、予測前半燃焼期間は実測前半燃焼期間に略一致しており、式（３）によって前半燃焼期間ａが高い精度で算出されていることが解る。
【００８８】
以上、述べたように前半燃焼期間ａを、負荷率、空燃比、ＥＧＲ率、油水温の影響を受けないものとして、熱発生率最大時筒内体積Ｖ@dQpeakとエンジン回転速度Ｎｅとに基づいて推定できる。これら熱発生率最大時筒内体積Ｖ@dQpeakとエンジン回転速度Ｎｅは前述した如く気筒内の乱れに相関のあるパラメータとなっている。言い換えると、負荷率、空燃比、ＥＧＲ率、油水温は、気筒内の乱れに殆ど相関が無いことから前半燃焼期間ａに影響を与えないものと想定される。そして、これら負荷率、空燃比、ＥＧＲ率、油水温を考慮することなしに、気筒内の乱れに相関のあるパラメータである熱発生率最大時筒内体積Ｖ@dQpeakとエンジン回転速度Ｎｅとに基づいて前半燃焼期間ａを推定できることから、エンジンの種々の運転条件における前半燃焼期間ａを決定するための工数が大幅に低減できることになる。
【００８９】
なお、前述したように、前半燃焼期間は負荷率の影響を受けない。この負荷率は、燃料噴射量を制御するためのパラメータの一つであり、燃料噴射量は気筒内の燃料密度に影響を及ぼす制御パラメータである。このため、前半燃焼期間は、気筒内の燃料密度に依ることなく推定される。具体的には、前述したように、熱発生率最大時筒内体積Ｖ@dQpeakおよびエンジン回転速度Ｎｅといった気筒内の乱れに影響を与えるパラメータに基づいて前半燃焼期間は推定される。これに対し、後述するように、熱発生率傾きは気筒内の燃料密度に基づいて推定される。このように、本実施形態で推定される前半燃焼期間と熱発生率傾きとは、互いに独立した（従属関係にない）値として推定されるものとなっている。
【００９０】
−熱発生率傾き推定部−
熱発生率傾き推定部４は、前述した如く、前半燃焼期間ａにおける熱発生率の傾きｂ／ａ（以下、熱発生率傾きという）を推定する部分である。
【００９１】
この熱発生率傾き推定部４において行われる熱発生率傾きｂ／ａの推定処理は以下のとおりである。
【００９２】
この熱発生率傾きｂ／ａ［Ｊ／ＣＡ2］は、基本的には以下の式（４）を利用して推定される（この式が熱発生率傾き推定モデルに相当する）。
【００９３】
【数４】
2015162983.tif 000006
【００９４】
ρfuel@dQpeakは、前記熱発生率最大時期ｄＱｐｅａｋＡにおける燃料密度（筒内燃料量［ｍｏｌ］／熱発生率最大時期の筒内体積［Ｌ］）であり、以下では、熱発生率最大時燃料密度ともいう。Ｃ3は実験等に基づいて同定される係数である。
【００９５】
この式（４）は、エンジン回転速度が固定されており、空燃比が理論空燃比であり、ＥＧＲ率が「０」であり、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸気バルブの開閉タイミングが固定されていることを条件に成立する式となっている。なお、エンジン回転速度、空燃比、ＥＧＲ率、エンジンの油水温等の影響については後述する。
【００９６】
この式（４）によって熱発生率傾きｂ／ａが算出可能である根拠について以下に説明する。
【００９７】
図１８（ａ）〜図１８（ｄ）それぞれは、負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡが互いに一致するように点火時期ＳＡを調整した各熱発生率波形を重ねて表示したものである。図１８（ａ）〜図１８（ｄ）の順に点火時期は遅角側に変化しており、また、各図において負荷率は、ＫＬ１，ＫＬ２，ＫＬ３の順に大きくなっている。例えば図１８において、ＫＬ１は負荷率２０％、ＫＬ２は負荷率３０％，ＫＬ３は負荷率４０％である。
【００９８】
これら図１８（ａ）〜図１８（ｄ）に表れているように、熱発生率傾きｂ／ａは負荷率および点火時期ＳＡの影響を受けている。具体的には、それぞれ点火時期ＳＡの異なる図１８（ａ）〜図１８（ｄ）のいずれにおいても、負荷率が大きいほど熱発生率傾きｂ／ａは大きくなっている。このように熱発生率傾きｂ／ａが負荷率の影響を受ける要因として、負荷率に応じて気筒内の燃料密度が変化することが考えられる。つまり、負荷率が高いほど筒内燃料量が多くなるので、気筒内の燃料密度も高くなって、混合気の燃焼速度が高くなるものと考えられる。
【００９９】
また、図１８（ａ）〜図１８（ｄ）の順に点火時期ＳＡが遅角側に移行するに連れて、熱発生率傾きｂ／ａが小さくなっている。図１９（ａ）および図１９（ｂ）それぞれは、点火時期ＳＡの変化による影響を調べるために、点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。これら図１９（ａ）および図１９（ｂ）においてはそれぞれ負荷率が異なっているが、そのいずれにおいても熱発生率傾きｂ／ａは、点火時期ＳＡの遅角側への移行に連れて小さくなる傾向がある。
【０１００】
このように熱発生率傾きｂ／ａが点火時期ＳＡの影響を受ける要因も、前記した負荷率の場合と同じく気筒内の燃料密度によるものと考えられる。すなわち、ピストンが圧縮上死点（ＴＤＣ）付近にあるときは、クランク角度の変化に伴う筒内体積の変化は小さいが、膨張行程においてＴＤＣから離れてゆくに連れて（例えばＡＴＤＣ１０°ＣＡくらいからは）遅角側ほど筒内体積が大きくなってゆき、これに連れて気筒内の燃料密度は低下してゆく。
【０１０１】
そして、前記図１９（ａ）および図１９（ｂ）に表れているように熱発生率波形は、点火時期ＳＡの遅角に伴い全体として遅角側に移動するとともに、着火時期ＦＡ（波形の始点）がＴＤＣ以降になると、徐々に熱発生率波形の傾きも小さくなってゆく。この結果、着火時期ＦＡ（波形の始点）から熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂ（波形の頂点）までを結んだ直線（図中には一点鎖線で示す）の傾き、即ち熱発生率傾きｂ／ａも遅角側に向かって徐々に小さくなっている。
【０１０２】
そのように点火時期ＳＡの遅角（即ち着火時期ＦＡの遅角）が熱発生率傾きｂ／ａに及ぼす影響は、この熱発生率傾きｂ／ａと熱発生率最大時燃料密度ρfuel@dQpeakとの関係に顕著に表れている。すなわち、図１９（ａ）および図１９（ｂ）に表れているように、点火時期ＳＡの遅角に伴い熱発生率最大時期ｄＱｐｅａｋＡが遅角側に移行し、この熱発生率最大時期ｄＱｐｅａｋＡにおける筒内体積（熱発生率最大時筒内体積Ｖ@dQpeak）が徐々に大きくなってゆくと、これに応じて熱発生率最大時燃料密度ρfuel@dQpeakが小さくなってゆき、これに対応して熱発生率傾きｂ／ａが小さくなってゆく。
【０１０３】
本発明者は、そのように熱発生率最大時燃料密度ρfuel@dQpeakの変化に対応して熱発生率傾きｂ／ａが変化する様子を調べた。この実験の結果を図２０（ａ）〜図２０（ｄ）のグラフに示す。これらの各図において負荷率は、「○」「×」「＋」「△」「□」「◇」「▽」「☆」の順に大きくなっている。例えば図２０において、「○」は負荷率１５％、「×」は負荷率２０％、「＋」は負荷率２５％、「△」は負荷率３０％、「□」は負荷率３５％、「◇」は負荷率４０％、「▽」は負荷率４５％、「☆」は負荷率５０％である。
【０１０４】
また、図２０（ａ）〜図２０（ｄ）の順にエンジン回転速度Ｎｅが高くなっており、例えば、図２０（ａ）は８００ｒｐｍ、図２０（ｂ）は１２００ｒｐｍ、図２０（ｃ）は２０００ｒｐｍ、図２０（ｄ）は３２００ｒｐｍである。
【０１０５】
図２０（ａ）〜図２０（ｄ）それぞれに表れているように、エンジン回転速度を固定すれば、熱発生率最大時燃料密度ρfuel@dQpeakと熱発生率傾きｂ／ａとの関係は、負荷率や点火時期ＳＡが異なっていても概ね１本の直線として表すことができ、両者の間には高い相関（具体的には概略比例関係）のあることが解る。つまり、熱発生率傾きｂ／ａに対するエンジンの負荷率および点火時期ＳＡの影響は、熱発生率最大時燃料密度ρfuel@dQpeakという一つのパラメータによって、まとめて表現されるものである。
【０１０６】
本発明の発明者は、この新たな知見に基づいて前記式（４）を導き出した。
【０１０７】
前記式（４）の変数である熱発生率最大時燃料密度ρfuel@dQpeakは、前述したように筒内燃料量を熱発生率最大時筒内体積Ｖ@dQpeakによって除算して、求めることができる。熱発生率最大時筒内体積Ｖ@dQpeakを求める手順は、前半燃焼期間推定部３の説明において上述したとおりである。また、筒内燃料量は、熱発生率波形算出装置１の入力情報として与えられる。
【０１０８】
このようにして、熱発生率波形の特性値の一つである熱発生率傾きｂ／ａを、基本的には熱発生率最大時燃料密度ρfuel@dQpeakの一次関数（この実施形態では、一例として比例関数）として算出することができる。言い換えると、負荷率および点火時期ＳＡを考慮せず、主に熱発生率最大時燃料密度ρfuel@dQpeakに基づいて熱発生率傾きｂ／ａを推定できることから、エンジンの種々の運転条件における熱発生率傾きｂ／ａの決定にかかる工数が低減される。
【０１０９】
（燃焼効率およびエンジン回転速度に基づく補正）
ここで、前述した熱発生率傾きｂ／ａと熱発生率最大時燃料密度ρfuel@dQpeakとの相関において、この熱発生率最大時燃料密度ρfuel@dQpeakを算出する筒内燃料量としては、未燃分を除いて気筒内で実際に燃焼する量を用いることが好ましい。すなわち、エンジンの暖機運転中のように油水温が低い場合には燃焼効率が低下するし、それ以外にも負荷率や熱発生率最大時期ｄＱｐｅａｋＡが変化した場合は、燃焼効率が変化する可能性がある。そこで、燃焼効率に基づく補正係数ｋを前記式（４）の右辺に乗算することにより、以下の式（５）が得られる。
【０１１０】
【数５】
2015162983.tif 000007
【０１１１】
なお、前述の如く熱発生率最大時燃料密度ρfuel@dQpeakは、筒内燃料量を熱発生率最大時筒内体積Ｖ@dQpeakによって除算して得られるから、式（５）は以下の式（６）のように表すこともできる。
【０１１２】
【数６】
2015162983.tif 000008
【０１１３】
また、式（６）における「筒内燃料量×ｋ」は、燃焼期間の全期間において発生する熱発生量（総熱発生量Ｑall）として置き換えることもできるから、以下の式（７）のように表すこともできる。
【０１１４】
【数７】
2015162983.tif 000009
【０１１５】
このような燃焼効率の影響を示す一例として、上述した図１３のようにエンジンの油水温が異なる場合の熱発生率波形を見ると、図２６に示すように、熱発生率傾きｂ／ａが油水温の影響を受けていることが解る。すなわち、図中の油水温はＴ０，Ｔ１，Ｔ２の順に低くなっており、例えばＴ０は、エンジンの暖機運転の完了後に相当する８８℃である。また、例えばＴ１は４０℃、Ｔ２は２０℃である。
【０１１６】
図２６に表れているように油水温が低いほど、熱発生率傾きｂ／ａは小さくなっており、これは、油水温が低いほど燃料の未燃分や冷損が増大して、燃焼効率が低下することによると考えられる。この燃焼効率の影響が総熱発生量Ｑallに含まれていることに着目し、前記の式（７）のように総熱発生量Ｑallを熱発生率最大時筒内体積Ｖ@dQpeakで除算した値と、熱発生率傾きｂ／ａとの関係を調べた結果を図２７に示す。
【０１１７】
図２７においては油水温を、「○」「×」「＋」「△」「□」の順に高くしており、例えば「○」は２０℃、「×」は４０℃、「＋」は６０℃、「△」は８０℃、「□」は８８℃である。この図２７に表れているように、Ｑall／Ｖ@dQpeakと熱発生率傾きｂ／ａとの関係は、油水温が大きく異なっていても概ね１本の直線として表すことができ、両者の間には高い相関（この例では概略比例関係）のあることが解る。つまり、燃焼効率に基づく補正係数ｋを用いることによって、油水温の変化などによる影響も好適に反映させて、熱発生率波形ｂ／ａを算出することができる。
【０１１８】
一方、エンジン回転速度Ｎｅの影響については、前記の図２０（ａ）〜図２０（ｄ）に表れているように、熱発生率最大時燃料密度ρfuel@dQpeakと熱発生率傾きｂ／ａとの相関関係がエンジン回転速度Ｎｅによって変化している。具体的には図２０（ａ）〜図２０（ｄ）の順にエンジン回転速度Ｎｅが高くなるに連れて、熱発生率傾きｂ／ａは小さくなる傾向がある。そこで、このエンジン回転速度Ｎｅに基づく補正係数Ｎｅεを前記式（４）の右辺に乗算して、以下の式（８）を用いることもできる。
【０１１９】
【数８】
2015162983.tif 000010
【０１２０】
εは実験等に基づいて同定される値であるが、前記のようにエンジン回転速度Ｎｅが高くなるに連れて、熱発生率傾きｂ／ａは小さくなる傾向があるので、εは一般的には負の値となる。
【０１２１】
前記エンジン回転速度Ｎｅの影響については、前半燃焼期間推定部３の説明において上述したように、エンジン回転速度Ｎｅが高いほど時間当たりのクランク角度変化は大きくなるので、その分、熱発生率傾きｂ／ａは小さくなるはずであるが、前記図２０（ａ）および図２０（ｄ）を比較すると、エンジン回転速度Ｎｅが４倍になっても熱発生率傾きｂ／ａは１／４にはなっていない。
【０１２２】
これは、前半燃焼期間推定部３の説明において上述したように、エンジン回転速度Ｎｅが高いほど気筒内の乱れが強くなることに起因して、燃焼速度が高くなっていることによると考えられる。つまり、エンジン回転速度Ｎｅが高いほど燃焼速度が高くなるため、時間当たりのクランク角度変化が大きくなるほどには、熱発生率傾きｂ／ａが小さくならないと考えられる。
【０１２３】
（熱発生率最大時筒内体積に基づく補正）
ところで、前述のように熱発生率傾きｂ／ａが熱発生率最大時燃料密度ρfuel@dQpeakに基づいて算出できる理由は、以下のように考えられる。まず、気筒内の燃焼場の状態を鑑みると熱発生率傾きｂ／ａは、乱流燃焼速度や乱流火炎面面積の影響を受ける。図２１は、気筒内の燃焼場の状態を模式的に示した図であって、図２１（ａ）は気筒軸線に沿って見た平面図であり、図２１（ｂ）は側面図である（図中の符号６はピストンを示す）。また、この図２１において、実線の斜線を付した領域が既燃部であり、破線の斜線を付した領域が未燃部である。
【０１２４】
図２１には模式的に矢印で示すように、既燃部の外縁部分において未燃部に向かって進行する火炎面Ｆの面積が乱流火炎面面積であり、その進行する速度（気筒内において外周側に進む火炎面Ｆの速度）が乱流燃焼速度である。これら乱流火炎面面積と乱流燃焼速度とを乗算したものが単位時間当たりの燃焼場の体積であるから、これに燃料密度を乗算すれば熱発生率が得られる。よって、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率（最大熱発生率ｂ）は、以下の式（９）で与えられる。
【０１２５】
【数９】
2015162983.tif 000011
【０１２６】
Ａｔ@dQpeakは、熱発生率最大時期ｄＱｐｅａｋＡにおける乱流火炎面面積であり、Ｓｔ@dQpeakは、熱発生率最大時期ｄＱｐｅａｋＡにおける乱流燃焼速度である。前記の式（９）から熱発生率傾きｂ／ａは、以下の式（１０）として表される。
【０１２７】
【数１０】
2015162983.tif 000012
【０１２８】
ここで一般的に、乱流燃焼速度Ｓｔは気筒内の乱れ強さｕ’に概ね比例するとみなしてよい。前述したように乱れが強いときには前半燃焼期間ａが短くなることから、乱流燃焼速度Ｓｔは概略的には前半燃焼期間ａに反比例すると見なしてよい。一方、ＡＴＤＣ、即ち気筒の膨張行程における乱流火炎面面積Ａｔは、筒内体積Ｖとの間に相関があり、遅角側ほど気筒内において外周側に進んだ火炎の面積が大きくなるから、乱流火炎面面積Ａｔは概略的には筒内体積Ｖの指数関数に比例すると見なしてよい。よって、前記の式（１０）から以下の式（１１）が得られる。なお、η、Ｃ4は実験等に基づいて同定される。
【０１２９】
【数１１】
2015162983.tif 000013
【０１３０】
また、前半燃焼期間ａについては、上述した式（３）においてエンジン回転速度Ｎｅを固定すれば、  ａ  ＝  Ｃ×Ｖ@dQpeakα  となるから、この「ａ」を前記の式（１１）に代入して、以下の式（１２）を得る。
【０１３１】
【数１２】
2015162983.tif 000014
【０１３２】
この式（１２）において、ηの値は概ね「１」とみなしてよく、また、αの値は概ね「１／２」とみなしてよいので、Ｖ@dQpeakη-2αの値は概略的に「１」とみなすことができ、こうすれば式（１２）から前記の式（４）が得られるのである。
【０１３３】
但し、前記のＶ@dQpeakη-2α、即ち熱発生率最大時筒内体積の指数関数は、前述したエンジン回転速度Ｎｅと同じく気筒内の乱れの影響を表すものと考えられるので、より高い精度で熱発生率傾きｂ／ａを求めるのであれば、このＶ@dQpeakη-2αも残すことが好ましい。すなわち、指数「η−２α」を「κ」に置き換えて、熱発生率最大時期の筒内体積Ｖ@dQpeakに基づく補正係数として式（８）の右辺に乗算すれば、熱発生率傾きｂ／ａを算出する以下の式（１３）が得られる。なお、κは実験等に基づいて同定される。
【０１３４】
【数１３】
2015162983.tif 000015
【０１３５】
この式（１３）は、式（１０）における「Ａｔ@dQpeak×Ｓｔ@dQpeak／ａ」をＶ@dQpeakκとしてまとめたものと言うことができる。また、その式（１３）においても、上述したように燃焼効率の影響を反映させるべく右辺に補正係数ｋを乗算し、以下の式（１４）によって熱発生率傾きｂ／ａを算出するようにしてもよい。
【０１３６】
【数１４】
2015162983.tif 000016
【０１３７】
また、熱発生率最大時燃料密度ρfuel@dQpeakを筒内燃料量および熱発生率最大時筒内体積Ｖ@dQpeakによって表し、以下の式（１５）のように表すこともできる。この式（１５）における「筒内燃料量×ｋ」を総熱発生量Ｑallで置き換えて、以下の式（１６）のように表すこともできる。
【０１３８】
【数１５】
2015162983.tif 000017
【０１３９】
【数１６】
2015162983.tif 000018
【０１４０】
（空燃比およびＥＧＲ率に基づく補正）
次に、熱発生率傾きｂ／ａに対する空燃比およびＥＧＲ率の影響について説明する。図２２は、空燃比の変化による熱発生率波形の変化を示す実験結果の図であって、図２２（ａ）、図２２（ｂ）それぞれは、図１２と同じく、熱発生率最大時期ｄＱｐｅａｋＡが互いに一致するように点火時期ＳＡを調整しながら、空燃比の変化によって得られる熱発生率波形を重ねて表示している。
【０１４１】
詳しくは図２２（ａ）は、理論空燃比を超える（リーンとなる）ように空燃比を変化させた場合であり、図２２（ｂ）は反対に理論空燃比未満となる（リッチとなる）ように空燃比を変化させた場合である。例えば図２２（ａ）および図２２（ｂ）においてＡＦ０は理論空燃比（Ａ／Ｆ＝１４．６）であり、図２２（ａ）においてＡＦ１はＡ／Ｆ＝１９．５である。一方、図２２（ｂ）において、ＡＦ２はＡ／Ｆ＝１２．０である。
【０１４２】
これら図２２（ａ）、図２２（ｂ）に表れているように、熱発生率傾きｂ／ａは空燃比の影響を受けていて、空燃比が大きい（リーンな）ほど小さな値になっている。このように熱発生率傾きｂ／ａが空燃比の影響を受ける要因として、理論空燃比よりも大きな（リーンな）空燃比では、燃料と酸素との反応（燃焼）に寄与しない空気中の窒素および酸素の分子が火炎面の伝播を邪魔することが挙げられる。この場合、空燃比が大きいほど余剰な空気が多いことから、燃焼速度が低下する。
【０１４３】
一方で、空燃比が理論空燃比よりもやや小さな（リッチな）状態では、燃料の分子と酸素との邂逅率が高くなることによって、混合気の燃焼速度が高くなると考えられるが、さらに空燃比が小さく（リッチに）なると燃焼状態は急激に悪化し、やがては失火に至る。図２２（ｂ）の例では、Ａ／Ｆ＝１２．０よりも小さく（リッチに）なると、燃焼状態が悪化すると考えられる。
【０１４４】
本発明者は、そのように空燃比の変化に対応して熱発生率傾きｂ／ａが変化する様子を調べた。この実験の結果を図２３のグラフに示す。この図２３においては空燃比を、「○」「×」「＋」「△」「□」の順に大きくしており、例えば「○」はＡ／Ｆ＝１２．５、「×」はＡ／Ｆ＝１３．５、「＋」はＡ／Ｆ＝１４．６、「△」はＡ／Ｆ＝１５．５、「□」はＡ／Ｆ＝１６．５、である。
【０１４５】
図２３に表れているように、空燃比に基づく補正係数Ｑall／ｎairを熱発生率最大時燃料密度ρfuel@dQpeakに乗算すれば、この補正した燃料密度（ρfuel@dQpeak×Ｑall／ｎair）と熱発生率傾きｂ／ａとの関係は、空燃比が異なっていても概ね１本の直線として表すことができる。つまり、熱発生率傾きｂ／ａに対する空燃比の影響は、前記のように補正係数Ｑall／ｎairを乗算することによって好適に反映させることができる。
【０１４６】
補正係数Ｑall／ｎairは総熱発生量Ｑallを気筒内の空気量ｎair［ｍｏｌ］で除算したものであり、総熱発生量Ｑallは筒内燃料量ｎfuel［ｍｏｌ］×燃焼効率ｋとして表されるから、Ｑall／ｎair＝筒内燃料量ｎfuel×燃焼効率ｋ／筒内空気量ｎair＝ｋ／ＡＦであって、空燃比に基づくものである。
【０１４７】
この空燃比に基づく補正係数Ｑall／ｎairを上述した式（１３）の右辺に乗算して、空燃比の影響を反映させた以下の式（１７）を得る。なお、変数Ｑall，ｎairについては指数関数Ｑallζ，ｎairζとして、実験結果との合わせ込みを行ってもよい。また、そうして空燃比に基づく補正係数を乗算する代わりに、空燃比の影響については実験によって作成した補正マップによって反映させるようにしてもよい。
【０１４８】
【数１７】
2015162983.tif 000019
【０１４９】
なお、前記の式（１７）において、ρfuel@dQpeak＝筒内燃料量ｎfuel／Ｖ@dQpeakと表せば、以下の式（１８）が得られる。そして、この式（１８）において、ｎfuel／ｎair＝１／ＡＦと表せば、以下の式（１９）が得られる。
【０１５０】
【数１８】
2015162983.tif 000020
【０１５１】
【数１９】
2015162983.tif 000021
【０１５２】
続いて図２４には、ＥＧＲ率の変化による熱発生率波形への影響を示す。この図２４は、図１１と同じく、熱発生率最大時期ｄＱｐｅａｋＡが互いに一致するように点火時期ＳＡを調整しながら、ＥＧＲ率の変化によって得られる熱発生率波形を重ねて表示したもので、例えば図中のＥＧＲ０は０％であり、ＥＧＲ１は２０％である。
【０１５３】
この図２４に表れているように、熱発生率傾きｂ／ａはＥＧＲ率の影響を受けていて、ＥＧＲ率が大きいほど熱発生率傾きｂ／ａは小さくなっている。つまり、ＥＧＲ率が大きいほど燃焼速度が低下している。このように熱発生率傾きｂ／ａがＥＧＲ率の影響を受ける要因としては、燃料と酸素との反応（燃焼）に寄与しないＥＧＲガスの分子が火炎面の伝播を邪魔することが考えられる。
【０１５４】
そのようにＥＧＲ率の変化に対応して熱発生率傾きｂ／ａが変化する様子を調べた実験の結果を図２５に示す。この図２５においてはＥＧＲ率を、「○」「×」「＋」「△」「□」の順に大きくしており、例えば「○」はＥＧＲ率ｅｋegr＝０％、「×」はｅｋegr＝１０％、「＋」はｅｋegr＝１５％、「△」はｅｋegr＝２０％、「□」はｅｋegr＝２５％である。
【０１５５】
この図２５では、ＥＧＲ率に基づく補正係数（100−ekegr）／100を熱発生率最大時燃料密度ρfuel@dQpeakに乗算して補正しており、この補正後の燃料密度（ρfuel@dQpeak×（100−ekegr）／100と熱発生率傾きｂ／ａとの関係は、ＥＧＲ率が異なっていても概ね１本の直線として表すことができる。つまり、熱発生率傾きｂ／ａに対するＥＧＲ率の影響は、前記のように補正係数（100−ekegr）／100を乗算することによって、好適に反映させることができる。
【０１５６】
よって、ＥＧＲ率に基づく補正係数（100−ekegr）／100を上述した式（１３）の右辺に乗算して、ＥＧＲ率の影響を反映させた以下の式（２０）を得る。なお、この式（２０）においてもｅｋegrを指数関数ｅｋegrξとして、実験結果との合わせ込みを行ってもよい。また、そうしてＥＧＲ率に基づく補正係数（100−ekegr）／100を乗算する代わりに、ＥＧＲ率の影響については実験によって作成した補正マップによって反映させるようにしてもよい。
【０１５７】
【数２０】
2015162983.tif 000022
【０１５８】
以上、説明した空燃比およびＥＧＲ率に基づく熱発生率傾きｂ／ａの補正についてまとめると、前記式（１９）および式（２０）から以下の式（２１）が得られる。
【０１５９】
【数２１】
2015162983.tif 000023
【０１６０】
この式（２１）において、空燃比に基づく補正係数Ｑall／ｎairとＥＧＲ率に基づく補正係数（100−ekegr）／100とをまとめて、以下の式（２２）のように表してもよい。
【０１６１】
【数２２】
2015162983.tif 000024
【０１６２】
Ｇは、気筒内の空気量ｎair［ｍｏｌ］とＥＧＲガス量ｎEGR［ｍｏｌ］とを合算したものであり、Ｆは、筒内燃料量ｎfuel［ｍｏｌ］であるから、Ｇ／Ｆ＝（ｎair＋ｎEGR）／ｎfuelである。この補正係数Ｇ／Ｆは、余剰な空気中の窒素および酸素とＥＧＲガスとが火炎面の伝播を邪魔することの影響を表している。
【０１６３】
以上、説明したように、基本的には熱発生率最大時燃料密度ρfuel@dQpeakに基づいて熱発生率傾きｂ／ａを算出しつつ、必要に応じて燃焼効率、エンジン回転速度、熱発生率最大時筒内体積、空燃比、ＥＧＲ率などに基づく補正を行うことによって、それらの影響を好適に反映させることができ、熱発生率傾きｂ／ａをより高い精度で推定することが可能になる。
【０１６４】
−熱発生量推定部−
熱発生量推定部５は、前述した如く、燃焼期間の全期間において発生した熱発生量（総熱発生量Ｑall）を推定する部分である。
【０１６５】
以下、この熱発生量推定部における総熱発生量Ｑallの推定動作について説明する。
【０１６６】
まず、前半燃焼期間ａにおける熱発生量Ｑ１は以下の式（２３）によって算出される。
【０１６７】
【数２３】
2015162983.tif 000025
【０１６８】
一方、上述したように総熱発生量Ｑallは、筒内燃料量×ｋ（燃焼効率）として表すことができ（この式が熱発生量推定モデルに相当する）、燃焼効率ｋは、例えば暖機運転中のように油水温が低い場合に低下すると共に、負荷率やエンジン回転速度などの変化によっても変化する。そこで、本実施形態では、種々のエンジンの実験結果のデータベースを利用して燃焼効率ｋの値を、油水温、負荷率およびエンジン回転速度に対応づけて決定するマップを予め決定しておく。そして、この燃焼効率ｋの値を用いて、筒内燃料量から総熱発生量Ｑallが算出される。
【０１６９】
図２を参照して上述したように、熱発生率波形を作成するためには、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂおよび後半燃焼期間ｃを求める必要がある。この後半燃焼期間ｃにおける熱発生量Ｑ２は以下の式（２４）によって算出される。
【０１７０】
【数２４】
2015162983.tif 000026
【０１７１】
そして、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂは以下の式（２５）により、また、後半燃焼期間ｃは以下の式（２６）によりそれぞれ求められる。
【０１７２】
【数２５】
2015162983.tif 000027
【０１７３】
【数２６】
2015162983.tif 000028
【０１７４】
以上のようにして熱発生率波形算出装置１では、着火遅れ推定モデルを使用した着火遅れ推定部２での着火遅れ期間τの推定、前半燃焼期間推定モデルを使用した前半燃焼期間推定部３での前半燃焼期間ａの推定、熱発生率傾き推定モデルを使用した熱発生率傾き推定部４での熱発生率傾きｂ／ａの推定、熱発生量推定モデルを使用した熱発生量推定部５での熱発生量Ｑallの推定が行われると共に、最大熱発生率ｂおよび後半燃焼期間ｃが算出される。そして、これらの算出値を用いて、図２のように熱発生率波形を近似する三角波形を作成し、この三角波形を熱発生率波形として出力する。この出力された熱発生率波形を利用して、エンジン設計に際してのシステムの検討、制御の検討、適合値の検討が行われることになる。
【０１７５】
以上、説明したように本実施形態では、エンジンの熱発生率波形に近似した三角波形を作成する際に、その特性値の一つとして前半燃焼期間ａを用いている。そして、この前半燃焼期間ａをエンジンの負荷率、ＥＧＲ率、空燃比および油水温のいずれにも依らないものとし、熱発生率最大時筒内体積Ｖ@dQpeakとエンジン回転速度Ｎｅとに基づいて算出するようにしている。このため、前半燃焼期間ａの決定にかかる工数が大幅に低減できる。
【０１７６】
また、本実施形態では、同じく熱発生率波形（三角波形）の特性値の一つとして、着火時期から熱発生率最大時期ｄＱｐｅａｋＡまでの熱発生率傾きｂ／ａを用い、この熱発生率傾きｂ／ａを熱発生率最大時燃料密度ρfuel@dQpeakから算出するようにしている。このため、エンジンの負荷率および点火時期の両方に基づいて算出するのに比べて、熱発生率傾きｂ／ａの決定にかかる工数が低減できる。
【０１７７】
そして、それら前半燃焼期間ａや熱発生率傾きｂ／ａの算出値を用いて熱発生率波形（三角波形）を作成することにより、従来技術（例えば、ウィーベ関数を利用し、複数のパラメータを種々の運転条件に対して同定するもの）に比べて、十分な精度を保証しながらも工数を大幅に低減できる。よって、この熱発生率波形（三角波形）を利用して、エンジン設計に際しての種々の検討などを効率的に行うことができ、開発コストを削減できる。
【０１７８】
また、前記のように熱発生率最大時筒内体積Ｖ@dQpeakおよびエンジン回転速度Ｎｅから前半燃焼期間ａを算出し、熱発生率最大時燃料密度ρfuel@dQpeakから熱発生率傾きｂ／ａを算出して、熱発生率波形を作成していることから、この熱発生率波形は、気筒内での燃焼状態における物理現象に従って作成されたものとなっている。この点で、単に波形形状を合わせ込むために、形状パラメータなどの各種パラメータを数学的に合わせ込むといったウィーベ関数を利用した熱発生率波形の作成手法に比べて、本実施形態に係る熱発生率波形算出装置１によって算出された熱発生率波形は高い信頼性を得ることができる。
【０１７９】
さらに、本実施形態では、熱発生率波形の全体を作成しなくても、前記のように熱発生率最大時筒内体積Ｖ@dQpeakとエンジン回転速度Ｎｅとから前半燃焼期間ａを算出することができ、また、熱発生率最大時燃料密度ρfuel@dQpeakから熱発生率傾きｂ／ａを算出することができる。このため、混合気の着火後の初期の燃焼速度（熱発生率傾き）や前半燃焼期間ａについて所要の精度を確保しながら、従来よりも簡易に推定し、評価することができる。
【０１８０】
また、前述した如く、本実施形態で推定される前半燃焼期間ａと熱発生率傾きｂ／ａとは、互いに独立した（従属関係にない）値として推定されるものとなっている。このため、前半燃焼期間ａは、熱発生率傾きｂ／ａに推定誤差が含まれる場合にその誤差の影響を受けることの無い値として推定され、熱発生率傾きｂ／ａは、前半燃焼期間ａに推定誤差が含まれる場合にその誤差の影響を受けることの無い値として推定される。その結果、これら値の推定精度を高く確保することができる。
【０１８１】
−他の実施形態−
以上、説明した実施形態は、自動車用のガソリンエンジンを対象とした熱発生率波形を算出する熱発生率波形算出装置に本発明を適用した場合について説明した。本発明はこれに限らず、自動車用以外の火花点火機関に対しても適用が可能である。また、ガソリンエンジンにも特に限定されるものではなく、例えばガスエンジンに対しても適用が可能である。
【０１８２】
また、前記の実施形態で説明した熱発生率波形算出装置において実施される熱発生率波形算出方法も本発明の技術的思想の範疇である。
【０１８３】
前記の実施形態では、例えば式（４）のように、熱発生率最大時燃料密度ρfuel@dQpeakに基づいて熱発生率傾きｂ／ａを算出するようにしていた。本発明はこれに限らず、予め設定した所定時期（例えば熱発生率が最大熱発生率の所定割合に達する時期）における燃料密度に基づいて熱発生率傾きｂ／ａを算出するようにしてもよい。この場合、予め設定された所定の熱発生率に応じて、前記式（４）における各係数それぞれは実験等に基づいて同定されることになる。
【０１８４】
すなわち、例えば、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの間で熱発生率が増加している期間（熱発生率増加期間）において、熱発生率が最大熱発生率の８０％（適宜設定が可能）に達した時期における燃料密度に基づいて、熱発生率傾きｂ／ａを算出するようにしてもよい。但し、熱発生率最大時燃料密度ρfuel@dQpeakと熱発生率傾きｂ／ａとの相関が最も高いので、前記「予め設定された所定の熱発生率に達する時期」としては熱発生率最大時が好ましい。
【０１８５】
また、前記の実施形態では、基本的に熱発生率最大時燃料密度ρfuel@dQpeakに基づいて算出した熱発生率傾きｂ／ａを燃焼効率、エンジン回転速度、熱発生率最大時筒内体積、空燃比、ＥＧＲ率などに基づいて補正するようにしているが、これらの補正はいずれも必要に応じて行えばよく、補正は一切、しなくてもよい。
【０１８６】
さらに、前記の実施形態における混合気の着火遅れ期間τや前半燃焼期間ａなどの算出の仕方も一例に過ぎず、これらに限定されることはない。例えば、着火遅れ期間τや前半燃焼期間ａなどは実験、シミュレーションによって決定するようにしてもよい。
【０１８７】
また、前記の実施形態に係る熱発生率波形算出装置１は、三角波形を出力するものであった。本発明はこれに限らず、作成した三角波形に対して所定のフィルタ処理を行って熱発生率波形を作成し、これを出力するものであってもよい。
【産業上の利用可能性】
【０１８８】
本発明によると、火花点火式の内燃機関における熱発生率波形の作成にかかる工数を低減でき、コストの削減が可能になるので、例えば自動車用の内燃機関に適用可能である。
【符号の説明】
【０１８９】
１            熱発生率波形算出装置
ａ            前半燃焼期間（着火時期から熱発生率最大時期までの期間）
ｂ            最大熱発生率（熱発生率最大時期における熱発生率）
ｂ／ａ        熱発生率傾き
ａ＋ｃ        燃焼期間（混合気の着火から燃焼終了までのクランク角度期間）
ｄＱｐｅａｋＡ    熱発生率最大時期
Ｖ@dQpeak         熱発生率最大時筒内体積
ρfuel@dQpeak     熱発生率最大時燃料密度
ｋ            燃焼効率に基づく補正係数
Ｎｅε         エンジン回転速度に基づく補正係数
Ｖ@dQpeakη      熱発生率最大時筒内体積に基づく補正係数
Ｑall／ｎair    空燃比に基づく補正係数
(100-ｅｋegr)/100   ＥＧＲ率に基づく補正係数
【図１】
2015162983.tif 000029
【図２】
2015162983.tif 000030
【図３】
2015162983.tif 000031
【図４】
2015162983.tif 000032
【図５】
2015162983.tif 000033
【図６】
2015162983.tif 000034
【図７】
2015162983.tif 000035
【図８】
2015162983.tif 000036
【図９】
2015162983.tif 000037
【図１０】
2015162983.tif 000038
【図１１】
2015162983.tif 000039
【図１２】
2015162983.tif 000040
【図１３】
2015162983.tif 000041
【図１４】
2015162983.tif 000042
【図１５】
2015162983.tif 000043
【図１６】
2015162983.tif 000044
【図１７】
2015162983.tif 000045
【図１８】
2015162983.tif 000046
【図１９】
2015162983.tif 000047
【図２０】
2015162983.tif 000048
【図２１】
2015162983.tif 000049
【図２２】
2015162983.tif 000050
【図２３】
2015162983.tif 000051
【図２４】
2015162983.tif 000052
【図２５】
2015162983.tif 000053
【図２６】
2015162983.tif 000054
【図２７】
2015162983.tif 000055
【国際調査報告】
2015162983.tif 000056
2015162983.tif 000057
2015162983.tif 000058
2015162983.tif 000059