JP6256599B2 - 内燃機関の熱発生率波形算出装置および熱発生率波形算出方法 - Google Patents

Description

本発明は、火花点火式の内燃機関における熱発生率波形を算出するための装置、およびその算出方法に関連し、特に、混合気が着火してから熱発生率が最大になるまでの期間（本明細書では、この期間を「前半燃焼期間」と呼ぶ）に着目して、熱発生率波形を得るための技術に係る。

従来より、内燃機関の燃焼状態を表現するために、気筒内の熱発生率をウィーベ（Wiebe）関数によって近似することが行われている。ウィーベ関数は、複数のパラメータを特定することによって、熱発生率波形を好適に表現することができるものであって、内燃機関の燃焼による熱発生率や質量燃焼割合などを推定するために利用されている。

例えば、特許文献１に記載されているウィーベ関数パラメータの決定方法では、熱発生率が最大となるクランク角度での燃焼割合に基づいて、所定の式によりウィーベ関数の形状パラメータｍを同定するようにしている。そして、ｋ、ａ／θ_p ^m+1、θ_bといった他のパラメータについてもそれぞれ所定の式により同定して、実際の熱発生パターンに高い精度で適合するように、ウィーベ関数を決定することができる。

前記特許文献１には、そのように複数のパラメータｍ，ｋ，ａ／θ_p ^m+1，θ_bを同定し、ウィーベ関数を決定する作業を種々の運転条件に対して行うことにより、これらのパラメータと内燃機関の運転パラメータ（負荷率、回転速度、空燃比、点火時期など）との関係を把握することができる、と記載されている。また、そうして把握した関係を利用すれば、内燃機関のあらゆる運転条件に対してウィーベ関数を決定し、当該内燃機関の燃焼状態を精度良く表すことができると記載されている。

特開２００７−１７７６５４号公報

しかしながら前記特許文献１には、ウィーベ関数のパラメータｍ、ｋ、ａ／θ_p ^m+1、θ_bと内燃機関の運転パラメータとの関係を特定する具体的な方法について開示されていない。このため、実際には殆ど全ての運転条件についてパラメータｍ、ｋ、ａ／θ_p ^m+1、θ_bを同定し、それぞれの運転条件についてウィーベ関数を決定しなくてはならない。つまり、前記従来の方法には、熱発生率波形を作成するための工数をさらに低減し、コストを削減する余地が残されている。

また、前記の方法では、パラメータｍ、ｋ、ａ／θ_p ^m+1、θ_bをそれぞれ同定し、ウィーベ関数を決定して初めて熱発生率波形全体を表現することができ、これに基づいて燃焼状態を評価できるものである。よって、熱発生率波形全体を表現することなく、例えば前半燃焼期間（混合気が着火してから熱発生率が最大になるまでの期間）だけを簡易な方法で推定し、評価することはできない。

本発明はかかる諸点に鑑みてなされたものであり、その目的とするところは、燃焼状態を表す指標の一つである前半燃焼期間に着目し、熱発生率波形を作成（算出）するための工数を低減すると共に、例えば前半燃焼期間について、所要の精度を確保しながら簡易に推定し、評価できるようにすることにある。

−発明の解決原理−
本発明の発明者は、火花点火によって混合気が着火してから熱発生率が最大になるまでの期間である前半燃焼期間が、機関負荷率、ＥＧＲ率、空燃比および油水温の影響を受けず、熱発生率が最大となる時期における物理量と高い相関を有しているという新規な知見を得た。

かかる新規な知見に基づいて、本発明の解決原理は、熱発生率波形の特性値の一つとして前記前半燃焼期間を用い、この前半燃焼期間を、熱発生率が最大となる時期における物理量に基づいて推定するようにしたことにある。

−解決手段−
具体的に、本発明は、火花点火式の内燃機関における熱発生率波形を算出するための装置を対象とするもので、混合気の燃焼期間のうち着火時期から熱発生率が最大になる熱発生率最大時期までの期間を、前記熱発生率波形の特性値の一つである前半燃焼期間として規定する。そして、前記前半燃焼期間を、前記熱発生率最大時期における物理量に基づいて推定し、この推定した前半燃焼期間を用いて前記熱発生率波形を算出する構成とし、且つ、前記前半燃焼期間を推定するための前記熱発生率最大時期における物理量は、熱発生率最大時期における筒内体積および機関回転速度であって、下記の式（３）により、これら熱発生率最大時期における筒内体積および機関回転速度によって前記前半燃焼期間を推定する構成とする。

（ａは前半燃焼期間、Ｖ _@dQpeak は熱発生率最大時期における筒内体積、Ｎｅは機関回転速度、Ｃおよびαは実験に基づいた同定によって得られる係数、βは気筒内でのタンブル比に応じタンブル比が大きいほど大きな値として与えられる係数）

前記の特定事項によると、内燃機関において気筒内の混合気の燃焼による熱発生率の波形を算出するに際して、その熱発生率波形の特性値の一つとして、混合気の着火時期から熱発生率が最大になる熱発生率最大時期までの期間である前半燃焼期間を用いている。そして、本発明の発明者は、この前半燃焼期間が、内燃機関の負荷率、ＥＧＲ率、空燃比および油水温の運転条件の影響を受けず、熱発生率最大時期における物理量の影響を受けることを新たな知見として見出した。このため、この前半燃焼期間に基づいて熱発生率波形を算出するようにすれば、前記各運転条件を考慮する必要がなくなる。従って、負荷率、ＥＧＲ率、空燃比および油水温の運転条件を考慮して（これら運転条件を変化させて）熱発生率波形を作成するといった従来技術（前記ウィーベ関数を利用し、複数のパラメータ（形状パラメータなど）を種々の運転条件に対して同定するもの）に比べて、十分な精度を保証しながらも熱発生率波形を作成する工数を低減できる。

しかも、熱発生率波形全体を作成することなく、前記のように熱発生率最大時期における物理量に基づいて前半燃焼期間だけを推定できるので、所要の精度を確保しながら簡易に前半燃焼期間を推定し、その評価を行うことができる。

前半燃焼期間は気筒内の乱れの影響を大きく受けると想定される。つまり、気筒内の乱れが強いほど火炎伝播が急速になって、前半燃焼期間は短くなると想定される。そして、気筒内の乱れは、筒内体積に応じて変化する。つまり、前記熱発生率最大時期における筒内体積が大きいほど（燃焼行程においてピストンが下死点に近い位置にあるほど）気筒内の乱れは弱くなる。そして、気筒内の乱れが弱くなると、火炎伝播が緩慢になって前半燃焼期間は長くなる。このため、前記熱発生率最大時期が、ピストンが圧縮上死点に達するタイミング（ＴＤＣ）よりも遅角側にある場合、熱発生率最大時期における筒内体積が大きいほど気筒内の乱れは弱くなり前半燃焼期間は長くなる。一方、熱発生率最大時期における筒内体積が小さいほど気筒内の乱れは強くなり、火炎伝播が急速になって前半燃焼期間は短くなる。このように熱発生率最大時期における筒内体積は、気筒内の乱れに相関のあるパラメータとなっている。このため、本解決手段では、前記前半燃焼期間を、前記熱発生率最大時期における筒内体積に基づいて算出するようにしている。これにより、気筒内の乱れの影響を反映した前半燃焼期間を算出することができ、前半燃焼期間の推定精度が十分に確保され、この前半燃焼期間を用いて算出された熱発生率波形の信頼性も十分に確保されることになる。

気筒内の乱れは機関回転速度によっても変化する。つまり、機関回転速度が低いほど吸気系から気筒内に流入する空気の流速が低くなって気筒内の乱れは弱くなる。そして、気筒内の乱れが弱くなると、火炎伝播が緩慢になって前半燃焼期間は長くなる。逆に、機関回転速度が高いほど吸気系から気筒内に流入する空気の流速が高くなって気筒内の乱れは強くなる。そして、気筒内の乱れが強くなると、火炎伝播が急速になって前半燃焼期間は短くなる。このように機関回転速度も、気筒内の乱れに相関のあるパラメータとなっている。このため、本解決手段では、前記前半燃焼期間を、機関回転速度に基づく補正係数（機関回転速度の指数関数）を乗算して算出するようにしている。これにより、気筒内の乱れの影響をよりいっそう反映した前半燃焼期間を算出することができ、前半燃焼期間の推定精度が十分に確保され、この前半燃焼期間を用いて算出された熱発生率波形の信頼性も十分に確保されることになる。

前記機関回転速度に基づく補正係数としては、タンブル比に応じた値を指数とする機関回転速度の指数関数としている。前記タンブル比は、機関回転速度と共に気筒内の乱れに大きく影響を及ぼすものであるため、機関回転速度に基づく補正係数を、タンブル比に応じた値を指数とする機関回転速度の指数関数とすることによって、気筒内の乱れの影響をよりいっそう反映した前半燃焼期間を算出することができ、前半燃焼期間の推定精度が十分に確保される。

より具体的な前半燃焼期間の算出手法としては、仮想の熱発生率最大時期を設定し、この仮想の熱発生率最大時期における物理量に従って求められた仮想の前半燃焼期間が、実際の着火時期から仮想の熱発生率最大時期までの期間に一致したか否かを、仮想の熱発生率最大時期を変更しながら繰り返し計算することが挙げられる。つまり、仮想の熱発生率最大時期を設定し、実際の点火時期に従って求められた着火時期と前記仮想の熱発生率最大時期との間の期間である仮想の前半燃焼期間と、前記仮想の熱発生率最大時期における物理量に基づいて推定された前半燃焼期間とを比較し、これらが一致した場合における前記推定された前半燃焼期間が真の前半燃焼期間として算出される。そして、この真の前半燃焼期間を用いて前記熱発生率波形を算出する。これによれば、仮想の熱発生率最大時期を真の熱発生率最大時期に近付けていくことが可能になり、前半燃焼期間を推定するための熱発生率最大時期を正確に求めることができ、前半燃焼期間を高い精度で算出することができる。

前記のようにして算出した前半燃焼期間を用いて算出する熱発生率波形としては、例えば、混合気の着火から燃焼終了までのクランク角度期間を底辺とし、熱発生率最大時期における熱発生率を頂点とする三角波形が挙げられる。この三角波形によって熱発生率波形を近似すれば、着火時期から熱発生率最大時期までの熱発生率を表す三角形の底辺が前記前半燃焼期間として規定されることになる。

見方を変えると本発明は、火花点火式の内燃機関における熱発生率波形を算出する方法に係る。すなわち、まず、混合気の燃焼期間のうち着火時期から熱発生率が最大になる熱発生率最大時期までの期間を、前記熱発生率波形の特性値の一つである前半燃焼期間として規定する。そして、この前半燃焼期間を、前記熱発生率最大時期における物理量に基づいて推定し、この推定した前半燃焼期間を用いて前記熱発生率波形を算出し、且つ、前記前半燃焼期間を推定するための前記熱発生率最大時期における物理量は、熱発生率最大時期における筒内体積および機関回転速度であって、下記の式（３）により、これら熱発生率最大時期における筒内体積および機関回転速度によって前記前半燃焼期間を推定するものである。

（ａは前半燃焼期間、Ｖ _@dQpeak は熱発生率最大時期における筒内体積、Ｎｅは機関回転速度、Ｃおよびαは実験に基づいた同定によって得られる係数、βは気筒内でのタンブル比に応じタンブル比が大きいほど大きな値として与えられる係数）

本発明によると、内燃機関における熱発生率波形の特性値の一つとして、混合気の着火時期から熱発生率が最大になる熱発生率最大時期までの期間である前半燃焼期間を用い、この前半燃焼期間を熱発生率最大時期における物理量に基づいて推定するようにしたので、熱発生率波形の作成にかかる工数を低減できると共に、熱発生率波形の全体を作成することなく、前半燃焼期間について、所要の精度を確保しながら、簡易に推定し、評価できるようになる。

実施形態に係る熱発生率波形算出装置の構成、および、この熱発生率波形算出装置の入出力情報を表す図である。 熱発生率波形算出装置から出力される熱発生率波形の一例を示す図である。 熱発生率波形算出装置において行われる熱発生率波形の作成手順を示すフローチャート図である。 ＢＴＤＣ着火の場合における、点火時期ＳＡでの筒内の燃料密度ρ_fuel@SAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示す図である。 式（１）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示す図である。 ＡＴＤＣ着火の場合における、着火時期ＦＡでの筒内の燃料密度ρ_fuel@FAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示す図である。 式（２）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示す図である。 ＢＴＤＣ着火する場合の点火時期ＳＡと熱発生率波形とを示す図である。 ＡＴＤＣ着火する場合の点火時期ＳＡと熱発生率波形とを示し、図９（ａ）は点火時期ＳＡがＢＴＤＣである場合を、図９（ｂ）は点火時期ＳＡがＡＴＤＣである場合をそれぞれ示す図である。 負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。 ＥＧＲ率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。 空燃比のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。 油水温のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。 点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示した図である。 エンジン回転速度Ｎｅのみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。 あるエンジンに対して、式（３）で算出された予測前半燃焼期間と、実機において計測された実測前半燃焼期間との関係を検証した結果を示す図である。 他のエンジンに対して、式（３）で算出された予測前半燃焼期間と、実機において計測された実測前半燃焼期間との関係を検証した結果を示す図である。 負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示した図である。 点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示した図である。 熱発生率最大時燃料密度ρ_fuel@dQpeakと熱発生率傾きｂ／ａとの関係を、互いに異なるエンジン回転速度Ｎｅそれぞれに対して調べた実験の結果を示す図である。

以下、本発明の実施の形態を図面に基づいて説明する。本実施形態では、自動車用のガソリンエンジン（火花点火機関）を対象とした熱発生率波形を算出（作成）する熱発生率波形算出装置に本発明を適用した場合について説明する。

図１は、本実施形態に係る熱発生率波形算出装置１の構成、および、この熱発生率波形算出装置１の入出力情報を表す図である。この熱発生率波形算出装置１は、エンジンの状態量、制御パラメータの制御量および物理量の各種情報が入力される。これら入力情報としては、エンジン回転速度、負荷率、点火時期、ＥＧＲ率、空燃比、油水温、吸排気の各バルブの開閉タイミング（バルブタイミング）等が挙げられる。また、熱発生率波形算出装置１は、各入力情報に基づいて以下の各推定モデルを格納した推定部２〜５により熱発生率波形の各種特性値を推定して、この各種特性値を利用して作成した熱発生率波形を出力する。

−熱発生率波形の各特性値の推定部−
熱発生率波形算出装置１は、熱発生率波形の特性値として着火遅れ、前半燃焼期間、熱発生率傾きおよび熱発生量をそれぞれ推定するために、着火遅れ推定モデルを格納した着火遅れ推定部２、前半燃焼期間推定モデルを格納した前半燃焼期間推定部３、熱発生率傾き推定モデルを格納した熱発生率傾き推定部４、および、熱発生量推定モデルを格納した熱発生量推定部５を備えている。

着火遅れ推定部２は、エンジンの点火プラグ（点火栓）によって混合気に点火された時期（以下、点火時期という）、即ち点火プラグの電極間に火花放電が行われてから、この火花によって混合気に着火し初期火炎核が形成される時期（以下、着火時期という）までの期間（以下、着火遅れ期間という）を、着火遅れ推定モデルを使用して推定する部分である。この着火遅れ期間はクランク角度［ＣＡ］で表される。なお、前記着火時期の定義として、本実施形態では、点火時期以降に熱発生率（クランクシャフトの回転の単位クランク角度当たりの熱発生量）が１［Ｊ／ＣＡ］に達した時期とする。この値はこれに限らず適宜設定可能であり、例えば、点火時期以降の熱発生量が総熱発生量に対して所定割合（例えば５％）に達した時期を着火時期としてもよい。さらに、総熱発生量に対する熱発生量の割合が所定値に達した時期（例えば１０％に達した時期でのクランク角度位置）と前記熱発生量の割合が他の所定値に達した時期（例えば５０％に達した時期でのクランク角度位置）とに基づいて着火時期を定義するようにしてもよい。つまり、これらのクランク角度位置と熱発生量の割合とによって、熱発生率が上昇していく期間での熱発生率波形に近似した三角形（三角波形）を作成し、この三角波形に基づいて着火時期を定義するものである。また、熱発生率が上昇していく期間での一般的な熱発生率波形の形状を、前記クランク角度位置と熱発生量の割合との関係が成立するように当て嵌めて熱発生率波形を作成し、この熱発生率波形に基づいて着火時期を定義するようにしてもよい。前記各値はこれに限定されることなく適宜設定可能である。

前半燃焼期間推定部３は、混合気の燃焼期間のうち、前記着火時期から、火炎核の成長に伴って熱発生率が最大となるタイミング（点火時期から燃焼終了時期までの期間中で熱発生率が最大となるタイミング）までの期間である前半燃焼期間を、前半燃焼期間推定モデルを使用して推定する部分である。以下、この熱発生率が最大となるタイミングを熱発生率最大時期という。また、この熱発生率最大時期および前半燃焼期間はそれぞれクランク角度［ＣＡ］で表される。

熱発生率傾き推定部４は、前記前半燃焼期間、即ち着火時期から熱発生率最大時期までの期間における、クランク角度変化に対する平均的な熱発生率の増加率（熱発生率の傾き）を、熱発生率傾き推定モデルを使用して推定する部分である。すなわち、本実施形態では、図２を参照して以下に述べるように、熱発生率波形に近似した三角波形を作成するようにしており、熱発生率傾き推定部４は、その三角波形における着火時期から熱発生率最大時期までの熱発生率を表す斜辺の傾きを推定するものである。この熱発生率の傾きの単位としては［Ｊ／ＣＡ²］で表される。

熱発生量推定部５は、混合気の燃焼によって発生した熱発生量（燃焼期間の全期間において発生した熱発生量であって、点火時期から燃焼終了時期までの期間中における熱発生率の積算値）を、熱発生量推定モデルを使用して推定する部分である。この熱発生量の単位としては［Ｊ］で表される。

各推定部２〜５での推定動作によってそれぞれ、着火遅れ、前半燃焼期間、熱発生率傾き、および、熱発生量という熱発生率波形の特性値が求められ、これらの特性値を利用して熱発生率波形が作成される。そして、この作成された熱発生率波形が熱発生率波形算出装置１の出力となる。

このため、本実施形態に係る熱発生率波形算出装置１にあっては、図３に示すフローチャートの如く、着火遅れ推定部２における着火遅れ期間の推定動作（ステップＳＴ１）、前半燃焼期間推定部３における前半燃焼期間の推定動作（ステップＳＴ２）、熱発生率傾き推定部４における熱発生率傾きの推定動作（ステップＳＴ３）、熱発生量推定部５における熱発生量の推定動作（ステップＳＴ４）が順に行われた後、これら推定された特性値を利用した熱発生率波形の作成動作（ステップＳＴ５）が行われることになる。

図２は、各推定部２〜５で推定された特性値を利用して作成されて、熱発生率波形算出装置１から出力される熱発生率波形の一例を示している。この図２では、図中のタイミングＳＡが点火時期であり、図中のタイミングＦＡが着火時期である。このため、図中のτが着火遅れ期間となる。また、図中においてｄＱｐｅａｋＡが熱発生率最大時期であって、この熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率が図中のｂである。つまり、この熱発生率ｂが燃焼期間中における最大熱発生率となっている。また、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの期間である図中のａが前半燃焼期間となっている。このため、前半燃焼期間ａにおける熱発生率の傾きはｂ／ａとして表される。さらに、熱発生率最大時期ｄＱｐｅａｋＡから燃焼終了時期ＥＡまでの期間である図中のｃが後半燃焼期間となっている。また、図中のＱ１は前半燃焼期間ａにおける熱発生量であり、Ｑ２は後半燃焼期間ｃにおける熱発生量である。そして、燃焼期間の全期間において発生する熱発生量（総熱発生量Ｑ_all）は、これら熱発生量Ｑ１と熱発生量Ｑ２との和として表される。

言い換えると、本実施形態の熱発生率波形算出装置１は、混合気の着火から燃焼終了までのクランク角度期間（図中のＦＡからＥＡ）を底辺とし、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂを頂点とする三角波形によって熱発生率波形を近似するようにしている。この場合、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの熱発生率を表す三角形の底辺が前半燃焼期間ａとなる。本実施形態では、この熱発生率波形算出装置１の出力である熱発生率波形を利用し、エンジン設計に際してのシステムの検討、制御の検討、適合値の検討が行われる。

以下、各推定部２〜５での推定処理について具体的に説明する。

−着火遅れ推定部−
着火遅れ推定部２は、前述した如く、点火時期ＳＡから着火時期ＦＡまでの期間である着火遅れ期間τを推定する部分である。

この着火遅れ推定部２において行われる着火遅れ期間τの推定処理は以下のとおりである。

この着火遅れ期間τは、以下の式（１）および式（２）の何れか一方を利用して推定される（これらの式が着火遅れ推定モデルに相当する）。

ρ_fuel@SAは点火時期ＳＡにおける筒内の燃料密度（筒内燃料量［ｍｏｌ］／点火時期の筒内体積［Ｌ］）である。ρ_fuel@FAは着火時期ＦＡにおける筒内の燃料密度（筒内燃料量［ｍｏｌ］／着火時期の筒内体積［Ｌ］）である。Ｎｅはエンジン回転速度である。Ｃ₁，Ｃ₂，χ，δ，φ，ψはそれぞれ実験等に基づいて同定される係数である。

これら式（１）および式（２）は、空燃比が理論空燃比であり、ＥＧＲ率が「０」であり、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸気バルブの開閉タイミングが固定されていることを条件に成立する式となっている。

式（１）は、ピストンが圧縮上死点に達するタイミング（ＴＤＣ）よりも進角側（ＢＴＤＣ）で混合気が着火した（以下、ＢＴＤＣ着火という）場合の着火遅れ期間τの算出式である。また、式（２）は、ピストンが圧縮上死点に達したタイミング（ＴＤＣ）よりも遅角側（ＡＴＤＣ）で混合気が着火した（以下、ＡＴＤＣ着火という）場合の着火遅れ期間τの算出式である。

これらの式に示すように、着火遅れ期間τは、所定タイミングでの筒内の燃料密度ρ_fuelおよびエンジン回転速度Ｎｅを変数とする演算式によって算出される。

これら演算式によって着火遅れ期間τが算出可能である根拠について以下に説明する。

図４は、ＢＴＤＣ着火の場合における、点火時期ＳＡでの筒内の燃料密度ρ_fuel@SAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示すグラフである。この実験は、空燃比を理論空燃比とし、ＥＧＲ率を「０」とし、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸気バルブの開閉タイミングを固定して行われたものである。また、この図４では、「○」「△」「□」「◇」「×」「＋」「▽」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「△」は１０００ｒｐｍ、「□」は１２００ｒｐｍ、「◇」は１６００ｒｐｍ、「×」は２４００ｒｐｍ、「＋」は３２００ｒｐｍ、「▽」は３６００ｒｐｍである。

この図４に示すように、ＢＴＤＣ着火した場合、点火時期ＳＡにおける筒内の燃料密度ρ_fuel@SAと着火遅れ期間τとの間にはエンジン回転速度Ｎｅ毎に相関がある。つまり、これらの相関を概ね１本の曲線で表すことができる。図４では、エンジン回転速度Ｎｅが１０００ｒｐｍの場合および２４００ｒｐｍの場合のそれぞれについて、点火時期ＳＡにおける筒内の燃料密度ρ_fuel@SAと着火遅れ期間τとの相関を１本の曲線で表している。

図４に示すように、点火時期ＳＡにおける筒内の燃料密度ρ_fuel@SAが高いほど着火遅れ期間τは短くなっている。これは、燃料密度ρ_fuel@SAが高いほど、点火プラグ周辺の燃料分子の数が多くなり、点火プラグの点火後における火炎核の成長が急速に行われるためであると考えられる。また、エンジン回転速度Ｎｅは着火遅れ期間τに影響を与えている。つまり、エンジン回転速度Ｎｅが高いほど着火遅れ期間τは短くなっている。これは、エンジン回転速度Ｎｅが高いほど気筒内での混合気流れの乱れ（以下、単に乱れという）が強くなって前記火炎核の成長が急速に行われるためであると考えられる。このように、点火時期ＳＡにおける筒内の燃料密度ρ_fuel@SAおよびエンジン回転速度Ｎｅは着火遅れ期間τに影響を与えるパラメータとなっている。

図５は、式（１）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示すグラフである。この予測着火遅れ期間を求めるに当たっては、式（１）におけるＣ₁，χ，δの各係数をエンジン運転条件に応じて同定することにより得られた予測式を使用している。この図５では、「○」「△」「□」「◇」「×」「＋」「▽」「☆」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「△」は１０００ｒｐｍ、「□」は１２００ｒｐｍ、「◇」は１６００ｒｐｍ、「×」は２０００ｒｐｍ、「＋」は２４００ｒｐｍ、「▽」は３２００ｒｐｍ、「☆」は３６００ｒｐｍである。

この図５から明らかなように、予測着火遅れ期間は実測着火遅れ期間に略一致しており、式（１）によって、ＢＴＤＣ着火した場合の着火遅れ期間が高い精度で算出されていることが解る。

図６は、ＡＴＤＣ着火の場合における、着火時期ＦＡでの筒内の燃料密度ρ_fuel@FAの変化に対する着火遅れ期間τの変化を実験によって計測した結果を示すグラフである。この実験は、エンジン回転速度を固定し、空燃比を理論空燃比とし、ＥＧＲ率を「０」とし、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸気バルブの開閉タイミングを固定して行われたものである。また、この図６では、「○」「×」「＋」「△」の順で負荷率が高くなっている。例えば「○」は負荷率２０％、「×」は負荷率３０％、「＋」は負荷率４０％、「△」は負荷率５０％である。

この図６に示すように、ＡＴＤＣ着火した場合、着火時期ＦＡにおける筒内の燃料密度ρ_fuel@FAと着火遅れ期間τとの間には負荷率に依らず（負荷率に関わりなく）相関がある。つまり、これらの相関を概ね１本の曲線で表すことができる。

図６に示すように、着火時期ＦＡにおける筒内の燃料密度ρ_fuel@FAが高いほど着火遅れ期間τは短くなっている。これは、前述したように、燃料密度ρ_fuel@FAが高いほど、点火プラグ周辺の燃料分子の数が多くなり、点火プラグの点火後における火炎核の成長が急速に行われるためであると考えられる。このように、着火時期ＦＡにおける筒内の燃料密度ρ_fuel@FAは着火遅れ期間τに影響を与えるパラメータとなっている。また、前述した場合と同様に、エンジン回転速度Ｎｅも着火遅れ期間τに影響を与えるパラメータになっていると想定される。

図７は、式（２）で算出された予測着火遅れ期間と、実機において計測された実測着火遅れ期間との関係を検証した結果を示すグラフである。この予測着火遅れ期間を求めるに当たっては、式（２）におけるＣ₂，φ，ψの各係数をエンジン運転条件に応じて同定することにより得られた予測式を使用している。この図７では、「○」「×」「＋」「△」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「×」は１２００ｒｐｍ、「＋」は３６００ｒｐｍ、「△」は４８００ｒｐｍである。

この図７から明らかなように、予測着火遅れ期間は実測着火遅れ期間に略一致しており、式（２）によって、ＡＴＤＣ着火した場合の着火遅れ期間が高い精度で算出されていることが解る。

本発明の発明者は、これらの新たな知見に基づいて前記式（１）および式（２）を導き出した。

以下、着火遅れ期間τの算出に当たり、着火時期で層別する理由について説明する。つまり、ＢＴＤＣ着火した場合と、ＡＴＤＣ着火した場合とを層別し、それぞれに対して異なる演算式（前記式（１）および式（２））を用いて着火遅れ期間τを算出する理由について説明する。

まず、ＢＴＤＣ着火する場合、図８（点火時期ＳＡと熱発生率波形とを示す図）に示すように、点火時期ＳＡも、ピストンが圧縮上死点に達するタイミングよりも進角側（ＢＴＤＣ）である。この場合、点火時期ＳＡを迎えた後、ピストンは圧縮上死点に向かって移動することになる。つまり、筒内体積が小さくなっていき、これに伴って燃料密度ρ_fuelは大きくなっていく。このため、燃料密度ρ_fuelとしては、着火時期ＦＡでの燃料密度ρ_fuel@FAよりも点火時期ＳＡでの燃料密度ρ_fuel@SAの方が小さくなっている。そして、着火遅れ期間の最大値（想定される着火遅れ期間のうち最も長い着火遅れ期間）に相関のある値である点火時期ＳＡでの燃料密度ρ_fuel@SAに対して、予め同定された各種係数を乗算することにより、着火遅れ期間τを高い精度で求めることができる。

一方、ＡＴＤＣ着火する場合、図９（点火時期ＳＡと熱発生率波形とを示す図）に示すように、点火時期ＳＡは、ピストンが圧縮上死点に達するタイミングよりも進角側（ＢＴＤＣ）である場合（図９（ａ）を参照）と、遅角側（ＡＴＤＣ）である場合（図９（ｂ）を参照）とがある。これらの場合、着火時期ＦＡを迎えた後、ピストンは下死点に向かって移動することになる。つまり、筒内体積が大きくなっていき、これに伴って燃料密度ρ_fuelは小さくなっていく。このため、燃料密度ρ_fuelとしては、点火時期ＳＡでの燃料密度ρ_fuel@SAよりも着火時期ＦＡでの燃料密度ρ_fuel@FAの方が小さくなっている可能性が高い。そして、着火遅れ期間の最大値（想定される着火遅れ期間のうち最も長い着火遅れ期間）に相関のある値である着火時期ＦＡでの燃料密度ρ_fuel@FAに対して、予め同定された各種係数を乗算することにより、着火遅れ期間τを高い精度で求めることができる。

また、これら式（１）および式（２）の何れを使用するかを判断する手順（着火時期がＢＴＤＣ着火であるかＡＴＤＣ着火であるかを判断する手順）および着火遅れ期間（後述する真の着火遅れ期間）を算出する手順は以下のとおりである。まず、仮想の着火時期を設定し、この仮想の着火時期における筒内体積を求める。この筒内体積は仮想の着火時期に対応するクランク角度位置（ピストンの位置）によって幾何学的に求めることが可能であるので、仮想の着火時期から筒内体積は一義的に決まる。そして、この筒内体積と燃料噴射量とから燃料密度を求める。そして、仮想の着火時期をＢＴＤＣ着火として設定した場合には、この仮想の着火時期における燃料密度およびエンジン回転速度を式（１）に代入して推定着火遅れ期間を算出する。一方、仮想の着火時期をＡＴＤＣ着火として設定した場合には、この仮想の着火時期における燃料密度およびエンジン回転速度を式（２）に代入して推定着火遅れ期間を算出する。そして、仮想の着火時期に対し、前記算出された推定着火遅れ期間だけ進角した時期を仮想の点火時期とする。ここで、この仮想の点火時期と実際の点火時期（入力情報である点火時期）とを比較する。仮想の点火時期と実際の点火時期とが一致しない場合、前記仮想の着火時期を変更する。例えば仮想の着火時期を遅角側に変更する。そして、再度、この仮想の着火時期における燃料密度およびエンジン回転速度を式（１）または式（２）に代入して（仮想の着火時期をＢＴＤＣ着火として設定した場合には式（１）に代入し、仮想の着火時期をＡＴＤＣ着火として設定した場合には式（２）に代入する）推定着火遅れ期間を算出して、仮想の点火時期を求め、これと実際の点火時期（入力情報である点火時期）とを比較する。この動作を繰り返し、仮想の点火時期と実際の点火時期とが一致した場合における仮想の着火時期が真の着火時期として得られることになる。そして、この際（真の着火時期が得られた際）、式（１）または式（２）において算出された推定着火遅れ期間が真の着火遅れ期間として得られることになる。また、前記真の着火時期がＢＴＤＣであった場合（ＢＴＤＣ着火であった場合）に前記求められた着火時期を改めて式（１）に代入して着火遅れ期間τの算出を行い、前記真の着火時期がＡＴＤＣであった場合（ＡＴＤＣ着火であった場合）に前記求められた着火時期を改めて式（２）に代入して着火遅れ期間τの算出を行うようにしてもよい。

以上の動作を言い換えると以下のとおりとなる。前記実際の点火時期と仮想の着火時期との間の期間（仮想の着火時期で着火したとした場合の仮想の着火遅れ期間）と、式（１）または式（２）で算出された（推定された）推定着火遅れ期間とを比較し、これらが一致しない場合には、前記仮想の着火時期を変更する。そして、再度、式（１）または式（２）によって推定着火遅れ期間を算出した後、前記実際の点火時期と仮想の着火時期との間の期間（仮想の着火遅れ期間）と、式（１）または式（２）で算出された推定着火遅れ期間とを比較する。この動作を繰り返し、これらが一致した場合（仮想の着火遅れ期間と推定着火遅れ期間とが一致した場合）における推定着火遅れ期間が真の着火遅れ期間として得られることになる。

以上のような着火遅れ推定部２による着火遅れ期間τの推定により、エンジン運転領域の全領域に対して着火遅れ期間τの推定が可能になる。

以上のようにして着火遅れ期間τが算出されると、前記点火時期ＳＡに着火遅れ期間τを加算することで着火時期ＦＡを求めることができる。

−前半燃焼期間推定部−
前半燃焼期間推定部３は、前述した如く、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの期間である前半燃焼期間ａを推定する部分である。

この前半燃焼期間推定部３において行われる前半燃焼期間ａの推定処理は以下のとおりである。

この前半燃焼期間ａ［ＣＡ］は、以下の式（３）を利用して推定される（この式が前半燃焼期間推定モデルに相当する）。

Ｖ_@dQpeakは前記熱発生率最大時期ｄＱｐｅａｋＡにおける物理量としての筒内体積［Ｌ］であり、以下では、熱発生率最大時筒内体積ともいう。Ｎｅはエンジン回転速度（機関回転速度）である。

この式（３）は、吸気バルブの開閉タイミングが固定されていることを条件に成立する式となっている。また、この式（３）は、負荷率、ＥＧＲ率、空燃比、油水温の影響を受けることなく成立するものとなっている。つまり、式（３）は、前半燃焼期間ａが、負荷率、ＥＧＲ率、空燃比、油水温の影響を受けないことに基づいて成立している。

この式（３）によって前半燃焼期間ａが算出可能である根拠について以下に説明する。

図１０〜図１３それぞれは、互いに異なるエンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示したものである。図１０は、負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。図１１は、ＥＧＲ率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。図１２は、空燃比のみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。また、図１３は、エンジンの暖機運転の途中のように油水温のみが互いに異なる場合に得られる熱発生率波形を重ねて表示したものである。

これら図１０〜図１３に示すように、負荷率、ＥＧＲ率、空燃比、油水温の何れが変化しても前半燃焼期間ａは一定に維持されている。つまり、前半燃焼期間ａは、負荷率、ＥＧＲ率、空燃比、油水温の影響を受けないものであることが解る。

一方、図１４は、点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。この図１４から解るように、点火時期ＳＡが遅角されるほど前半燃焼期間ａは長くなっている。

また、図１５は、エンジン回転速度Ｎｅのみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡを互いに一致させるように点火時期ＳＡを調整した各熱発生率波形を重ねて表示したものである。エンジン回転速度Ｎｅが高いほど単位時間［ｍｓ］当たりにおけるクランクの回転角度［ＣＡ］は大きくなるので、その分、前半燃焼期間ａは長くなる（クランク角度軸上で長くなる）はずであるが、図１５に示すものでは、エンジン回転速度Ｎｅが異なっても前半燃焼期間ａは殆ど変化していない。これは、エンジン回転速度Ｎｅが高いほど前半燃焼期間ａが短くなる要因が存在しているものと考えられる。つまり、エンジン回転速度Ｎｅが高いほど単位時間当たりにおけるクランクの回転角度が大きくなることに起因して前半燃焼期間ａが長くなることとは別に、「他の要因」によって前半燃焼期間ａが短くなっていると想定できる。

このように前半燃焼期間ａは、点火時期ＳＡおよびエンジン回転速度Ｎｅの影響を受けるものであることが解る。

前半燃焼期間ａが、点火時期ＳＡおよびエンジン回転速度Ｎｅの影響を受ける要因としては、点火時期ＳＡおよびエンジン回転速度Ｎｅが、気筒内の乱れに影響を与えることが考えられる。

つまり、前記熱発生率最大時期ｄＱｐｅａｋＡがＴＤＣよりも遅角側にある場合を考えると、点火時期ＳＡが遅角側に移行するほど着火時期ＦＡおよび熱発生率最大時期ｄＱｐｅａｋＡも遅角側に移行し、この熱発生率最大時期ｄＱｐｅａｋＡにおける筒内体積（熱発生率最大時筒内体積Ｖ_@dQpeak）が大きくなると共に、気筒内の乱れは弱くなる。そして、気筒内の乱れが弱くなると、火炎伝播が緩慢になって前半燃焼期間ａは長くなる。逆に、点火時期ＳＡが進角側に移行するほど着火時期ＦＡおよび熱発生率最大時期ｄＱｐｅａｋＡも進角側に移行し、熱発生率最大時筒内体積Ｖ_@dQpeakが小さくなると共に、気筒内の乱れは強くなる。これにより火炎伝播が急速になって、前半燃焼期間ａは短くなる。

また、エンジン回転速度Ｎｅが低いほど吸気系から気筒内に流入する空気の流速が低くなって気筒内の乱れは弱くなる。そして、気筒内の乱れが弱くなると、火炎伝播が緩慢になって前半燃焼期間ａは長くなる。逆に、エンジン回転速度Ｎｅが高いほど吸気系から気筒内に流入する空気の流速が高くなって気筒内の乱れは強くなる。そして、気筒内の乱れが強くなると、火炎伝播が急速になって前半燃焼期間ａは短くなる。前述した「他の要因（前半燃焼期間ａを短くする要因）」は、このエンジン回転速度Ｎｅが高いほど気筒内の乱れが強くなることに起因して火炎伝播が急速になっていることである。

本発明の発明者は、この新たな知見に基づいて前記式（３）を導き出した。そして、この式（３）においては、制御量である点火時期ＳＡに相関のある物理量として筒内体積、特に、熱発生率最大時筒内体積Ｖ_@dQpeakを変数として用いている。つまり、前述した如く、点火時期ＳＡが遅角側に移行するほど熱発生率最大時期ｄＱｐｅａｋＡも遅角側に移行し、筒内体積Ｖ_@dQpeakが大きくなることから、点火時期ＳＡに相関のある物理量として、熱発生率最大時筒内体積Ｖ_@dQpeakを変数として用いている。

前記式（３）の変数である熱発生率最大時筒内体積Ｖ_@dQpeakを求める手順および前半燃焼期間ａを算出する手順は以下のとおりである。まず、仮想の熱発生率最大時期を設定し、この仮想の熱発生率最大時期における筒内体積を求める。この筒内体積は仮想の熱発生率最大時期に対応するクランク角度位置（ピストンの位置）によって幾何学的に求めることが可能であるので、仮想の熱発生率最大時期から筒内体積は一義的に決まる。そして、この仮想の熱発生率最大時期における筒内体積およびエンジン回転速度を式（３）に代入して推定前半燃焼期間を算出する。そして、前記仮想の熱発生率最大時期に対し、この算出した推定前半燃焼期間だけ進角した時期を仮想の着火時期とする。前述した着火遅れ推定部２では着火遅れ期間τが算出されているため、前記点火時期ＳＡに着火遅れ期間τを加算することで着火時期ＦＡが算出できる。ここで、前記仮想の着火時期と算出された着火時期ＦＡとを比較する。仮想の着火時期と算出された着火時期ＦＡとが一致しない場合、前記仮想の熱発生率最大時期を変更する。例えば仮想の熱発生率最大時期を遅角側に変更する。そして、再度、この仮想の熱発生率最大時期における筒内体積およびエンジン回転速度を式（３）に代入して推定前半燃焼期間を算出して、仮想の着火時期を求め、これと算出された着火時期ＦＡ（前記点火時期ＳＡに、着火遅れ推定部２で算出された着火遅れ期間τを加算することで得られる着火時期ＦＡ）とを比較する。この動作を繰り返し、仮想の着火時期と算出された着火時期ＦＡとが一致した場合における仮想の熱発生率最大時期が真の熱発生率最大時期ｄＱｐｅａｋＡとして得られることになる。そして、この際（真の熱発生率最大時期ｄＱｐｅａｋＡが得られた際）、式（３）において算出された推定前半燃焼期間が真の前半燃焼期間として得られることになる。また、前記真の熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率最大時筒内体積Ｖ_@dQpeakを幾何学的に求め、これを改めて式（３）に代入して前半燃焼期間ａの算出を行うようにしてもよい。

以上の動作を言い換えると以下のとおりとなる。前記着火時期ＦＡ（実際の点火時期に従って求められた着火時期）と仮想の熱発生率最大時期との間の期間（仮想の前半燃焼期間）と、式（３）で算出された（推定された）推定前半燃焼期間（仮想の熱発生率最大時期における物理量に基づいて推定された前半燃焼期間）とを比較し、これらが一致しない場合には、前記仮想の熱発生率最大時期を変更する。そして、再度、式（３）によって推定前半燃焼期間を算出した後、前記着火時期ＦＡと仮想の熱発生率最大時期との間の期間（仮想の前半燃焼期間）と、式（３）で算出された推定前半燃焼期間とを比較する。この動作を繰り返し、これらが一致した場合（仮想の前半燃焼期間と推定前半燃焼期間とが一致した場合）における推定前半燃焼期間が真の前半燃焼期間ａとして得られることになる。

式（３）における各係数について具体的に説明すると、Ｃおよびαは実験等に基づいた同定によって得られる。また、βは気筒内でのタンブル比に応じた値となっており、タンブル比が大きいほど大きな値として与えられる。なお、βは実験等に基づいた同定によって設定するようにしてもよい。また、これら係数は、吸気バルブの開閉タイミングの変化に対して同定することも可能である。このように、式（３）は、熱発生率最大時筒内体積Ｖ_@dQpeakに基づき、且つタンブル比に応じた値βを指数とするエンジン回転速度Ｎｅの指数関数（補正係数）を乗算することにより、前半燃焼期間ａを算出するものとなっている。

図１６および図１７は、互いに異なるエンジンに対して、式（３）で算出された予測前半燃焼期間と、実機において計測された実測前半燃焼期間との関係を検証した結果を示すグラフである。この予測前半燃焼期間を求めるに当たっては、式（３）における係数Ｃをエンジン運転条件に応じて同定することにより得られた予測式を使用している。図１６では、「○」「△」「□」「◇」「×」「＋」「▽」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「△」は１０００ｒｐｍ、「□」は１２００ｒｐｍ、「◇」は１６００ｒｐｍ、「×」は２４００ｒｐｍ、「＋」は３２００ｒｐｍ、「▽」は３６００ｒｐｍである。また、図１７では、「○」「×」「＋」「△」「□」の順でエンジン回転速度Ｎｅが高くなっている。例えば「○」は８００ｒｐｍ、「×」は１２００ｒｐｍ、「＋」は２４００ｒｐｍ、「△」は３６００ｒｐｍ、「□」は４８００ｒｐｍである。

これら図１６および図１７から明らかなように、予測前半燃焼期間は実測前半燃焼期間に略一致しており、式（３）によって前半燃焼期間ａが高い精度で算出されていることが解る。

以上、述べたように前半燃焼期間ａを、負荷率、空燃比、ＥＧＲ率、油水温の影響を受けないものとして、熱発生率最大時筒内体積Ｖ_@dQpeakとエンジン回転速度Ｎｅとに基づいて推定できる。これら熱発生率最大時筒内体積Ｖ_@dQpeakとエンジン回転速度Ｎｅは前述した如く気筒内の乱れに相関のあるパラメータとなっている。言い換えると、負荷率、空燃比、ＥＧＲ率、油水温は、気筒内の乱れに殆ど相関が無いことから前半燃焼期間ａに影響を与えないものと想定される。そして、これら負荷率、空燃比、ＥＧＲ率、油水温を考慮することなしに、気筒内の乱れに相関のあるパラメータである熱発生率最大時筒内体積Ｖ_@dQpeakとエンジン回転速度Ｎｅとに基づいて前半燃焼期間ａを推定できることから、エンジンの種々の運転条件における前半燃焼期間ａを決定するための工数が大幅に低減できることになる。

なお、前述したように、前半燃焼期間は負荷率の影響を受けない。この負荷率は、燃料噴射量を制御するためのパラメータの一つであり、燃料噴射量は気筒内の燃料密度に影響を及ぼす制御パラメータである。このため、前半燃焼期間は、気筒内の燃料密度に依ることなく推定される。具体的には、前述したように、熱発生率最大時筒内体積Ｖ_@dQpeakおよびエンジン回転速度Ｎｅといった気筒内の乱れに影響を与えるパラメータに基づいて前半燃焼期間は推定される。これに対し、後述するように、熱発生率傾きは気筒内の燃料密度に基づいて推定される。このように、本実施形態で推定される前半燃焼期間と熱発生率傾きとは、互いに独立した（従属関係にない）値として推定されるものとなっている。

−熱発生率傾き推定部−
熱発生率傾き推定部４は、前述した如く、前半燃焼期間ａにおける熱発生率の傾きｂ／ａ（以下、熱発生率傾きという）を推定する部分である。

この熱発生率傾き推定部４において行われる熱発生率傾きｂ／ａの推定処理は以下のとおりである。

この熱発生率傾きｂ／ａ［Ｊ／ＣＡ²］は、基本的には以下の式（４）を利用して推定される（この式が熱発生率傾き推定モデルに相当する）。

ρ_fuel@dQpeakは、前記熱発生率最大時期ｄＱｐｅａｋＡにおける燃料密度（筒内燃料量［ｍｏｌ］／熱発生率最大時期の筒内体積［Ｌ］）であり、以下では、熱発生率最大時燃料密度ともいう。Ｃ₃は実験等に基づいて同定される係数である。

この式（４）は、エンジン回転速度が固定されており、空燃比が理論空燃比であり、ＥＧＲ率が「０」であり、エンジンの暖機運転が完了しており（油水温が所定値以上であり）、吸気バルブの開閉タイミングが固定されていることを条件に成立する式となっている。なお、エンジン回転速度、空燃比、ＥＧＲ率、エンジンの油水温等の影響については後述する。

この式（４）によって熱発生率傾きｂ／ａが算出可能である根拠について以下に説明する。

図１８（ａ）〜図１８（ｄ）それぞれは、負荷率のみが互いに異なる各エンジン運転状態において得られる熱発生率波形であって、熱発生率最大時期ｄＱｐｅａｋＡが互いに一致するように点火時期ＳＡを調整した各熱発生率波形を重ねて表示したものである。図１８（ａ）〜図１８（ｄ）の順に点火時期は遅角側に変化しており、また、各図において負荷率は、ＫＬ１，ＫＬ２，ＫＬ３の順に大きくなっている。例えば図１８において、ＫＬ１は負荷率２０％、ＫＬ２は負荷率３０％，ＫＬ３は負荷率４０％である。

これら図１８（ａ）〜図１８（ｄ）に表れているように、熱発生率傾きｂ／ａは負荷率および点火時期ＳＡの影響を受けている。具体的には、それぞれ点火時期ＳＡの異なる図１８（ａ）〜図１８（ｄ）のいずれにおいても、負荷率が大きいほど熱発生率傾きｂ／ａは大きくなっている。このように熱発生率傾きｂ／ａが負荷率の影響を受ける要因として、負荷率に応じて気筒内の燃料密度が変化することが考えられる。つまり、負荷率が高いほど筒内燃料量が多くなるので、気筒内の燃料密度も高くなって、混合気の燃焼速度が高くなるものと考えられる。

また、図１８（ａ）〜図１８（ｄ）の順に点火時期ＳＡが遅角側に移行するに連れて、熱発生率傾きｂ／ａが小さくなっている。図１９（ａ）および図１９（ｂ）それぞれは、点火時期ＳＡの変化による影響を調べるために、点火時期ＳＡのみが互いに異なる各エンジン運転状態において得られる熱発生率波形を重ねて表示したものである。これら図１９（ａ）および図１９（ｂ）においてはそれぞれ負荷率が異なっているが、そのいずれにおいても熱発生率傾きｂ／ａは、点火時期ＳＡの遅角側への移行に連れて小さくなる傾向がある。

このように熱発生率傾きｂ／ａが点火時期ＳＡの影響を受ける要因も、前記した負荷率の場合と同じく気筒内の燃料密度によるものと考えられる。すなわち、ピストンが圧縮上死点（ＴＤＣ）付近にあるときは、クランク角度の変化に伴う筒内体積の変化は小さいが、膨張行程においてＴＤＣから離れてゆくに連れて（例えばＡＴＤＣ１０°ＣＡくらいからは）遅角側ほど筒内体積が大きくなってゆき、これに連れて気筒内の燃料密度は低下してゆく。

そして、前記図１９（ａ）および図１９（ｂ）に表れているように熱発生率波形は、点火時期ＳＡの遅角に伴い全体として遅角側に移動するとともに、着火時期ＦＡ（波形の始点）がＴＤＣ以降になると、徐々に熱発生率波形の傾きも小さくなってゆく。この結果、着火時期ＦＡ（波形の始点）から熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂ（波形の頂点）までを結んだ直線（図中には一点鎖線で示す）の傾き、即ち熱発生率傾きｂ／ａも遅角側に向かって徐々に小さくなっている。

そのように点火時期ＳＡの遅角（即ち着火時期ＦＡの遅角）が熱発生率傾きｂ／ａに及ぼす影響は、この熱発生率傾きｂ／ａと熱発生率最大時燃料密度ρ_fuel@dQpeakとの関係に顕著に表れている。すなわち、図１９（ａ）および図１９（ｂ）に表れているように、点火時期ＳＡの遅角に伴い熱発生率最大時期ｄＱｐｅａｋＡが遅角側に移行し、この熱発生率最大時期ｄＱｐｅａｋＡにおける筒内体積（熱発生率最大時筒内体積Ｖ_@dQpeak）が徐々に大きくなってゆくと、これに応じて熱発生率最大時燃料密度ρ_fuel@dQpeakが小さくなってゆき、これに対応して熱発生率傾きｂ／ａが小さくなってゆく。

本発明者は、そのように熱発生率最大時燃料密度ρ_fuel@dQpeakの変化に対応して熱発生率傾きｂ／ａが変化する様子を調べた。この実験の結果を図２０（ａ）〜図２０（ｄ）のグラフに示す。これらの各図において負荷率は、「○」「×」「＋」「△」「□」「◇」「▽」「☆」の順に大きくなっている。例えば図２０において、「○」は負荷率１５％、「×」は負荷率２０％、「＋」は負荷率２５％、「△」は負荷率３０％、「□」は負荷率３５％、「◇」は負荷率４０％、「▽」は負荷率４５％、「☆」は負荷率５０％である。

また、図２０（ａ）〜図２０（ｄ）の順にエンジン回転速度Ｎｅが高くなっており、例えば、図２０（ａ）は８００ｒｐｍ、図２０（ｂ）は１２００ｒｐｍ、図２０（ｃ）は２０００ｒｐｍ、図２０（ｄ）は３２００ｒｐｍである。

図２０（ａ）〜図２０（ｄ）それぞれに表れているように、エンジン回転速度を固定すれば、熱発生率最大時燃料密度ρ_fuel@dQpeakと熱発生率傾きｂ／ａとの関係は、負荷率や点火時期ＳＡが異なっていても概ね１本の直線として表すことができ、両者の間には高い相関（具体的には概略比例関係）のあることが解る。つまり、熱発生率傾きｂ／ａに対するエンジンの負荷率および点火時期ＳＡの影響は、熱発生率最大時燃料密度ρ_fuel@dQpeakという一つのパラメータによって、まとめて表現されるものである。

本発明の発明者は、この新たな知見に基づいて前記式（４）を導き出した。

前記式（４）の変数である熱発生率最大時燃料密度ρ_fuel@dQpeakは、前述したように筒内燃料量を熱発生率最大時筒内体積Ｖ_@dQpeakによって除算して、求めることができる。熱発生率最大時筒内体積Ｖ_@dQpeakを求める手順は、前半燃焼期間推定部３の説明において上述したとおりである。また、筒内燃料量は、熱発生率波形算出装置１の入力情報として与えられる。

このようにして、熱発生率波形の特性値の一つである熱発生率傾きｂ／ａを、基本的には熱発生率最大時燃料密度ρ_fuel@dQpeakの一次関数（この実施形態では、一例として比例関数）として算出することができる。言い換えると、負荷率および点火時期ＳＡを考慮せず、主に熱発生率最大時燃料密度ρ_fuel@dQpeakに基づいて熱発生率傾きｂ／ａを推定できることから、エンジンの種々の運転条件における熱発生率傾きｂ／ａの決定にかかる工数が低減される。

−熱発生量推定部−
熱発生量推定部５は、前述した如く、燃焼期間の全期間において発生した熱発生量（総熱発生量Ｑ_all）を推定する部分である。

以下、この熱発生量推定部における総熱発生量Ｑ_allの推定動作について説明する。

まず、前半燃焼期間ａにおける熱発生量Ｑ１は以下の式（５）によって算出される。

一方、総熱発生量Ｑ_allは、筒内燃料量×ｋ（燃焼効率）として表すことができる（この式が熱発生量推定モデルに相当する）。燃焼効率ｋは、例えば暖機運転中のように油水温が低い場合に低下すると共に、負荷率やエンジン回転速度などの変化によっても変化する。そこで、本実施形態では、種々のエンジンの実験結果のデータベースを利用して燃焼効率ｋの値を、油水温、負荷率およびエンジン回転速度に対応づけて決定するマップを予め決定しておく。そして、この燃焼効率ｋの値を用いて、筒内燃料量から総熱発生量Ｑ_allが算出される。

図２を参照して上述したように、熱発生率波形を作成するためには、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂおよび後半燃焼期間ｃを求める必要がある。この後半燃焼期間ｃにおける熱発生量Ｑ２は以下の式（６）によって算出される。

そして、熱発生率最大時期ｄＱｐｅａｋＡにおける熱発生率ｂは以下の式（７）により、また、後半燃焼期間ｃは以下の式（８）によりそれぞれ求められる。

以上のようにして熱発生率波形算出装置１では、着火遅れ推定モデルを使用した着火遅れ推定部２での着火遅れ期間τの推定、前半燃焼期間推定モデルを使用した前半燃焼期間推定部３での前半燃焼期間ａの推定、熱発生率傾き推定モデルを使用した熱発生率傾き推定部４での熱発生率傾きｂ／ａの推定、熱発生量推定モデルを使用した熱発生量推定部５での熱発生量Ｑ_allの推定が行われると共に、最大熱発生率ｂおよび後半燃焼期間ｃが算出される。そして、これらの算出値を用いて、図２のように熱発生率波形を近似する三角波形を作成し、この三角波形を熱発生率波形として出力する。この出力された熱発生率波形を利用して、エンジン設計に際してのシステムの検討、制御の検討、適合値の検討が行われることになる。

以上説明したように本実施形態では、エンジンの熱発生率波形に近似した三角波形を作成するに際し、その波形の特性値の一つである前半燃焼期間ａを用いている。そして、この前半燃焼期間ａは、エンジンの負荷率、ＥＧＲ率、空燃比および油水温の運転条件の影響を受けないことが新たな知見として見出されている。つまり、前半燃焼期間ａを、エンジンの負荷率、ＥＧＲ率、空燃比および油水温のいずれにも依らないものとして、熱発生率最大時筒内体積Ｖ_@dQpeakとエンジン回転速度Ｎｅ（より具体的にはタンブル比に応じた値βを指数とするエンジン回転速度Ｎｅの指数関数）とに基づいて算出している。そして、本実施形態にあっては、この前半燃焼期間ａに基づいて熱発生率波形を作成するようにしている。このため、前記各運転条件を考慮する必要がない。従って、負荷率、ＥＧＲ率、空燃比および油水温の運転条件を考慮して（これら運転条件を変化させて）熱発生率波形を作成するといった従来技術（例えば、ウィーベ関数を利用し、複数のパラメータ（形状パラメータなど）を種々の運転条件に対して同定するもの）に比べて、十分な精度を保証しながらも熱発生率波形を作成する工数の低減を図ることができる。

よって、熱発生率波形を利用して、エンジン設計に際しての種々の検討などを効率的に行うことによって、開発コストを削減できる。

また、熱発生率最大時筒内体積Ｖ_@dQpeakとエンジン回転速度Ｎｅとに基づいて前半燃焼期間ａを算出し、この前半燃焼期間ａに基づいて熱発生率波形を作成していることから、この熱発生率波形は、気筒内での燃焼状態における物理現象に従って作成されたものとなっている。この点で、単に波形形状を合わせ込むために、形状パラメータなどの各種パラメータを数学的に合わせ込むといったウィーベ関数を利用した熱発生率波形の作成手法に比べて、本実施形態に係る熱発生率波形算出装置１によって作成された熱発生率波形は高い信頼性を得ることができる。

さらに、本実施形態では、熱発生率波形全体を作成しなくても、前記のように熱発生率最大時筒内体積Ｖ_@dQpeakとエンジン回転速度Ｎｅとから前半燃焼期間ａを算出することができるので、この前半燃焼期間ａについて所要の精度を確保しながら、簡易に推定し、評価することができる。

また、前述した如く、本実施形態で推定される前半燃焼期間ａと熱発生率傾きｂ／ａとは、互いに独立した（従属関係にない）値として推定されるものとなっている。このため、前半燃焼期間ａは、熱発生率傾きｂ／ａに推定誤差が含まれる場合にその誤差の影響を受けることの無い値として推定され、熱発生率傾きｂ／ａは、前半燃焼期間ａに推定誤差が含まれる場合にその誤差の影響を受けることの無い値として推定される。その結果、これら値の推定精度を高く確保することができる。

−他の実施形態−
以上説明した実施形態は、自動車用のガソリンエンジンを対象とした熱発生率波形を作成する熱発生率波形算出装置に本発明を適用した場合について説明した。本発明はこれに限らず、自動車用以外の火花点火機関に対しても適用が可能である。また、ガソリンエンジンにも特に限定されるものではなく、例えばガスエンジンに対しても適用が可能である。

また、前記の実施形態で説明した熱発生率波形算出装置において実施される熱発生率波形算出方法も本発明の技術的思想の範疇である。

前記の実施形態では、混合気の着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの期間における熱発生率の平均的な増加率を熱発生率傾きｂ／ａとし、式（４）のように、この熱発生率傾きｂ／ａを熱発生率最大時燃料密度ρ_fuel@dQpeakの一次関数として算出するようにしているが、これにも限定されない。

すなわち、着火時期ＦＡから熱発生率最大時期ｄＱｐｅａｋＡまでの間で熱発生率が増加している間（熱発生率増加期間）で、例えば、着火時期から熱発生率最大時期ｄＱｐｅａｋＡの少し手前の所定時期までの期間における熱発生率の増加率を熱発生率傾きとし、この熱発生率傾きを、前記所定時の燃料密度に基づいて推定するようにしてもよい。

また、前記の実施形態では、熱発生率最大時筒内体積Ｖ_@dQpeakおよび前半燃焼期間ａを求める手法として、仮想の熱発生率最大時期を設定し、前記式（３）の繰り返し計算によって算出するようにしていた。本発明はこれに限らず、実機での試験において熱発生率最大時期をセンシングし、これに基づいて熱発生率最大時を設定するようにしたり、熱発生率波形算出装置１への入力信号として所望の熱発生率最大時期を入力したりすることによって、熱発生率最大時筒内体積Ｖ_@dQpeakおよび前半燃焼期間ａを求めるようにしてもよい。

さらに、前記の実施形態における混合気の着火遅れ期間τの算出の仕方も一例に過ぎず、これらに限定されることはない。例えば、着火遅れ期間τを、実験、シミュレーションによって決定するようにしてもよい。

また、前記の実施形態に係る熱発生率波形算出装置１は、三角波形を出力するものであった。本発明はこれに限らず、作成した三角波形に対して所定のフィルタ処理を行って熱発生率波形を作成し、この熱発生率波形を出力するものであってもよい。

また、前記の実施形態では、熱発生率最大時期における物理量として筒内体積を挙げたが、気筒内の乱れを規定する（気筒内の乱れに影響を与える）他の物理量であってもよい。

また、前記の実施形態では、前半燃焼期間ａを、エンジンの負荷率、ＥＧＲ率、空燃比および油水温のいずれにも依らないものとして算出していたが、これら運転条件のうち少なくとも一つに依らないものとして算出するようにしてもよい。

本発明によると、火花点火式の内燃機関における熱発生率波形の作成にかかる工数を低減でき、コストの削減が可能になるので、例えば自動車用の内燃機関に適用可能である。

１ 熱発生率波形算出装置
ＳＡ 点火時期
ＦＡ 混合気の着火時期
ａ 前半燃焼期間（着火時期から熱発生率最大時期までの期間）
ｂ 熱発生率最大時期における熱発生率
ｄＱｐｅａｋＡ 熱発生率最大時期
Ｖ_@dQpeak 熱発生率最大時筒内体積
Ｎｅ^β エンジン回転速度に基づく補正係数
Ｖ_@dQpeak ^α 熱発生率最大時筒内体積に基づく補正係数

Claims (4)

1. 火花点火式の内燃機関における熱発生率波形を算出するための装置であって、
   混合気の燃焼期間のうち着火時期から熱発生率が最大になる熱発生率最大時期までの期間を、前記熱発生率波形の特性値の一つである前半燃焼期間として規定し、
   前記前半燃焼期間を、前記熱発生率最大時期における物理量に基づいて推定し、この推定した前半燃焼期間を用いて前記熱発生率波形を算出する構成となっており、且つ、
   前記前半燃焼期間を推定するための前記熱発生率最大時期における物理量は、熱発生率最大時期における筒内体積および機関回転速度であって、下記の式（３）により、これら熱発生率最大時期における筒内体積および機関回転速度によって前記前半燃焼期間が推定される構成となっていることを特徴とする内燃機関の熱発生率波形算出装置。
   

   

   （ａは前半燃焼期間、Ｖ _@dQpeak は熱発生率最大時期における筒内体積、Ｎｅは機関回転速度、Ｃおよびαは実験に基づいた同定によって得られる係数、βは気筒内でのタンブル比に応じタンブル比が大きいほど大きな値として与えられる係数）
2. 請求項１に記載の内燃機関の熱発生率波形算出装置において、
   前記前半燃焼期間は、
   仮想の熱発生率最大時期を設定し、実際の点火時期に従って求められた着火時期と前記仮想の熱発生率最大時期との間の期間である仮想の前半燃焼期間と、前記仮想の熱発生率最大時期における物理量に基づいて推定された前半燃焼期間とを比較し、これらが一致した場合における前記推定された前半燃焼期間が真の前半燃焼期間として算出され、
   この真の前半燃焼期間を用いて前記熱発生率波形を算出する構成となっている、内燃機関の熱発生率波形算出装置。
3. 請求項１または２に記載の内燃機関の熱発生率波形算出装置において、
   混合気の着火から燃焼終了までのクランク角度期間を底辺とし、熱発生率最大時期における熱発生率を頂点とする三角波形によって熱発生率波形を近似し、
   前記三角波形において着火時期から熱発生率最大時期までの熱発生率を表す三角形の底辺を前記前半燃焼期間として規定する構成となっている、内燃機関の熱発生率波形算出装置。
4. 火花点火式の内燃機関における熱発生率波形を算出する方法であって、
   混合気の燃焼期間のうち着火時期から熱発生率が最大になる熱発生率最大時期までの期間を、前記熱発生率波形の特性値の一つである前半燃焼期間として規定し、
   前記前半燃焼期間を、前記熱発生率最大時期における物理量に基づいて推定し、この推定した前半燃焼期間を用いて前記熱発生率波形を算出し、且つ、
   前記前半燃焼期間を推定するための前記熱発生率最大時期における物理量は、熱発生率最大時期における筒内体積および機関回転速度であって、下記の式（３）により、これら熱発生率最大時期における筒内体積および機関回転速度によって前記前半燃焼期間を推定することを特徴とする内燃機関の熱発生率波形算出方法。
   

   

   （ａは前半燃焼期間、Ｖ _@dQpeak は熱発生率最大時期における筒内体積、Ｎｅは機関回転速度、Ｃおよびαは実験に基づいた同定によって得られる係数、βは気筒内でのタンブル比に応じタンブル比が大きいほど大きな値として与えられる係数）

Images