(57)【要約】      （修正有）【課題】吸気弁の開弁特性が様々に変更され得る内燃機関の前半燃焼期間長を効率良く推定する。【解決手段】火花点火式内燃機関の機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁の開弁特性が基準特性である場合における「熱発生率最大時期θdQpeakと前半燃焼期間長ａrefとの関係（即ち、基準関係）」がａref＝ｆ(θdQpeak）なる関数で表される場合、この状態から吸気弁の開弁特性が特定特性へと変化したとき、任意の熱発生率最大時期θ２に対する前半燃焼期間長ａを「ａ＝ｆ(θ２）＋Δａ１」なる式に基づいて推定する。ここで、Δａ１は、吸気弁の開弁特性が特定特性である場合であって熱発生率最大時期θdQpeakが基準時期θ１であるときの前半燃焼期間長ad1と、吸気弁の開弁特性が基準特性である場合であって熱発生率最大時期θdQpeakが基準時期θ１であるときの前半燃焼期間長ar1と、の差分値である。【選択図】図１
2016160761.tif 000002
【特許請求の範囲】
【請求項１】
火花点火式内燃機関の回転速度が所定の基準回転速度であり且つ前記機関の吸気弁の開弁特性が所定の基準開弁特性に設定された基準状態における、前記機関の燃焼室内での燃料の着火時期から前記燃料の燃焼に伴う熱の発生率が最大となる熱発生率最大時期までの期間である前半燃焼期間の長さを表す前半燃焼期間長と、前記熱発生率最大時期と、の関係を基準関係として推定又は実測により取得し、前記機関の回転速度が前記基準回転速度であり且つ前記吸気弁の開弁特性が前記基準開弁特性と異なる特定開弁特性に設定された特定状態における前記熱発生率最大時期と前記前半燃焼期間長との関係を前記基準関係に基づいて推定する燃焼状態推定方法であって、
前記基準状態において前記熱発生率最大時期が所定の第１時期であるときの前記前半燃焼期間長を前記基準関係に基づいて又は実測により第１基準期間長として取得し、
前記特定状態において前記熱発生率最大時期が前記第１時期であるときの前記前半燃焼期間長を第１特定期間長として実測により取得し、
前記第１特定期間長から前記第１基準期間長を減じることにより第１差分値を取得するか、又は、前記第１基準期間長から前記第１特定期間長を減じることにより第２差分値を取得し、
前記基準状態において前記熱発生率最大時期が前記第１時期と異なる第２時期であるときの前記前半燃焼期間長を前記基準関係に基づいて第２基準期間長として取得し、
前記特定状態において前記熱発生率最大時期が前記第２時期であるときの前記前半燃焼期間長を、前記第２基準期間長に前記第１差分値を加えるか又は前記第２基準期間長から前記第２差分値を減じることにより推定する、
燃焼状態推定方法。
【請求項２】
請求項１に記載の燃焼状態推定方法において、
前記特定開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）と、前記吸気弁が開弁している期間をクランク角幅により表した吸気弁作用角（ＶＣＡＭ）と、のうちの少なくとも一方が、前記基準開弁特性と異なる、燃焼状態推定方法。
【請求項３】
請求項２に記載の燃焼状態推定方法において、
前記基準開弁特性は、前記吸気弁作用角（ＶＣＡＭ）が所定の基準作用角であり、且つ、前記吸気弁の開弁時期と前記吸気弁の閉弁時期との中心を示すクランク角と所定の基準クランク角とのクランク角差である吸気弁位相角（ＩＮＶＴ）が所定の基準位相角である特性であり、
前記特定開弁特性は、前記吸気弁作用角（ＶＣＡＭ）が前記基準作用角であり、且つ、前記吸気弁位相角（ＩＮＶＴ）が前記基準位相角と異なる特定位相角であることによって前記吸気弁の閉弁時期（ＩＶＣ）が前記基準開弁特性における前記吸気弁の閉弁時期（ＩＶＣ）と異なる特性である、
燃焼状態推定方法。
【請求項４】
請求項２に記載の燃焼状態推定方法において、
前記基準開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）が所定の基準閉弁時期であり、且つ、前記吸気弁作用角（ＶＣＡＭ）が所定の基準作用角である特性であり、
前記特定開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）が前記基準閉弁時期であり、且つ、前記吸気弁作用角（ＶＣＡＭ）が前記基準作用角と異なる特定作用角である特性である、
燃焼状態推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、内燃機関における燃料の燃焼に係る前半燃焼期間の長さ（前半燃焼期間長）を、燃焼状態を表すパラメータの一つとして推定する燃焼状態推定方法に関する。
【背景技術】
【０００２】
内燃機関における燃料の燃焼状態を推定するため、燃焼室内における燃料の燃焼により発生する熱の量（熱発生量、発熱量）Ｑの推移を推定する種々の手法が提案されている。特に、あるクランク角における熱発生率（ｄＱ／ｄθ）の推移（即ち、熱発生率波形）は内燃機関の燃焼状態を推定する上で重要である。なお、熱発生率ｄＱ／ｄθは、クランク軸の回転に対する単位クランク角あたりの熱発生量［Ｊ／ＣＡ］と定義される。
【０００３】
Ｗｉｅｂｅ関数は、種々の変数と熱発生率との関係を規定した関数であり、熱発生率を推定するための有効な関数の一つとして良く知られている。Ｗｉｅｂｅ関数において使用される種々の変数は、例えば、形状パラメータ、効率及び着火遅れ期間等である。特許文献１は、これらの変数を「機関運転状態を表すパラメータ（例えば、点火時期、機関回転速度、負荷率及び空燃比等）」に基づいて決定する手法を開示している。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００７−１７７６５４号公報
【発明の概要】
【０００５】
上記特許文献１は、更に、Ｗｉｅｂｅ関数に基づいて前半燃焼期間を推定可能であることも開示している。前半燃焼期間とは、「火花点火式内燃機関における燃料の着火時期」から「その燃料の燃焼に伴う熱の発生率（熱発生率ｄＱ／ｄθ）が最大となる時期（即ち、熱発生率最大時期）」までの期間である。前半燃焼期間の長さ（クランク角幅）は「前半燃焼期間長」とも称呼される。前半燃焼期間長は内燃機関の燃焼状態を推定する上で重要なパラメータの一つである。
【０００６】
しかしながら、前半燃焼期間長は、吸気弁の開弁特性（例えば、後述する「吸気弁位相角、吸気弁作用角、吸気弁開弁時期及び吸気弁閉弁時期」等により決まる、吸気弁が開弁している期間に係る特性値）に応じて大きく変化する。これは、吸気弁の開弁特性が変更されると筒内（燃焼室内）の気流の乱れの強さが変化し、それにより、燃焼速度が変化することに起因すると推察される。
【０００７】
そのため、仮に、Ｗｉｅｂｅ関数を利用して前半燃焼期間長を推定したとしても、吸気弁の開弁特性が変化した場合には、再度、Ｗｉｅｂｅ関数の各変数を種々の機関運転状態パラメータに対して決定し直す必要がある。更に、実測により前半燃焼期間長を取得する場合、吸気弁の開弁特性が変更される毎に各種の機関運転状態に対して膨大な実験を行わなければならない。加えて、本明細書において開示する「本出願人が提案しているモデル（モデル式）」を用いて前半燃焼期間長を推定したとしても、同様に、吸気弁の開弁特性が変化する毎に各種のパラメータを同定し直す必要がある。即ち、前半燃焼期間長を吸気弁の種々の開弁特性に対して精度良く取得するには、極めて多くの労力及び時間を費やす必要があった。
【０００８】
本発明は上述した課題に対処するためになされた。即ち、本発明の目的の一つは、吸気弁の開弁特性が様々に変更され得る内燃機関の前半燃焼期間長を効率良く推定することができる燃焼状態推定方法を提供することにある。
【０００９】
発明者は、吸気弁の開弁特性が変化した場合に前半燃焼期間長がどのように変化するかにつき詳細な検討を行った。以下、その検討結果について説明する。
【００１０】
図１は、機関回転速度が所定速度（一定値）であり、且つ、吸気弁作用角ＶＣＡＭが所定角（一定値）である第１の場合に、種々の吸気弁位相角ＩＮＶＴに対して「熱発生率最大時期θdQpeak及び前半燃焼期間長ａ」を実測した結果の概略を示すグラフである。
【００１１】
なお、熱発生率最大時期θdQpeakはクランク角θ［deg］により表され、前半燃焼期間長ａはクランク角幅［ＣＡ］により表される。単位、［deg］及び［ＣＡ］は、何れも、クランク軸が１回転する角度を３６０°としたときのクランク角１°を示す。一般に、［deg］は特定のクランク角を指し示す場合に使用され、［ＣＡ］はクランク角の幅を示す場合に使用される。
【００１２】
ここで、吸気弁作用角ＶＣＡＭ及び吸気弁位相角ＩＮＶＴについて、図２の（Ａ）及び（Ｂ）を参照して説明する。図２の（Ａ）はクランク角θと吸気弁リフト量との関係を示し、図２の（Ｂ）は「吸気上死点から圧縮上死点までを円により示すとともに、その円に吸気弁が開弁している期間（即ち、吸気弁開弁時期ＩＶＯから吸気弁閉弁時期ＩＶＣまでの吸気弁開弁期間）を重ねて描いた図である。図２の（Ａ）及び（Ｂ）において、実線Ｃｒにより示された開弁特性（以下、便宜上「Ａ特性」と称呼する。）と破線Ｃｘにより示された開弁特性（以下、便宜上「Ｂ特性」と称呼する。）とは、吸気弁作用角ＶＣＡＭが互いに同じであり、吸気弁位相角ＩＮＶＴが互いに相違している。なお、吸気弁閉弁時期ＩＶＣは、吸気下死点後のクランク角［ABDC deg］により表す。
【００１３】
図２の（Ａ）及び（Ｂ）から理解されるように、吸気弁作用角ＶＣＡＭは「吸気弁開弁時期ＩＶＯから吸気弁閉弁時期ＩＶＣまでの期間に対応するクランク角幅」である。図２の（Ａ）及び（Ｂ）において、Ａ特性の吸気弁作用角ＶＣＡＭ及びＢ特性の吸気弁作用角ＶＣＡＭは、何れも所定値（＝ＶＣ１）である。
【００１４】
吸気弁位相角ＩＮＶＴは、吸気弁開弁時期ＩＶＯと吸気弁閉弁時期ＩＶＣとの中心が位置するクランク角（即ち、吸気弁開弁期間の中心を示すクランク角）が「基準クランク角」に対してどれだけ進角させられているかを示すクランク角幅である。ここで、基準クランク角は、「吸気弁開弁期間が吸気弁制御装置の作動範囲内で最も遅角側に設定された場合における、吸気弁開弁時期ＩＶＯと吸気弁閉弁時期ＩＶＣの中心が位置するクランク角」である。従って、Ａ特性における吸気弁位相角ＩＮＶＴが「０」であるとし（即ち、Ａ特性が、吸気弁開弁期間の中心が最も遅角側の値（＝ｖｔｒ）に設定されている特性であるとし）、Ｂ特性の吸気弁開弁期間の中心が値ｖｔｘであるとすると、Ｂ特性の吸気弁位相角ＩＮＶＴは、｜ｖｔｒ−ｖｔｘ｜となる。
【００１５】
従って、図２の（Ａ）及び（Ｂ）に示したように、吸気弁作用角ＶＣＡＭが一定値ＶＣ１である場合、吸気弁位相角ＩＮＶＴがクランク角θｓであるとき、吸気弁開弁時期ＩＶＯは最も遅角側の吸気弁開弁時期ＩＶＯｒから吸気弁開弁時期ＩＶＯｘへとクランク角θｓだけ進角し、吸気弁閉弁時期ＩＶＣは最も遅角側の吸気弁閉弁時期ＩＶＣｒから吸気弁閉弁時期ＩＶＣｘへとクランク角θｓだけ進角する。
【００１６】
再び、図１を参照すると、
実線Ｃｒは吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒである場合における熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示している。ここでは、基準位相角ＩＮＶＴｒは「０」より大きい所定値に設定されている。
一点鎖線Ｃｓは吸気弁位相角ＩＮＶＴが「基準位相角ＩＮＶＴｒと異なる位相角（例えば、基準位相角ＩＮＶＴｒよりも小さい遅角側の値）ＩＮＶＴｓ」である場合における熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示している。
二点鎖線Ｃｄは吸気弁位相角ＩＮＶＴが「基準位相角ＩＮＶＴｒ及び位相角ＩＮＶＴｓの何れとも異なる位相角（例えば、基準位相角ＩＮＶＴｒよりも大きい進角側の値）ＩＮＶＴｄ」である場合における熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示している。
【００１７】
このグラフから理解されるように、発明者は、機関回転速度が所定速度に維持され且つ吸気弁作用角ＶＣＡＭが所定角に維持されている場合、吸気弁位相角ＩＮＶＴが第１位相角（ＩＮＶＴｒ）から第２位相角（ＩＮＶＴｓ又はＩＮＶＴｄ）へと変化すると、熱発生率最大時期θdQpeakに関わらず前半燃焼期間長ａが「所定値（例えば、Δａ１又はΔａ２）」だけ変化することを見出した。換言すると、発明者は、機関回転速度が所定速度に維持され且つ吸気弁作用角ＶＣＡＭが所定角に維持されている場合、吸気弁位相角ＩＮＶＴが変化すると、熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示す曲線が前半燃焼期間長ａの方向（増大及び減少方向の何れか、即ち、図１においてＹ軸方向）に平行移動することを見出した。
【００１８】
より具体的に述べると、吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒである場合、熱発生率最大時期θdQpeakが基準値θ１のときの前半燃焼期間長ａがａr1であり（点Ｐｒ１を参照。）、熱発生率最大時期θdQpeakが任意値θ２のときの前半燃焼期間長ａがａr2であった（点Ｐｒ２を参照。）。更に、吸気弁位相角ＩＮＶＴが「基準位相角ＩＮＶＴｒと異なる所定位相角（ＩＮＶＴｒより小さい遅角側のＩＮＶＴｓ又はＩＮＶＴｒより大きい進角側のＩＮＶＴｄ）」である場合、熱発生率最大時期θdQpeakが基準値θ１のときの前半燃焼期間長ａがａd1であり（点Ｐａ１を参照。）、熱発生率最大時期θdQpeakが任意値θ２のときの前半燃焼期間長ａがａd2であった（点Ｐａ２を参照。）。そして、発明者は下記の（１Ａ）式が任意値θ２の総てについて成立することを見出した。更に、発明者は、後述するように、この（１Ａ）式が、他の機関回転速度ＮＥ及び他の吸気弁作用角ＶＣＡＭに対しても、それらが変化しない限り同様に成立することを見出した。
（ａd1−ａr1）＝（ａd2−ａr2）＝Δａ１      …（１Ａ）
【００１９】
従って、発明者は次の知見１を得た。
機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭ０であり且つ吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒである場合における「熱発生率最大時期θdQpeakと前半燃焼期間長ａrefとの関係（即ち、基準関係）」がａref＝ｆ(θdQpeak）なる関数で表されるとする。
このとき、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭ０であり、且つ、
吸気弁位相角ＩＮＶＴが所定位相角ＩＮＶＴｓであり且つ熱発生率最大時期θdQpeakが任意値（任意の基準値）θ１であるときの前半燃焼期間長ａｓ（θ１）がａd1であるとすると、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭ０であり、且つ、
吸気弁位相角ＩＮＶＴが所定位相角ＩＮＶＴｓであり且つ熱発生率最大時期θdQpeakが任意値θ２であるときの前半燃焼期間長ａｓ（θ２）は、次の（１Ｂ）式により求めることができる。
ａｓ（θ２）＝ｆ(θ２）＋Δａ１
＝ｆ(θ２）＋（ａd1−ａr1）
＝ｆ(θ２）−（ａr1−ａd1）
＝ｆ(θ２）＋{ａd1−ｆ（θ１）}    …（１Ｂ）
なお、値（ａd1−ａr1）は便宜上、「第１差分値」とも称呼される。値（ａr1−ａd1）は便宜上、「第２差分値」とも称呼される。従って、上記（１Ｂ）式は、ｆ（θ２）に第１差分値を加えることにより前半燃焼期間長ａｓ（θ２）を推定する式であり、ｆ（θ２）から第２差分値を減じることにより前半燃焼期間長ａｓ（θ２）を推定する式であると表現することができる。
【００２０】
それ故、上記基準関係（ａref＝ｆ(θdQpeak））が実測により又はモデル等を用いた推定に基づいて得られていて、更に、上記前半燃焼期間長ａd1が実測等により得られていれば、吸気弁位相角ＩＮＶＴが位相角ＩＮＶＴｓである場合の「熱発生率最大時期θdQpeakと前半燃焼期間長ａｓとの関係」は上記（１Ｂ）式から推定することができる。
【００２１】
一方、図３は、機関回転速度が所定速度（一定値）であり、且つ、吸気弁の閉弁時期ＩＶＣが所定クランク角（一定値）である場合に、種々の吸気弁作用角ＶＣＡＭに対して「熱発生率最大時期θdQpeak及び前半燃焼期間長ａ」を実測した結果の概略を示すグラフである。
【００２２】
ここで、図３に示した結果の前提条件の理解を容易にするため、図４の（Ａ）に、吸気弁閉弁時期ＩＶＣを一定のクランク角ＩＶＣｒに維持し、吸気弁作用角ＶＣＡＭを作用角ＶＣｒに設定した場合（実線Ｃｒを参照。）及び吸気弁作用角ＶＣＡＭを「作用角ＶＣｒよりも小さな作用角ＶＣｓ」に設定した場合（破線Ｃｘを参照。）における「クランク角θと吸気弁リフト量との関係」を示す。更に、図２の（Ｂ）と同様の図４の（Ｂ）に、図４の（Ａ）の実線Ｃｒ及び破線Ｃｘに対応した吸気弁の開弁期間を実線Ｃｒ及び破線Ｃｘによりそれぞれ示す。この図４の（Ａ）及び（Ｂ）から理解されるように、吸気弁の閉弁時期ＩＶＣが所定クランク角（一定値）である場合に、吸気弁作用角ＶＣＡＭが変化すると、少なくとも吸気弁の開弁時期ＩＶＯが変化し、実際には、吸気弁リフト量及び吸気弁位相角ＩＮＶＴも変化する。
【００２３】
再び、図３を参照すると、
実線Ｌｒは吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒである場合における熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示している。
一点鎖線Ｌｓは吸気弁作用角ＶＣＡＭが「基準作用角ＶＣＡＭｒと異なる作用角（例えば、基準作用角ＶＣＡＭｒよりも小さい作用角）ＶＣＡＭｓ」である場合における熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示している。
二点鎖線Ｌｄは吸気弁作用角ＶＣＡＭが「基準作用角ＶＣＡＭｒ及び作用角ＶＣＡＭｓの何れとも異なる作用角（例えば、基準作用角ＶＣＡＭｒよりも大きい作用角）ＶＣＡＭｄ」である場合における熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示している。
【００２４】
このグラフから理解されるように、発明者は、機関回転速度が所定速度に維持され且つ吸気弁閉弁時期ＩＶＣが所定クランク角に維持されている場合、吸気弁作用角ＶＣＡＭが第１作用角（ＶＣＡＭｒ）から第２作用角（ＶＣＡＭｓ又はＶＣＡＭｄ）へと変化すると、熱発生率最大時期θdQpeakに関わらず前半燃焼期間長ａが「所定値（例えば、Δａ３又はΔａ４）」だけ変化することを見出した。換言すると、発明者は、機関回転速度が所定速度に維持され且つ吸気弁閉弁時期ＩＶＣが所定クランク角に維持されている場合、吸気弁作用角ＶＣＡＭが変化すると、熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を示す曲線が前半燃焼期間長ａの方向（増大及び減少方向の何れか、即ち、図３においてＹ軸方向）に平行移動することを見出した。
【００２５】
より具体的に述べると、吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒである場合、熱発生率最大時期θdQpeakが基準値θ３のときの前半燃焼期間長ａがａr3であり（点Ｐｒ３を参照。）、熱発生率最大時期θdQpeakが任意値θ4のときの前半燃焼期間長ａがａr4であった（点Ｐｒ４を参照。）。更に、吸気弁作用角ＶＣＡＭが「基準作用角ＶＣＡＭｒと異なる所定作用角（ＶＣＡＭｒよりも大きい作用角ＶＣＡＭｄ又は基準作用角ＶＣＡＭｒよりも小さい作用角ＶＣＡＭｓ）」である場合、熱発生率最大時期θdQpeakが基準値θ３のときの前半燃焼期間長ａがａd3であり（点Ｐａ３を参照。）、熱発生率最大時期θdQpeakが任意値θ４のときの前半燃焼期間長ａがａd4であった（点Ｐａ４を参照。）。そして、発明者は下記の（２Ａ）式が任意値θ４の総てについて成立することを見出した。更に、発明者は、後述するように、この（２Ａ）式は、他の機関回転速度ＮＥ及び他の吸気弁閉弁時期ＩＶＣに対しても、それらが変化しない限り同様に成立することを見出した。
（ａd3−ａr3）＝（ａd4−ａr4）＝Δａ３      …（２Ａ）
【００２６】
従って、発明者は次の知見２を得た。
機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁閉弁時期ＩＶＣが任意閉弁時期ＩＶＣ０であり且つ吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒである場合における「熱発生率最大時期θdQpeakと前半燃焼期間長ａrefとの関係（即ち、基準関係）」がａref＝ｇ(θdQpeak）なる関数で表されるとする。
このとき、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁閉弁時期ＩＶＣが任意閉弁時期ＩＶＣ０であり、且つ、
吸気弁作用角ＶＣＡＭが所定作用角ＶＣＡＭｓであり且つ熱発生率最大時期θdQpeakが任意値（任意の基準値）θ３であるときの前半燃焼期間長ａｓ（θ３）がａd3であるとすると、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁閉弁時期ＩＶＣが任意閉弁時期ＩＶＣ０であり、且つ、
吸気弁作用角ＶＣＡＭが所定作用角ＶＣＡＭｓであり且つ熱発生率最大時期θdQpeakが任意値θ４であるときの前半燃焼期間長ａｓ（θ４）は、次の（２Ｂ）式により求めることができる。
ａｓ(θ４）＝ｇ(θ４）＋Δａ３
＝ｇ(θ４）＋（ａd3−ａr3）
＝ｇ(θ４）−（ａr3−ａd3）
＝ｇ(θ４）＋｛ａd3−ｇ（θ３）｝    …（２Ｂ）
なお、値（ａd3−ａr3）は便宜上、「第１差分値」とも称呼される。値（ａr3−ａd3）は便宜上、「第２差分値」とも称呼される。従って、上記（２Ｂ）式は、ｇ(θ４）に第１差分値を加えることにより前半燃焼期間長ａｓ（θ４）を推定する式であり、ｇ（θ４）から第２差分値を減じることにより前半燃焼期間長ａｓ（θ４）を推定する式であると表現することができる。
【００２７】
それ故、上記基準関係（ａref＝g(θdQpeak））が実測により又はモデル等を用いた推定に基づいて得られていて、更に、上記前半燃焼期間長ａd3が実測等により得られていれば、吸気弁作用角ＶＣＡＭが作用角ＶＣＡＭｓである場合の「熱発生率最大時期θdQpeakと前半燃焼期間長ａｓとの関係」は上記（２Ｂ）式から推定することができる。
【００２８】
ところで、図２及び図４を参照しながら行った説明及び下記の関係式（＊１）及び（＊２）から理解されるように、吸気弁閉弁時期（ＩＶＣ）、吸気弁開弁時期（ＩＶＯ）、吸気弁作用角（ＶＣＡＭ）及び吸気弁位相角（ＩＮＶＴ）の間には、これらのうちの二つが決まれば残りの二つは自動的に決まるという関係がある。なお、ここでは、吸気弁開弁時期ＩＶＯも、吸気下死点後のクランク角［ABDC deg］により表している。
（ＩＶＣ−ＩＶＯ）＝ＶＣＡＭ                  …（＊１）
（ＩＶＣ＋ＩＶＯ）／２＝Ｖｔｒ−ＩＮＶＴ      …（＊２）
【００２９】
従って、例えば、ある吸気弁の開弁特性における「吸気弁の閉弁時期（ＩＶＣ）及び吸気弁作用角（ＶＣＡＭ）」が、他の吸気弁の開弁特性における「吸気弁の閉弁時期（ＩＶＣ）及び吸気弁作用角（ＶＣＡＭ）」のそれぞれと同じであれば、それら両方の吸気弁の開弁特性は互いに同じであると言うことができる。換言すれば、前記特定開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）と、前記吸気弁作用角（ＶＣＡＭ）と、のうちの少なくとも一方が、前記基準開弁特性と異なる特性であり得る。或いは、前記特定開弁特性は、前記吸気弁位相角（ＩＮＶＴ）と、前記吸気弁作用角（ＶＣＡＭ）と、のうちの少なくとも一方が、前記基準開弁特性と異なる特性であり得る。
【００３０】
更に、吸気弁カムシャフト及びカムにより吸気弁を開閉させる通常の内燃機関においては、吸気弁のリフト量（ＩＶＬift）の最大値（最大リフト量）は吸気弁作用角（ＶＣＡＭ）が決まればある特定の値になる。よって、前記特定開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）と、前記吸気弁の最大リフト量と、のうちの少なくとも一方が、前記基準開弁特性と異なる特性であり得る。或いは、前記特定開弁特性は、前記吸気弁位相角（ＩＮＶＴ）と、前記吸気弁の最大リフト量と、のうちの少なくとも一方が、前記基準開弁特性と異なる特性であり得る。
【００３１】
ここで、上記知見１は、下記の場合１と下記の場合２との間で「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を表す曲線」が前半燃焼期間長ａが増大又は減少する方向に平行移動することを意味している。
場合１
吸気弁位相角ＩＮＶＴ＝Ａ１
吸気弁作用角ＶＣＡＭ＝Ｂ１
吸気弁閉弁時期ＩＶＣ＝Ｃ１
場合２
吸気弁位相角ＩＮＶＴ＝Ａ２
吸気弁作用角ＶＣＡＭ＝Ｂ１
吸気弁閉弁時期ＩＶＣ＝Ｃ２
【００３２】
上記知見２は、吸気弁作用角ＶＣＡＭを適切な値に設定すれば、上記場合２と下記の場合３との間で「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を表す曲線」が前半燃焼期間長ａが増大又は減少する方向に平行移動することを意味している。
場合３
吸気弁位相角ＩＮＶＴ＝Ａ１
吸気弁作用角ＶＣＡＭ＝Ｂ２
吸気弁閉弁時期ＩＶＣ＝Ｃ２
【００３３】
従って、上記場合１と上記場合３との間でも（即ち、吸気弁位相角ＩＮＶＴが一定値Ａ１に維持されていて、吸気弁作用角ＶＣＡＭが変化した場合（その結果として、吸気弁閉弁時期ＩＶＣが変化した場合）にも、「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を表す曲線」が前半燃焼期間長ａが増大又は減少する方向に平行移動することを意味している。
【００３４】
以上から、発明者は、機関の回転速度が一定速度である場合、「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を表す曲線」は、前半燃焼期間長ａが増大又は減少する方向に「基準開弁特性及び特定開弁特性によって決まる一定量」だけ平行移動するという知見３を得た。
【００３５】
係る知見１乃至３に基づいてなされた本発明の燃焼状態推定方法は、
火花点火式内燃機関の回転速度が所定の基準回転速度であり且つ前記機関の吸気弁の開弁特性が所定の基準開弁特性に設定された基準状態における「熱発生率最大時期と前半燃焼期間長との関係」を、例えば、Ｗｅｉｂｅ関数及び後述する本明細書に開示されるモデルに基づく推定により、或いは、実測により、基準関係として取得し、
前記機関の回転速度が前記基準回転速度であり且つ前記吸気弁の開弁特性が「前記基準開弁特性と異なる特定開弁特性」に設定された特定状態における「熱発生率最大時期と前半燃焼期間長との関係」を前記基準関係に基づいて推定する燃焼状態推定方法であり、以下に述べるステップを含む。但し、以下に述べるステップは、矛盾が生じない範囲において実行順序を問わない。
【００３６】
前記基準状態において前記熱発生率最大時期が所定の第１時期であるときの前記前半燃焼期間長を前記基準関係に基づいて又は実測により第１基準期間長として取得するステップ（第１ステップ）；
前記特定状態において前記熱発生率最大時期が前記第１時期であるときの前記前半燃焼期間長を第１特定期間長として実測により取得するステップ（第２ステップ）；
前記第１特定期間長から前記第１基準期間長を減じることにより第１差分値（例えば、Δａ１＝ａd1−ａr1、又は、Δａ３＝ａd3−ａr3）を取得するか、又は、前記第１基準期間長から前記第１特定期間長を減じることにより第２差分値（例えば、−Δａ１＝ａr1−ａd1、又は、−Δａ３＝ａr3−ａd3）を取得するステップ（第３ステップ）；
前記基準状態において前記熱発生率最大時期が前記第１時期と異なる第２時期であるときの前記前半燃焼期間長を前記基準関係に基づいて第２基準期間長として取得するステップ（第４ステップ）；及び、
前記特定状態において前記熱発生率最大時期が前記第２時期であるときの前記前半燃焼期間長を、前記第２基準期間長に前記第１差分値を加えるか（例えば、ｆ(θ２）＋（ａd1−ａr1））、又は、前記第２基準期間長から前記第２差分値を減じる（例えば、ｆ(θ２）−（ａr1−ａd1））ことにより推定するステップ（第５ステップ）。
【００３７】
この方法によれば、機関の回転速度が基準回転速度である場合、「基準関係」及び「吸気弁の開弁特性が特定開弁特性となった場合の第１特定期間長」を取得することができれば、吸気弁の開弁特性が特定開弁特性へと変化した場合の熱発生率最大時期と前半燃焼期間長との関係を直ちに推定することができる。従って、種々の吸気弁の開弁特性における「熱発生率最大時期に対する前半燃焼期間長」を極めて少ない時間及び労力により推定することができる。
【００３８】
前記基準開弁特性及び前記特定開弁特性について、より具体的な例を挙げると、図１及び図２から理解されるように、
前記基準開弁特性は、前記吸気弁作用角（ＶＣＡＭ）が所定の基準作用角であり、且つ、前記吸気弁の開弁時期と前記吸気弁の閉弁時期との中心を示すクランク角と所定の基準クランク角とのクランク角差である吸気弁位相角（ＩＮＶＴ）が所定の基準位相角である特性であり、
前記特定開弁特性は、前記吸気弁作用角（ＶＣＡＭ）が前記基準作用角であり、且つ、前記吸気弁位相角（ＩＮＶＴ）が前記基準位相角と異なる特定位相角である特性であり得る。
【００３９】
同様に、図３及び図４から理解されるように、
前記基準開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）が所定の基準閉弁時期であり、且つ、前記吸気弁作用角（ＶＣＡＭ）が所定の基準作用角である特性であり、
前記特定開弁特性は、前記吸気弁の閉弁時期（ＩＶＣ）が前記基準閉弁時期であり、且つ、前記吸気弁作用角（ＶＣＡＭ）が前記基準作用角と異なる特定作用角である特性であり得る。
【００４０】
本発明の他の目的、他の特徴及び付随する利点は、以下の図面を参照しつつ記述される本発明の各実施形態についての説明から容易に理解されるであろう。尚、上記説明においては、発明の理解を助けるために、基準関係を表す関数、第１差分値及び第２差分値等に具体的な記号・符号を用いたが、これは何ら本発明を限定するものではない。
【図面の簡単な説明】
【００４１】
【図１】図１は、種々の吸気弁位相角に対する「熱発生率最大時期と前半燃焼期間長との関係」を示したグラフである。
【図２】図２は、吸気弁の開弁特性を規定するパラメータを説明するための図である。
【図３】図３は、種々の吸気弁作用角に対する「熱発生率最大時期と前半燃焼期間長との関係」を示したグラフである。
【図４】図４は、吸気弁の開弁特性を規定するパラメータを説明するための図である。
【図５】図５は、本発明の実施形態に係る装置（シミュレータ）の構成図である。
【図６】図６は、図５に示した装置（本実施装置）が実行する工程（即ち、本実施方法）を示したフローチャートである。
【図７】図７は、本実施装置が推定する熱発生率波形の一例を示した図である。
【図８】図８は、本実施装置のＣＰＵが実行するルーチンを示したフローチャートである。
【図９】図９は、ＢＴＤＣ着火の場合における、点火時期での筒内燃料密度と実測された着火遅れ期間長との関係を示したグラフである。
【図１０】図１０は、ＢＴＤＣ着火の場合において本実施装置により予測された着火遅れ期間長と実測された着火遅れ期間長との関係を示したグラフである。
【図１１】図１１は、ＡＴＤＣ着火の場合における、着火時期での筒内燃料密度と実測された着火遅れ期間長との関係を示したグラフである。
【図１２】図１２は、ＡＴＤＣ着火の場合において本実施装置により予測された着火遅れ期間長と実測された着火遅れ期間長との関係を示したグラフである。
【図１３】図１３は、ＢＴＤＣ着火の場合における熱発生率波形を概念的に示した図である。
【図１４】図１４は、ＡＴＤＣ着火の場合における熱発生率波形を概念的に示した図である。
【図１５】図１５は、本実施装置のＣＰＵが実行するルーチンを示したフローチャートである。
【図１６】図１６は、種々の負荷率に対する熱発生率波形を概念的に示した図である。
【図１７】図１７は、種々のＥＧＲ率に対する熱発生率波形を概念的に示した図である。
【図１８】図１８は、種々の空燃比に対する熱発生率波形を概念的に示した図である。
【図１９】図１９は、種々の冷却水温に対する熱発生率波形を概念的に示した図である。
【図２０】図２０は、種々の点火時期に対する熱発生率波形を概念的に示した図である。
【図２１】図２１は、種々の機関回転速度に対する熱発生率波形を概念的に示した図である。
【図２２】図２２は、第１の型式の内燃機関について、本実施装置により予測された前半燃焼期間長と実測された前半燃焼期間長との関係を示したグラフである。
【図２３】図２３は、第２の型式の内燃機関について、本実施装置により予測された前半燃焼期間長と実測された前半燃焼期間長との関係を示したグラフである。
【図２４】図２４は、本実施装置のＣＰＵが実行するルーチンを示したフローチャートである。
【図２５】図２５は、点火時期を除く機関運転状態パラメータのうち負荷率のみが互いに異なる複数の機関運転状態において取得された熱発生率波形であって、熱発生率最大時期を互いに一致させるように点火時期を調整した熱発生率波形を重ねて表示したグラフである。
【図２６】図２６は、点火時期を変化させた場合の熱発生率波形を示したグラフである。
【図２７】図２７は、種々の機関回転速度における、熱発生率最大時期での筒内燃料密度と熱発生率傾きとの関係を示したグラフである。
【図２８】図２８は、本実施装置のＣＰＵが実行するルーチンを示したフローチャートである。
【図２９】図２９は、機関回転速度及び吸気弁作用角をそれぞれ一定に維持した場合における種々の吸気弁位相角に対する「実測された熱発生率最大時期と実測された前半燃焼期間長との関係」を示したグラフである。
【図３０】図３０の（Ａ）乃至（Ｃ）のそれぞれは、機関回転速度及び吸気弁作用角をそれぞれ一定に維持した場合における種々の吸気弁位相角に対する「実測された熱発生率最大時期と実測された前半燃焼期間長との関係」を示したグラフである。
【図３１】図３１の（Ａ）乃至（Ｃ）のそれぞれは、機関回転速度及び吸気弁作用角をそれぞれ一定に維持した場合における種々の吸気弁位相角に対する「実測された熱発生率最大時期と実測された前半燃焼期間長との関係」を示したグラフである。
【図３２】図３２は、機関回転速度及び吸気弁閉弁時期をそれぞれ一定に維持した場合における種々の吸気弁作用角に対する「実測された熱発生率最大時期と実測された前半燃焼期間長との関係」を示したグラフである。
【図３３】図３３の（Ａ）乃至（Ｃ）のそれぞれは、機関回転速度及び吸気弁閉弁時期をそれぞれ一定に維持した場合における種々の吸気弁作用角に対する「実測された熱発生率最大時期と実測された前半燃焼期間長との関係」を示したグラフである。
【図３４】図３４の（Ａ）乃至（Ｃ）のそれぞれは、機関回転速度及び吸気弁閉弁時期をそれぞれ一定に維持した場合における種々の吸気弁作用角に対する「実測された熱発生率最大時期と実測された前半燃焼期間長との関係」を示したグラフである。
【図３５】図３５は、本実施方法の適用例を示したルーチンである。
【図３６】図３６は、本実施方法の適用例を示したルーチンである。
【発明を実施するための形態】
【００４２】
以下、本発明の実施形態に係る燃焼状態推定方法（以下、単に「本実施方法」とも称呼する。）について図面を参照しながら説明する。本実施方法を使用するシミュレーション装置は、以下、単に「本実施装置」とも称呼する。本実施装置は、ＣＰＵ及びメモリ等含むコンピュータにより実現される。
【００４３】
図５に概略構成を示した本実施装置（シミュレーション装置）１０は、本実施方法を使用して、周知の「ガソリン燃料・ポート噴射式・ピストン往復動型・火花点火式・多気筒内燃機関（以下、単に「機関」とも称呼する。）」の燃焼状態を推定する。より具体的に述べると、本実施装置１０は、対象とする内燃機関の運転状態を表すパラメータ（以下、「機関運転状態パラメータ」とも称呼する。）ＥＰを入力し、その入力した機関運転状態パラメータに基づいて演算を行うことにより、クランク角θに対する熱発生率ｄＱ／ｄθを模擬的に示す波形（以下「熱発生率波形」とも称呼する）ＱＷを生成する。熱発生率波形ＱＷは、機関における燃料（混合気）の燃焼状態を良好に表す。
【００４４】
機関運転状態パラメータＥＰは、機関回転速度ＮＥ、負荷率ＫＬ、点火時期ＳＡ、ＥＧＲ率Ｇegr、機関に供給される混合気の空燃比Ａ／Ｆ、機関の冷却水温ＴＨＷ、吸気弁位相角ＩＮＶＴ及び吸気弁作用角ＶＣＡＭ等を含む。図５に示したように、本実施装置１０は、熱発生率波形ＱＷが実質的に三角形であるとして熱発生率波形を推定する。
【００４５】
熱発生率波形ＱＷを推定するために、本実施装置１０は、着火遅れ期間推定部（着火遅れ期間モデル）１１、前半燃焼期間推定部（前半燃焼期間推定モデル）１２、熱発生率傾き推定部（熱発生率傾きモデル）１３及び熱発生量推定部（熱発生量モデル）１４を備えている。これらの推定部（モデル）１１〜１４の機能は実際にはＣＰＵがプログラム（インストラクション）を実行することにより実現される。従って、図６に示したように、本実施装置１０は、ＣＰＵによって推定部１１〜１４の機能を順に実現する（即ち、ステップ６１０乃至ステップ６４０の処理を実行する）ことにより、着火遅れ期間の長さ（着火遅れ期間長さ、着火遅れ期間長）τ、前半燃焼期間の長さ（前半燃焼期間長）ａ、熱発生率傾きｂ／ａ及び熱発生量（Ｑ１、Ｑ２及びＱall）等を推定する。以下、これらの推定部について個別に説明する。
【００４６】
１．着火遅れ期間推定部
（着火遅れ期間推定部の概要）
着火遅れ期間は、図７に示したように、点火プラグから火花が発生させられる点火時期ＳＡから、燃焼室内に形成された混合気の燃料が実際に燃焼を開始する着火時期ＦＡまでの期間である。即ち、着火遅れ期間は、点火プラグの電極間に火花放電が行われる点火時期ＳＡから、この火花によって生成された火炎核が成長して爆発的な燃焼が開始する着火時期ＦＡまでの期間である。着火遅れ期間長τは、その着火遅れ期間の長さである。着火遅れ期間長τの単位はクランク角［ＣＡ］であるが、モデルの都合上［ｍｓ］の単位を有する期間として計算される場合がある。
【００４７】
なお、着火時期ＦＡは、点火時期ＳＡ以降に熱発生率ｄＱ／ｄθが所定閾値ｄｑｄθｔｈ（例えば、１［Ｊ／ＣＡ］）に達した時期であるとして定義される。但し、着火時期ＦＡは、点火時期ＳＡ以降に熱発生率ｄＱ／ｄθが０［Ｊ／ＣＡ］よりも大きくなった時期であるとして定義されてもよく、点火時期ＳＡ以降の熱発生量が総熱発生量の所定割合（例えば、５％）の熱発生量に達した時期であるとして定義されてもよい。その他、着火時期ＦＡは、当業者が認識する一般的な定義に従って定められた時期であってもよい。
【００４８】
着火遅れ期間推定部１１は、下記の（１）式及び（２）式の何れか一方を用いて着火遅れ期間長τ［ｍｓ］を推定する。但し、（１）式及び（２）式により推定される着火遅れ期間長τは以下の総ての条件が満たされる場合の着火遅れ期間長τである。これらの式に基づいて着火遅れ期間長τが精度良く推定され得る根拠については後述する。
【００４９】
（条件τ１）燃焼に供される混合気の空燃比Ａ／Ｆは理論空燃比（例えば、１４．６）である。
（条件τ２）ＥＧＲ率Ｇegrは「０」である。即ち、外部ＥＧＲは実行されていない。
（条件τ３）冷却水温ＴＨＷは機関完全暖機終了を示す冷却水温閾値Ｔｔｈ（例えば、８０℃）以上である。この条件τ３は機関が完全暖機後であることという条件であるので、冷却水温ＴＨＷに代わり機関の潤滑油の温度が条件τ３の判定に用いられてもよい。
（条件τ４）吸気弁の開弁特性（吸気弁開弁時期ＩＶＯ、吸気弁閉弁時期ＩＶＣ、吸気弁作用角ＶＣＡＭ、吸気弁位相角ＩＮＶＴ及び吸気弁リフト量ＩＶＬiftの最大値等）が吸気弁の基準特性（吸気弁基準特性）に設定されている。即ち、吸気弁開弁時期ＩＶＯが基準開弁時期ＩＶＯｒに、吸気弁閉弁時期ＩＶＣが基準閉弁時期ＩＶＣｒに、吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒに、吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒに、吸気弁リフト量ＩＶＬiftの最大値が基準最大値に、それぞれ設定されている。
（条件τ５）排気弁の開弁特性（排気弁開弁時期ＥＶＯ、排気弁閉弁時期ＥＶＣ、排気弁位相角ＥＸＶＴ及び排気弁リフト量ＥＶＬiftの最大値等）が排気弁の基準特性に設定されている。即ち、排気弁開弁時期ＥＶＯ、排気弁閉弁時期ＥＶＣ、排気弁位相角ＥＸＶＴ及び排気弁リフト量ＥＶＬiftのそれぞれが、それぞれの基準値に設定されている。
【００５０】
着火遅れ期間推定部１１は、着火時期ＦＡが、圧縮上死点よりも前の時期（クランク角が圧縮上死点よりも進角側のクランク角である時期）である場合、下記の（１）式に基づいて着火遅れ期間長τ［ｍｓ］を推定する。この場合を「ＢＴＤＣ着火」とも称呼する。
【００５１】
【数１】
2016160761.tif 000003
（１）式において、
Ｃ１は、実験等により予め適合された定数であり、
ρfuel@SAは、点火時期ＳＡにおける筒内（燃焼室内）の燃料密度（＝筒内燃料量［ｍｏｌ］／点火時期ＳＡでの燃焼室容積［Ｌ］）であり、
χ及びδは、何れも実験等により予め適合された定数である。
ＮＥは、機関回転速度（本明細書において同じ。）である。
【００５２】
着火遅れ期間推定部１１は、着火時期ＦＡが、圧縮上死点よりも後の時点（クランク角が圧縮上死点よりも遅角側のクランク角である時点）である場合、下記の（２）式に基づいて着火遅れ期間長τ［ｍｓ］を推定する。この場合を「ＡＴＤＣ着火」とも称呼する。
【００５３】
【数２】
2016160761.tif 000004
（２）式において、
Ｃ２は、実験等により予め適合された定数であり、
ρfuel@FAは、着火時期ＦＡにおける筒内の燃料密度（＝筒内燃料量［ｍｏｌ］／着火時期ＦＡでの燃焼室容積［Ｌ］］であり、
φ及びψは、何れも実験等により予め適合された定数である。
【００５４】
（着火遅れ期間推定部の作動）
ＣＰＵは、着火遅れ期間推定部１１に相当する機能を実現するために図８にフローチャートにより示したルーチンに従って着火遅れ期間長τ［ＣＡ］を推定する。その結果、ＣＰＵは更に着火時期ＦＡを推定する。なお、図８に示したルーチンにおいて、クランク角は圧縮上死点前（ＢＴＤＣ）のクランク角（圧縮上死点から進角側になるほど大きくなる正の値）により表される。
【００５５】
ＣＰＵは、ステップ８００から処理を開始してステップ８０５に進み、点火時期ＳＡ、燃料噴射量Ｆinj及び機関回転速度ＮＥ等の種々の機関運転状態パラメータを入力（取得）する。
【００５６】
次に、ＣＰＵはステップ８１０に進み、仮想着火時期ＦＡｚを初期値ＦＡ０に設定する。仮想着火時期ＦＡｚは着火遅れ期間長τを推定するために仮想的（暫定的）に設定される着火時期である。初期値ＦＡ０は、圧縮上死点前（ＢＴＤＣ）の所定のクランク角であり、点火時期ＳＡよりも遅角側であって且つ想定される着火時期のうちの最も進角側の着火時期ＦＡadよりも進角側のクランク角に予め定められている（ＳＡ＞ＦＡ０＞ＦＡad)。
【００５７】
次に、ＣＰＵはステップ８１５に進み、仮想着火時期ＦＡｚが圧縮上死点前であるか否かを判定する。ステップ８１５に最初に進んだとき、仮想着火時期ＦＡｚは「圧縮上死点前のクランク角である初期値ＦＡ０」である。従って、ＣＰＵはステップ８１５にて「Ｙｅｓ」と判定し、以下に述べるステップ８２０乃至ステップ８４０の処理を順に行う。
【００５８】
ステップ８２０：ＣＰＵは、燃料噴射量Ｆinjを点火時期ＳＡにおける筒内容積（燃焼室容積）Ｖ@SAにより除することによって「点火時期ＳＡにおける筒内燃料密度ρfuel@SA」を算出する。なお、燃焼室容積（筒内容積）Ｖはクランク角が定まれば一意に計算されるので、燃焼室容積Ｖ@SAも直ちに計算される。
ステップ８２５：ＣＰＵは、上記（１）式と同じ式に筒内燃料密度ρfuel@SA及び機関回転速度ＮＥを代入することによって、仮想着火遅れ期間長τz［ｍｓ］を算出する。
【００５９】
ステップ８３０：ＣＰＵは、仮想着火遅れ期間長τzの単位［ｍｓ］を機関回転速度ＮＥに基づいてクランク角の単位［ＣＡ］に変換し、仮想着火遅れ期間長τz［ＣＡ］を算出する。
ステップ８３５：ＣＰＵは、仮想着火時期ＦＡｚ［ＣＡ］に仮想着火遅れ期間長τz［ＣＡ］を加えることにより、仮想点火時期ＳＡｚ［ＣＡ］を算出する（図７を参照。）。
ステップ８４０：ＣＰＵは、仮想点火時期ＳＡｚが点火時期ＳＡに一致しているか否かを判定する。具体的には、ＣＰＵは、仮想点火時期ＳＡｚが、「点火時期ＳＡから微小な正の所定クランク角dsaを減じたクランク角」よりも大きく且つ「点火時期ＳＡに前記所定クランク角dsaを加えたクランク角」よりも小さいか否かを判定する。
【００６０】
仮想点火時期ＳＡｚが点火時期ＳＡに一致している場合、着火遅れ期間長τを算出するために仮想的に設定された仮想着火時期ＦＡｚが実際の着火時期ＦＡに等しいと考えることができる。そこで、仮想点火時期ＳＡｚが点火時期ＳＡに一致している場合、ＣＰＵはステップ８４０にて「Ｙｅｓ」と判定し、以下に述べるステップ８４５及びステップ８５０の処理を順に行い、ステップ８９５に進んで本ルーチンを終了する。
【００６１】
ステップ８４５：ＣＰＵは、着火遅れ期間長τ［ＣＡ］として仮想着火遅れ期間長τz［ＣＡ］を採用する。即ち、ＣＰＵは、この時点で計算されている仮想着火遅れ期間長τz［ＣＡ］を推定された着火遅れ期間長τ［ＣＡ］として確定する。
ステップ８５０：ＣＰＵは、着火時期ＦＡとして仮想着火時期ＦＡｚを採用する。即ち、ＣＰＵは、この時点での仮想着火時期ＦＡｚを推定された着火時期（真の着火時期、実際の着火時期）ＦＡとして確定する。
【００６２】
一方、ＣＰＵがステップ８４０の処理を行う時点において、仮想点火時期ＳＡｚが点火時期ＳＡに一致していないとき、ＣＰＵはステップ８４０にて「Ｎｏ」と判定してステップ８５５に進み、仮想着火時期ＦＡｚを微小な所定値dfazだけ減少させる。即ち、仮想着火時期ＦＡｚを微小な所定値dfazだけ遅角側のクランク角へと変更する。その後、ＣＰＵはステップ８１５に進む。
【００６３】
このような処理が繰り返され、仮想着火時期ＦＡｚが初期値ＦＡ０から圧縮上死点へと次第に近づく期間において、仮想点火時期ＳＡｚが点火時期ＳＡに一致すれば、その時点での仮想着火遅れ期間長τzが着火遅れ期間長τとして取得され、その時点の仮想着火時期ＦＡｚが着火時期ＦＡとして取得される。
【００６４】
これに対し、仮想着火時期ＦＡｚが負の値（即ち、圧縮上死点後のクランク角）になるまで仮想点火時期ＳＡｚが点火時期ＳＡに一致しない場合、ＣＰＵはステップ８１５に進んだとき、そのステップ８１５にて「Ｎｏ」と判定し、以下に述べるステップ８６０及びステップ８６５の処理を順に行う。その後、ＣＰＵはステップ８３０以降に進む。
【００６５】
ステップ８６０：ＣＰＵは、燃料噴射量Ｆinjを仮想着火時期ＦＡｚおける燃焼室容積Ｖ@FAzにより除することによって「仮想着火時期ＦＡｚにおける筒内燃料密度ρfuel@FAz」を算出する。なお、燃焼室容積Ｖはクランク角が定まれば一意に計算されるので、燃焼室容積Ｖ@FAzも直ちに計算される。
ステップ８６５：ＣＰＵは、上記（２）式と同じ式に筒内燃料密度ρfuel@FAz及び機関回転速度ＮＥを代入することによって、仮想着火遅れ期間長τz［ｍｓ］を算出する。
【００６６】
以上の処理により、着火遅れ期間長τ［ＣＡ］と着火時期ＦＡが推定（取得）される。なお、着火時期ＦＡが圧縮上死点後のクランク角により表される場合、ステップ８５０にて推定された着火時期ＦＡの符号が反転される。即ち、着火時期ＦＡ（ATDC deg）は「−ＦＡ（BTDC deg)」である。
【００６７】
（着火遅れ期間推定部の正当性）
次に、上記（１）式及び（２）式に基づくことにより、着火遅れ期間長τを精度良く推定できる点について説明する。即ち、（１）式及び（２）式が着火遅れ期間モデルとして適切なモデルである点について説明する。
【００６８】
・上記（１）式について（ＢＴＤＣ着火の場合）
図９はＢＴＤＣ着火であり且つ上記条件τ１〜τ５が総て成立している場合における、筒内燃料密度ρfuel@SAと着火遅れ期間長τ［ｍｓ］との関係を種々の機関回転速度ＮＥに対して実測した結果を示すグラフである。
【００６９】
図９から理解されるように、ＢＴＤＣ着火の場合、筒内燃料密度ρfuel@SAと着火遅れ期間長τ［ｍｓ］との間には機関回転速度ＮＥ別に強い相関がある。その相関関係は上記（１）式の関数形式により表される。実際、図９においては、機関回転速度ＮＥが１２００［ｒｐｍ］である場合の筒内燃料密度ρfuel@SAと着火遅れ期間長τ［ｍｓ］との関係を（１）式に基づいて算出した曲線と、機関回転速度ＮＥが２４００［ｒｐｍ］である場合の筒内燃料密度ρfuel@SAと着火遅れ期間長τ［ｍｓ］との関係を（１）式に基づいて算出した曲線と、が破線により描かれている。
【００７０】
物理的な観点から図９のグラフについて説明を加えると、ＢＴＤＣ着火においては筒内燃料密度ρfuel@SAが高いほど着火遅れ期間長τ［ｍｓ］は短くなっている。これは、筒内燃料密度ρfuel@SAが高いほど、点火時期ＳＡにおける点火プラグ周辺の燃料分子の数が多くなるから、点火用火花の発生後における火炎核の成長が急速に行われるようになるためであると推測される。一方、機関回転速度ＮＥが高いほど着火遅れ期間長τ［ｍｓ］は短くなっている。これは、機関回転速度ＮＥが高いほど、筒内における気流の乱れの強さがより強くなるから、点火用火花の発生後における火炎核の成長が急速に行われるようになるためであると推測される。
【００７１】
図１０は、（１）式により算出（予測）された着火遅れ期間長（予測着火遅れ期間長）と、実測された着火遅れ期間長（実測着火遅れ期間長）と、の関係を検証した結果を示すグラフである。この図１０からも明からなように、（１）式により算出される予測着火遅れ期間長は実測着火遅れ期間長に精度良く一致している。即ち、（１）式がＢＴＤＣ着火における着火遅れ期間長τの推定に適した式（着火遅れ期間モデル）であることが理解される。
【００７２】
・上記（２）式について（ＡＴＤＣ着火の場合）
図１１はＡＴＤＣ着火であり且つ上記条件τ１〜τ５が総て成立している場合であって機関回転速度ＮＥが１２００［ｒｐｍ］であるときに、筒内燃料密度ρfuel@FAと着火遅れ期間長τ［ｍｓ］との関係を種々の負荷率ＫＬに対して実測した結果を示すグラフである。なお、負荷率ＫＬは、空気充填率あり、着目する気筒が一回の吸気行程において吸入する空気量をＭｃ［ｇ］、空気密度をρ［ｇ／Ｌ］）、機関の排気量をＬｖ［Ｌ］、機関の気筒数を「４」であるとしたとき、次式により算出される。
ＫＬ＝｛Ｍｃ／（ρ・Ｌｖ／４）｝・100（％）
【００７３】
図１１から理解されるように、ＡＴＤＣ着火の場合、筒内燃料密度ρfuel@FAと着火遅れ期間長τ［ｍｓ］との間には負荷率ＫＬに依らず（即ち、機関の負荷に関わらず）強い相関がある。その相関関係は上記（２）式の関数形式により表される。
【００７４】
物理的な観点から図１１のグラフについて説明を加えると、ＡＴＤＣ着火においては筒内燃料密度ρfuel@FAが高いほど着火遅れ期間長τ［ｍｓ］は短くなっている。これは、ＢＴＤＣ着火の場合と同様、筒内燃料密度ρfuel@FAが高いほど、着火時期ＦＡの直前の期間における点火プラグ周辺の燃料分子の数が多くなるから、火炎核の成長が急速に行われるようになるためであると推測される。更に、図示されていないが、機関回転速度ＮＥが高いほど、着火遅れ期間長τ［ｍｓ］は短くなった。これは、機関回転速度ＮＥが高いほど、筒内における気流の乱れの強さがより強くなるから、点火用火花の発生後における火炎核の成長が急速に行われるようになるためであると推測される。従って、機関回転速度ＮＥも着火遅れ期間長τ［ｍｓ］を決定するパラメータであると推測される。
【００７５】
図１２は、（２）式により算出（予測）された予測着火遅れ期間長と、実測着火遅れ期間長と、の関係を検証した結果を示すグラフである。この図１２からも明からなように、（２）式により算出される予測着火遅れ期間長は実測着火遅れ期間長に精度良く一致している。即ち、（２）式がＡＴＤＣ着火における着火遅れ期間長τの推定に適した式（着火遅れ期間モデル）であることが理解される。
【００７６】
・ＢＴＤＣ着火とＡＴＤＣ着火とで着火遅れ期間モデル（（１）式及び（２）式）を使い分ける理由
着火遅れ期間においては点火用火花により生成された火炎核が成長している。一方、着火遅れ期間において燃焼室容積は時々刻々変化しているから、火炎核の成長に相関が強い筒内燃料密度ρfuelも時々刻々変化している。従って、本来的には時々刻々変化する筒内燃料密度ρfuelを変数として有する着火遅れ期間のモデル式により着火遅れ期間の長さを推定することが好適であると考えられる。しかしながら、そのような着火遅れ期間のモデル式は複雑化する。そこで、発明者は、着火遅れ期間における燃料密度ρfuelの平均的な値を代表する値として「特定の時期での燃料密度ρfuel」を着火遅れ期間のモデル式の変数として採用することを検討した。
【００７７】
ところで、ＢＴＤＣ着火の場合、図１３に示したように、着火時期ＦＡのみならず点火時期ＳＡも当然に圧縮上死点前である。よって、ＢＴＤＣ着火においては、着火遅れ期間中に燃焼室容積が小さくなって行き（単調に減少し）、これに伴って筒内燃料密度ρfuelは高くなって行く（単調に増加する。）。
【００７８】
これに対し、ＡＴＤＣ着火の場合、図１４の（Ａ）及び（Ｂ）に示したように、点火時期ＳＡが圧縮上死点前である場合（図１４の（Ａ）を参照。）と、点火時期ＳＡが圧縮上死点後である場合（図１４の（Ｂ）を参照。）との２つの場合が存在する。
【００７９】
図１４の（Ｂ）に示した場合、点火時期ＳＡ以降においてピストンは膨張下死点に向かって移動して行く。つまり、着火遅れ期間において、燃焼室容積は大きくなって行き（単調に増加し）、これに伴って筒内燃料密度ρfuelは低くなって行く（単調に減少する）。また、図１４の（Ａ）に示した場合であっても、殆どの場合において点火時期ＳＡと圧縮上死点とのクランク角差の大きさの方が着火時期ＦＡと圧縮上死点とのクランク角差の大きさよりも小さい。更に、クランク角が上死点近傍にある場合、クランク角の変化に対する燃焼室容積の変化は非常に小さく、その後、クランク角が上死点近傍から離れるとクランク角の変化に対する燃焼室容積は急激に増加する。よって、ＡＴＤＣ着火においては、着火遅れ期間中に筒内燃料密度ρfuelは単調に減少すると考えて良い。
【００８０】
これらのことから、ＢＴＤＣ着火の場合（つまり、筒内燃料密度ρfuelが単調に増加する場合）、着火遅れ期間における燃料密度ρfuelの平均的な値を代表する値としては、「点火時期での燃料密度ρfuel@SA」が適切であると考えられる。更に、ＡＴＤＣ着火の場合（つまり、筒内燃料密度ρfuelが単調に減少する場合）、着火遅れ期間における燃料密度ρfuelの平均的な値を代表する値としては、「着火時期での燃料密度ρfuel@FA」が適切であると考えられる。以上が、（１）式と（２）式とを使い分ける理由である。
【００８１】
２．前半燃焼期間推定部
（前半燃焼期間推定部の概要）
前半燃焼期間は、図７に示したように、着火時期ＦＡから熱発生率ｄＱ／ｄθが最大となる時期（即ち、熱発生率最大時期θdQpeak）までの期間である。前半燃焼期間長ａは、その前半燃焼期間の長さである。熱発生率最大時期θdQpeakの単位は圧縮上死点後のクランク角［deg］であり、前半燃焼期間長ａの単位はクランク角（クランク角の幅）［ＣＡ］である。
【００８２】
前半燃焼期間推定部１２は、下記の（３）式を用いて前半燃焼期間長ａ［ＣＡ］を推定する。但し、（３）式により推定される前半燃焼期間長ａは以下の条件（条件ａ１及びａ２）が満たされる場合の前半燃焼期間長ａである。換言すると、（３）式は、負荷率ＫＬ、ＥＧＲ率Ｇegr、空燃比Ａ／Ｆ及び冷却水温ＴＨＷ（機関暖機状態）に関わらず成立する。（３）式に基づいて前半燃焼期間長ａが精度良く推定され得る根拠については後述する。更に、下記の（条件ａ１）が成立しない場合、即ち、吸気弁の開弁特性が吸気弁の基準特性（吸気弁基準特性）と相違した特性（特定特性：吸気弁特定特性）へと変化した場合の前半燃焼期間長ａの推定方法については後述する。
【００８３】
（条件ａ１）吸気弁の開弁特性（吸気弁開弁時期ＩＶＯ、吸気弁閉弁時期ＩＶＣ、吸気弁作用角ＶＣＡＭ、吸気弁位相角ＩＮＶＴ及び吸気弁リフト量ＩＶＬiftの最大値等）が吸気弁の基準特性（吸気弁基準特性）に設定されている。即ち、吸気弁開弁時期ＩＶＯが基準開弁時期ＩＶＯｒに、吸気弁閉弁時期ＩＶＣが基準閉弁時期ＩＶＣｒに、吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒに、吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒに、吸気弁リフト量ＩＶＬiftの最大値が基準最大値に、それぞれ設定されている。
（条件ａ２）排気弁の開弁特性（排気弁開弁時期ＥＶＯ、排気弁閉弁時期ＥＶＣ、排気弁位相角ＥＸＶＴ及び排気弁リフト量ＥＶＬiftの最大値等）が排気弁の基準特性に設定されている。即ち、排気弁開弁時期ＥＶＯ、排気弁閉弁時期ＥＶＣ、排気弁位相角ＥＸＶＴ及び排気弁リフト量ＥＶＬiftの最大値のそれぞれが、それぞれの基準値に設定されている。なお、実際には、この条件ａ２は必須ではなく、排気弁の開弁特性は任意の特性であってもよい。
【００８４】
【数３】
2016160761.tif 000005
（３）式において、
Ｃ３は、実験等により予め適合された定数であり、
Ｖ@θdQpeakは、熱発生率最大時期θdQpeakにおける燃焼室容積［Ｌ］であり、
α及びβは、何れも実験等により予め適合された定数である。
【００８５】
（前半燃焼期間推定部の作動）
ＣＰＵは、前半燃焼期間推定部１２に相当する機能を実現するために図１５にフローチャートにより示したルーチンに従って前半燃焼期間長ａ［ＣＡ］を推定する。なお、図１５に示したルーチンにおいて、クランク角は圧縮上死点後（ＡＴＤＣ）のクランク角（圧縮上死点から遅角側になるほど大きくなる正の値［ATDC deg］により表される。
【００８６】
ＣＰＵは、ステップ１５００から処理を開始し、以下に述べるステップ１５０５乃至ステップ１５３０の処理を順に行う。
【００８７】
ステップ１５０５：ＣＰＵは、先に推定された着火時期ＦＡ及び機関回転速度ＮＥ等の種々のパラメータを入力（取得）する。
ステップ１５１０：ＣＰＵは、仮想熱発生率最大時期θzdQpeakを初期値θ０に設定する。仮想熱発生率最大時期θzdQpeakは前半燃焼期間長ａを推定するために仮想的（暫定的）に設定される熱発生率最大時期である。ここでは、初期値θ０は、着火時期ＦＡに対して「小さな正の所定クランク角Δθ」だけ大きい（遅角側の）クランク角に予め定められている（ＦＡ＜θ０＝ＦＡ＋Δθ）。
【００８８】
ステップ１５１５：ＣＰＵは、仮想熱発生率最大時期θzdQpeakにおける燃焼室容積Ｖ@θzdQpeakを算出する。前述したように、燃焼室容積Ｖはクランク角が定まれば一意に計算されるので、燃焼室容積Ｖ@θzdQpeakも直ちに計算される。
ステップ１５２０：ＣＰＵは、上記（３）式と同じ式に燃焼室容積Ｖ@θzdQpeak及び機関回転速度ＮＥを代入することによって、仮想前半燃焼期間長ａz［ＣＡ］を算出する。
【００８９】
ステップ１５２５：ＣＰＵは、仮想熱発生率最大時期θzdQpeakから仮想前半燃焼期間長ａzを減じることにより、仮想着火時期ＦＡｘを算出する（図７を参照。）。
ステップ１５３０：ＣＰＵは、仮想着火時期ＦＡｘが着火時期ＦＡに一致しているか否かを判定する。具体的には、ＣＰＵは、仮想着火時期ＦＡｘが、「着火時期ＦＡから微小な正の所定クランク角dfaを減じたクランク角」よりも大きく且つ「着火時期ＦＡに前記クランク角dfaを加えたクランク角」よりも小さいか否かを判定する。
【００９０】
仮想着火時期ＦＡｘが着火時期ＦＡに一致している場合、仮想前半燃焼期間長ａzを算出するために仮想的に設定された仮想熱発生率最大時期θzdQpeakが実際の（真の）熱発生率最大時期θdQpeakに等しいと考えることができる。そこで、仮想着火時期ＦＡｘが着火時期ＦＡに一致している場合、ＣＰＵはステップ１５３０にて「Ｙｅｓ」と判定し、以下に述べるステップ１５３５及びステップ１５４０の処理を順に行い、ステップ１５９５に進んで本ルーチンを終了する。
【００９１】
ステップ１５３５：ＣＰＵは、熱発生率最大時期θdQpeakとして仮想熱発生率最大時期θzdQpeakを採用する。。即ち、ＣＰＵは、この時点での仮想熱発生率最大時期θzdQpeakを推定された（真の）熱発生率最大時期θdQpeakとして確定する。
ステップ１５４０：ＣＰＵは、前半燃焼期間長ａとして仮想前半燃焼期間長ａzを採用する。即ち、ＣＰＵは、この時点での仮想前半燃焼期間長ａzを推定された（真の）前半燃焼期間長ａとして確定する。
【００９２】
一方、ＣＰＵがステップ１５３０の処理を行う時点において、仮想着火時期ＦＡｘが着火時期ＦＡに一致していないとき、ＣＰＵはそのステップ１５３０にて「Ｎｏ」と判定して１５４５に進み、仮想熱発生率最大時期θzdQpeakを微小な正の所定値ddqpだけ増大させる。即ち、仮想熱発生率最大時期θzdQpeaを所定値ddqpだけ遅角側のクランク角へと変更する。その後、ＣＰＵはステップ１５１５に進む。
【００９３】
このような処理が繰り返され、仮想熱発生率最大時期θzdQpeakが次第に膨張下死点へと近づく期間において、仮想着火時期ＦＡｘが着火時期ＦＡに一致すれば、その時点での仮想熱発生率最大時期θzdQpeakが熱発生率最大時期θdQpeakとして採用され、その時点の仮想前半燃焼期間長ａzが前半燃焼期間長ａとして採用される。
【００９４】
（前半燃焼期間推定部の正当性）
次に、上記（３）式に基づくことにより、前半燃焼期間長ａを精度良く推定できる点について説明する。即ち、（３）式が前半燃焼期間モデルとして適切なモデルである点について説明する。
【００９５】
図１６は、点火時期ＳＡを除く機関運転状態パラメータのうち負荷率ＫＬのみが互いに異なる複数の機関運転状態のそれぞれにおいて取得された熱発生率波形であって、熱発生率最大時期θdQpeakを互いに一致させるように点火時期ＳＡを調整した熱発生率波形を重ねて表示したグラフである。
図１７は、点火時期ＳＡを除く機関運転状態パラメータのうちＥＧＲ率Ｇegrのみが互いに異なる複数の機関運転状態において取得された熱発生率波形であって、熱発生率最大時期θdQpeakを互いに一致させるように点火時期ＳＡを調整した熱発生率波形を重ねて表示したグラフである。
図１８は、点火時期ＳＡを除く機関運転状態パラメータのうち空燃比Ａ／Ｆのみが互いに異なる複数の機関運転状態において取得された熱発生率波形であって、熱発生率最大時期θdQpeakを互いに一致させるように点火時期ＳＡを調整した熱発生率波形を重ねて表示したグラフである。
図１９は、点火時期ＳＡを除く機関運転状態パラメータのうち冷却水温ＴＨＷのみが互いに異なる複数の機関運転状態（即ち、機関暖機状態のみが互いに異なる運転状態）において取得された熱発生率波形であって、熱発生率最大時期θdQpeakを互いに一致させるように点火時期ＳＡを調整した熱発生率波形を重ねて表示したグラフである。
【００９６】
これらの図１６乃至図１９から理解されるように、前半燃焼期間長ａは、負荷率ＫＬ、ＥＧＲ率Ｇegr、空燃比Ａ／Ｆ及び冷却水温ＴＨＷの何れが変化しても一定に維持されている。換言すると、前半燃焼期間長ａは、負荷率ＫＬ、ＥＧＲ率Ｇegr、空燃比Ａ／Ｆ及び冷却水温ＴＨＷ（即ち、機関暖機状態）に依存せず、これらの影響を受けないことが理解される。
【００９７】
これに対し、図２０は、機関運転状態パラメータのうち点火時期ＳＡのみが互いに異なる複数の機関運転状態において取得された熱発生率波形を重ねて表示したグラフである。図２０から、点火時期ＳＡが遅角されるほど、前半燃焼期間長ａは長くなることが理解される。
【００９８】
これは、点火時期ＳＡが、着火時期ＦＡから熱発生率最大時期θdQpeakに到るまでの間の筒内の気流の乱れの強さ（以下、「前半燃焼期間中の気流乱れの程度」とも称呼する。）に相関を有するためであると考えられる。即ち、点火時期ＳＡがより遅角側に移行されるほど、着火時期ＦＡ及び熱発生率最大時期θdQpeakも遅角側に移行する。その結果、前半燃焼期間中の気流乱れの程度が弱くなり、且つ、熱発生率最大時期θdQpeakにおける燃焼室容積Ｖ@θdQpeakが大きくなる。従って、点火時期ＳＡがより遅角側に移行されるほど、火炎伝搬速度が小さくなるから、前半燃焼期間長ａは長くなると推察される。
【００９９】
更に、図２１は、点火時期ＳＡを除く機関運転状態パラメータのうち機関回転速度ＮＥのみが互いに異なる複数の機関運転状態において取得された熱発生率波形であって、熱発生率最大時期θdQpeakを互いに一致させるように点火時期ＳＡを調整した熱発生率波形を重ねて表示したグラフである。図２１から理解されるように、機関回転速度ＮＥが高いほど前半燃焼期間長ａ［ＣＡ］は長くなっている。このことから、機関回転速度ＮＥが高いほど前半燃焼期間長ａ［ＣＡ］を長くする要因が存在していると推察される。
【０１００】
この要因について検討すると、機関回転速度ＮＥが前半燃焼期間中の気流乱れの程度に相関を有することが考えられる。一般に、機関回転速度ＮＥが高いほど吸気系から気筒内に流入する空気の流速が高くなるので、燃焼室内に形成される気流の乱れの程度は大きくなる。しかしながら、その気流の乱れの程度により前半燃焼期間長ａ［ＣＡ］が短くなる程度（割合）は機関回転速度ＮＥの増大程度（割合）と反比例の関係にない。即ち、例えば、機関回転速度ＮＥが２倍となったとしても、前半燃焼期間長ａ［ＣＡ］は（１／２）倍にはならず、（１／２）倍した値よりも長くなる。従って、機関回転速度ＮＥが高いほど前半燃焼期間長ａ［ＣＡ］は長くなり、機関回転速度ＮＥが低いほど前半燃焼期間長ａ［ＣＡ］は短くなると推察される。
【０１０１】
以上の検討の結果、発明者は、前半燃焼期間長ａに影響を与えるパラメータ（前半燃焼期間長ａを推定するための主要なパラメータ）として、点火時期ＳＡに相関のある物理量としての「熱発生率最大時期θdQpeakにおける燃焼室容積Ｖ@θdQpeak」と、「機関回転速度ＮＥ」と、を採用し、上記（３）式を得た。
【０１０２】
ここで、上記（３）式における各係数について説明を加える。（３）式における「Ｃ３及びα」は実験等に基いて同定される。（３）式におけるβはタンブル比が大きいほど多くなる値として決定されるが、実験等に基いて同定されてもよい。
【０１０３】
図２２及び図２３は、互いに異なる内燃機関に対して、（３）式を用いて算出された前半燃焼期間（予測前半燃焼期間）と、実際に測定された前半燃焼期間（実測前半燃焼期間）と、の関係を検証した結果を示すグラフである。これらの図２２及び図２３からも明からなように、（３）式により算出される予測前半燃焼期間は実測前半燃焼期間に精度良く一致している。即ち、（３）式が前半燃焼期間長ａの推定に適した式（前半燃焼期間モデル）であることが理解される。
【０１０４】
３．熱発生率傾き推定部
（熱発生率傾き推定部の概要）
熱発生率傾きｂ／ａ（＝Ｇ）は、図７に示したように、前半燃焼期間長ａにおける熱発生率ｄＱ／ｄθの傾き（前半燃焼期間における熱発生率の増加率の平均値）である。熱発生率傾きｂ／ａの単位は［Ｊ／ＣＡ２］である。
【０１０５】
熱発生率傾き推定部１３は、下記の（４）式を用いて熱発生率傾きｂ／ａを推定する。但し、ここで推定される熱発生率傾きｂ／ａは以下の総ての条件が満たされる場合の熱発生率傾きｂ／ａである。（４）式に基づいて熱発生率傾きｂ／ａが精度良く推定され得る根拠については後述する。
【０１０６】
（条件Ｇ１）燃焼に供される混合気の空燃比Ａ／Ｆは理論空燃比（例えば、１４．６）である。
（条件Ｇ２）ＥＧＲ率Ｇegrは「０」である。即ち、外部ＥＧＲは実行されていない。
（条件Ｇ３）冷却水温ＴＨＷは機関完全暖機終了を示す冷却水温閾値Ｔｔｈ（例えば、８０℃）以上である。この条件Ｇ３は機関が完全暖機後であることという条件であるので、冷却水温ＴＨＷに代わり機関の潤滑油の温度が条件Ｇ３の判定に用いられてもよい。
（条件Ｇ４）吸気弁の開弁特性（吸気弁開弁時期ＩＶＯ、吸気弁閉弁時期ＩＶＣ、吸気弁作用角ＶＣＡＭ、吸気弁位相角ＩＮＶＴ及び吸気弁リフト量ＩＶＬiftの最大値等）が吸気弁の基準特性（吸気弁基準特性）に設定されている。即ち、吸気弁開弁時期ＩＶＯが基準開弁時期ＩＶＯｒに、吸気弁閉弁時期ＩＶＣが基準閉弁時期ＩＶＣｒに、吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒに、吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒに、吸気弁リフト量ＩＶＬiftの最大値が基準最大値に、それぞれ設定されている。
（条件Ｇ５）排気弁の開弁特性（排気弁開弁時期ＥＶＯ、排気弁閉弁時期ＥＶＣ、排気弁位相角ＥＸＶＴ及び排気弁リフト量ＥＶＬiftの最大値等）が排気弁の基準特性に設定されている。即ち、排気弁開弁時期ＥＶＯ、排気弁閉弁時期ＥＶＣ、排気弁位相角ＥＸＶＴ及び排気弁リフト量ＥＶＬiftの最大値のそれぞれが、それぞれの基準値に設定されている。
【０１０７】
【数４】
2016160761.tif 000006
（４）式において、
Ｃ４及びγは、実験等により予め適合された定数であり、
ρfuel@θdQpeakは、燃焼が生じないとした場合における熱発生率最大時期θdQpeakでの筒内（燃焼室内）の燃料密度（＝筒内燃料量［ｍｏｌ］／熱発生率最大時期θdQpeakでの燃焼室容積［Ｌ］）である。
【０１０８】
（熱発生率傾き推定部の作動）
ＣＰＵは、熱発生率傾き推定部１３に相当する機能を実現するために図２４にフローチャートにより示したルーチンに従って熱発生率傾きｂ／ａを推定する。なお、図２４に示したルーチンにおいて、クランク角は圧縮上死点後（ＡＴＤＣ）のクランク角［ATDC deg］により表される。
【０１０９】
ＣＰＵは、ステップ２４００から処理を開始し、以下に述べるステップ２４１０乃至ステップ２４４０の処理を順に行い、ステップ２４９５に進んで本ルーチンを終了する。
【０１１０】
ステップ２４１０：ＣＰＵは、先に推定された「着火時期ＦＡ、熱発生率最大時期θdQpeak」及び燃料噴射量Ｆinj等の種々のパラメータを入力（取得）する。
ステップ２４２０：ＣＰＵは、熱発生率最大時期θdQpeakにおける燃焼室容積Ｖ@θdQpeakを算出する。
ステップ２４３０：ＣＰＵは、燃料噴射量Ｆinjを燃焼室容積Ｖ@θdQpeakにより除することによって「熱発生率最大時期θdQpeakにおける筒内燃料密度ρfuel@θdQpeak」を算出する。
ステップ２４４０：ＣＰＵは、上記（４）式と同じ式に筒内燃料密度ρfuel@θdQpeakを代入することによって、熱発生率傾きｂ／ａを算出する。以上により、熱発生率傾きｂ／ａが推定される。
【０１１１】
（熱発生率傾き推定部の正当性）
次に、上記（４）式に基づくことにより、熱発生率傾きｂ／ａを精度良く推定できる点について説明する。即ち、（４）式が熱発生率傾きモデルとして適切なモデルである点について説明する。
【０１１２】
図２５の（Ａ）乃至（Ｄ）のそれぞれは、点火時期ＳＡを除く機関運転状態パラメータのうち負荷率ＫＬのみが互いに異なる複数の機関運転状態において取得された熱発生率波形であって、熱発生率最大時期θdQpeakを互いに一致させるように点火時期ＳＡを調整した熱発生率波形を重ねて表示したグラフである。図２５の（Ａ）から（Ｄ）の順に点火時期ＳＡは遅角側に移行されている。更に、図２５の（Ａ）乃至（Ｄ）において、負荷率ＫＬは、ＫＬ１、ＫＬ２、ＫＬ３の順に大きい。
【０１１３】
更に、図２６の（Ａ）及び（Ｂ）のそれぞれは、点火時期ＳＡのみが互いに異なる複数の機関運転状態において取得された熱発生率波形を表示したグラフである。
【０１１４】
図２５の（Ａ）乃至（Ｄ）、並びに、図２６の（Ａ）及び（Ｂ）から理解されるように、熱発生率傾きｂ／ａは、負荷率ＫＬが大きいほど大きく、且つ、点火時期ＳＡが進角側であるほど大きい。
【０１１５】
負荷率ＫＬが大きいほど大きいほど熱発生率傾きｂ／ａが大きくなる理由は、負荷率ＫＬが高いほど燃料噴射量Ｆinjも大きくなり、従って、前半燃焼期間における筒内燃料密度が高くなり、その結果、混合気の燃焼速度（火炎伝搬速度）が大きくなるからであると推察される。
【０１１６】
点火時期ＳＡが進角側であるほど熱発生率傾きｂ／ａが大きくなる理由は、負荷率ＫＬの場合と同様、点火時期ＳＡが進角側であるほど前半燃焼期間における筒内燃料密度が高くなり、その結果、混合気の燃焼速度（火炎伝搬速度）が大きくなるからであると推察される。即ち、ピストンが圧縮上死点近傍にあるときには、クランク角の変化に伴う燃焼室容積の変化は小さいが、膨張行程においては（特に、クランク角が圧縮上死点後１０［deg］近傍より大きくなると）燃焼室容積が急激に増大して行くので、前半燃焼期間における筒内燃料密度が低くなり、その結果、混合気の燃焼速度（火炎伝搬速度）が小さくなると推察される。
【０１１７】
そこで、発明者は、前半燃焼期間における筒内燃料密度の代表値としての「熱発生率最大時期θdQpeakにおける筒内燃料密度ρfuel@θdQpeak」と、熱発生率傾きb／aと、の関係を「種々の負荷率ＫＬ及び種々の機関回転速度ＮＥ」に対して実験により調べた。その結果を、図２７の（Ａ）乃至（Ｄ）に示す。
【０１１８】
図２７の（Ａ）乃至（Ｄ）から、機関回転速度ＮＥが一定であれば、負荷率ＫＬ（従って、点火時期ＳＡ）が相違しても、筒内燃料密度ρfuel@θdQpeakと熱発生率傾きb／aとは概ね比例関係にあることが理解される。従って、（４）式により算出される熱発生率傾きは実測熱発生率傾きと精度良く一致することが理解される。即ち、（４）式が熱発生率傾きｂ／ａの推定に適した式（熱発生率傾きモデル）であることが理解される。
【０１１９】
４．熱発生量推定部
（熱発生量推定部の概要）
総熱発生量Ｑallは、図７に示したように、前半燃焼期間における熱発生量Ｑ１と、後半燃焼期間における熱発生量Ｑ２と、の和である。即ち、総熱発生量Ｑallは、一回の燃焼において発生した総熱量である。これらの熱発生量の単位は[Ｊ]である。なお、後半燃焼期間は、熱発生率最大時期θdQpeakから燃焼が終了する燃焼終了時期ＥＡまでの期間である。
【０１２０】
熱発生量推定部１４は、下記の（５）式乃至（９）式を用いて、熱発生率最大時期θdQpeakにおける熱発生率ｄＱ／ｄθ@θdQpeak（＝ｂ）、総熱発生量Ｑall、前半燃焼期間における熱発生量Ｑ１、後半燃焼期間における熱発生量Ｑ２及び後半燃焼期間長ｃを推定する。なお、（５）式には推定した「前半燃焼期間長ａ及び熱発生率傾きｂ／ａ」を適用する。また、（７）式におけるＫは燃焼効率である。燃焼効率Ｋは、負荷率ＫＬ、機関回転速度ＮＥ及び冷却水温ＴＨＷと燃焼効率Ｋとの関係を予め実験により測定し、その測定したデータに基づいてこれらの関係を定めるルックアップテーブルＭａｐＫ（ＫＬ，ＮＥ，ＴＨＷ）を作成しておき、そのテーブルＭａｐＫ（ＫＬ，ＮＥ，ＴＨＷ）に実際の「負荷率ＫＬ、機関回転速度ＮＥ及び冷却水温ＴＨＷ」を適用することにより求められる。Ｆinjは上述したように実際の燃料噴射量である。（６）式及び（９）式は図７に示した三角形の面積を求める式から容易に導き出せる。
【数５】
2016160761.tif 000007
【数６】
2016160761.tif 000008
【数７】
2016160761.tif 000009
【数８】
2016160761.tif 000010
【数９】
2016160761.tif 000011
【０１２１】
（熱発生量推定部の作動）
ＣＰＵは、熱発生量推定部１４に相当する機能を実現するために図２８にフローチャートにより示したルーチンに従って総熱発生量Ｑall及び後半燃焼期間等を推定する。即ち、ＣＰＵは、ステップ２８００から処理を開始し、以下に述べるステップ２８１０乃至ステップ２８７０の処理を順に行い、ステップ２８９５に進んで本ルーチンを終了する。
【０１２２】
ステップ２８１０：ＣＰＵは、先に推定された「熱発生率傾きｂ／ａ及び前半燃焼期間長ａ」並びに「冷却水温ＴＨＷ、負荷率ＫＬ、燃料噴射量Ｆinj及び機関回転速度ＮＥ」等の種々のパラメータを入力（取得）する。
ステップ２８２０：ＣＰＵは、上記（５）式により熱発生率最大時期θdQpeakにおける熱発生率ｂを算出する。
ステップ２８３０：ＣＰＵは、上記（６）式により前半燃焼期間における熱発生量Ｑ１を算出する。
ステップ２８４０：ＣＰＵは、ルックアップテーブルＭａｐＫ（ＫＬ，ＮＥ，ＴＨＷ）に実際の「負荷率ＫＬ、機関回転速度ＮＥ及び冷却水温ＴＨＷ」を適用することにより、燃焼効率Ｋを求める。
ステップ２８５０：ＣＰＵは、上記（７）式により総熱発生量Ｑallを算出する。
ステップ２８６０：ＣＰＵは、上記（８）式により後半燃焼期間における熱発生量Ｑ２を算出する。
ステップ２８７０：ＣＰＵは、上記（９）式により後半燃焼期間長ｃを算出する。
【０１２３】
以上から、本実施装置１０は本実施方法を用いて熱発生率波形ＱＷを推定することができる。
【０１２４】
（吸気弁の開弁特性が変化した場合の前半燃焼期間の推定）
ところで、上記（３）式の前提条件の一つである条件ａ１は、吸気弁の開弁特性が吸気弁の基準特性に設定されているときに成立する条件であった。換言すると、条件ａ１が成立しない場合（即ち、吸気弁開弁特性が吸気弁基準特性と相違する場合）、上記（３）式における値「Ｃ３、α及びβ」を決定し直さなければならない。その結果、各種の吸気弁の開弁特性に対する前半燃焼期間長ａを推定するためには、膨大な時間及び労力を必要としていた。
【０１２５】
これに対し、図１及び図２を参照しながら説明したように、本発明者は以下の知見１を得た。
（知見１）
機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭ０であり且つ吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒである場合における「熱発生率最大時期θdQpeakと前半燃焼期間長ａrefとの関係（即ち、基準関係）」がａref＝ｆ(θdQpeak）なる関数で表されるとする。
このとき、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭ０であり、且つ、
吸気弁位相角ＩＮＶＴが所定位相角ＩＮＶＴｓであり且つ熱発生率最大時期θdQpeakが任意値θ１であるときの前半燃焼期間長ａｓ（θ１）がａd1であるとすると、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭ０であり、且つ、
吸気弁位相角ＩＮＶＴが所定位相角ＩＮＶＴｓであり且つ熱発生率最大時期θdQpeakが任意値θ２であるときの前半燃焼期間長ａｓ（θ２）は、次の（１Ｂ）式に求めることができる。
ａｓ（θ２）＝ｆ(θ２）＋Δａ１
＝ｆ(θ２）＋（ａd1−ａr1）
＝ｆ(θ２）−（ａr1−ａd1）
＝ｆ(θ２）＋{ａd1−ｆ（θ１）}    …（１Ｂ）
【０１２６】
以下、上記知見１の正当性について説明する。
図２９は、機関回転速度ＮＥが一定速度（＝１２００［ｒｐｍ］）である場合において、吸気弁作用角ＶＣＡＭを一定の作用角（＝２６０［ＣＡ］）に固定し、吸気弁位相角ＩＮＶＴを種々の値に変化させた実験において「実測された熱発生率最大時期θdQpeak及び実測された前半燃焼期間長ａ」の関係をプロットしたグラフである。
【０１２７】
図２９に示したグラフから、熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を規定する曲線（破線）は、吸気弁位相角ＩＮＶＴと強い相関があり、且つ、吸気弁位相角ＩＮＶＴが変化すると前半燃焼期間長ａが増減する方向（紙面上下方向、即ち、Ｙ軸方向）に平行移動していることが理解される。
【０１２８】
図３０に示されたグラフのそれぞれは、機関回転速度ＮＥが１２００［ｒｐｍ］である場合において、吸気弁作用角ＶＣＡＭを種々の値に固定し、吸気弁位相角ＩＮＶＴを種々の値に変化させることによって吸気弁閉弁時期ＩＶＣを種々の値に変化させた実験において「実測された熱発生率最大時期θdQpeak及び実測された前半燃焼期間長ａ」の関係をプロットしたグラフである。
【０１２９】
なお、図３０に示された各グラフにおける凡例を表１に示す。
【表１】
2016160761.tif 000012
【０１３０】
図３０に示されたグラフは、機関回転速度ＮＥが１２００［ｒｐｍ］である場合において吸気弁作用角ＶＣＡＭが種々の値に固定されている場合、吸気弁位相角ＩＮＶＴが変化すると前半燃焼期間長ａが増減する方向に平行移動していることを示している。
【０１３１】
同様に、図３１に示されたグラフのそれぞれは、機関回転速度ＮＥが２４００［ｒｐｍ］である場合において、吸気弁作用角ＶＣＡＭを種々の値に固定し、吸気弁位相角ＩＮＶＴを種々の値に変化させることによって吸気弁閉弁時期ＩＶＣを種々の値に変化させた実験において「実測された熱発生率最大時期θdQpeak及び実測された前半燃焼期間長ａ」の関係をプロットしたグラフである。なお、図３１に示された各グラフにおける凡例も表１に示した凡例と同じである。
【０１３２】
図３１に示されたグラフは、機関回転速度ＮＥが２４００［ｒｐｍ］である場合において吸気弁作用角ＶＣＡＭが種々の値に固定されている場合、吸気弁位相角ＩＮＶＴが変化すると前半燃焼期間長ａが増減する方向に平行移動していることを示している。更に、図３０に示されたグラフと図３１に示されたグラフとから、機関回転速度ＮＥが相違しても、機関回転速度ＮＥが一定に維持されおり、且つ、吸気弁作用角ＶＣＡＭが一定に維持されている限り、吸気弁位相角ＩＮＶＴが変化すると前半燃焼期間長ａが増減する方向に平行移動することが理解される。
【０１３３】
以上から、図３０及び図３１に示されたグラフによっても、上記知見１の正当性が確認される。
【０１３４】
更に、図３及び図４を参照しながら説明したように、本発明者は以下の知見２も得た。
（知見２）
機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁閉弁時期ＩＶＣが任意閉弁時期ＩＶＣ０であり且つ吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒである場合における「熱発生率最大時期θdQpeakと前半燃焼期間長ａrefとの関係（即ち、基準関係）」がａref＝ｇ(θdQpeak）なる関数で表されるとする。
このとき、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁閉弁時期ＩＶＣが任意閉弁時期ＩＶＣ０であり、且つ、
吸気弁作用角ＶＣＡＭが所定作用角ＶＣＡＭｓであり且つ熱発生率最大時期θdQpeakが任意値（任意の基準値）θ３であるときの前半燃焼期間長ａｓ（θ３）がａd3であるとすると、
・機関回転速度ＮＥが任意速度ＮＥ０であり且つ吸気弁閉弁時期ＩＶＣが任意閉弁時期ＩＶＣ０であり、且つ、
吸気弁作用角ＶＣＡＭが所定作用角ＶＣＡＭｓであり且つ熱発生率最大時期θdQpeakが任意値θ４であるときの前半燃焼期間長ａｓ（θ４）は、次の（２Ｂ）式により求めることができる。
ａｓ(θ４）＝ｇ(θ４）＋Δａ３
＝ｇ(θ４）＋（ａd3−ａr3）
＝ｇ(θ４）−（ａr3−ａd3）
＝ｇ(θ４）＋｛ａd3−ｇ（θ３）｝    …（２Ｂ）
【０１３５】
以下、上記知見２の正当性について説明する。
図３２は、機関回転速度ＮＥが一定速度（＝２４００［ｒｐｍ］）である場合において、吸気弁閉弁時期ＩＶＣを所定時期（＝３０［ABDC deg］）に固定し、吸気弁作用角ＶＣＡＭを種々の値に変化させた実験において「実測された熱発生率最大時期θdQpeak及び実測された前半燃焼期間長ａ」の関係をプロットしたグラフである。
【０１３６】
図３２に示したグラフから、熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係を規定する曲線（破線）は、吸気弁作用角ＶＣＡＭと強い相関があり、且つ、吸気弁作用角ＶＣＡＭが変化すると前半燃焼期間長ａが増減する方向（紙面上下方向、即ち、Ｙ軸方向）に平行移動していることが理解される。
【０１３７】
図３３に示されたグラフのそれぞれは、機関回転速度ＮＥが１２００［ｒｐｍ］である場合において、吸気弁閉弁時期ＩＶＣを種々の値に固定し、吸気弁作用角ＶＣＡＭを種々の値に変化させた実験において「実測された熱発生率最大時期θdQpeak及び実測された前半燃焼期間長ａ」の関係をプロットしたグラフである。
【０１３８】
なお、図３３に示された各グラフにおける凡例を表２に示す。
【表２】
2016160761.tif 000013
【０１３９】
図３３に示されたグラフは、機関回転速度ＮＥが１２００［ｒｐｍ］である場合において吸気弁閉弁時期ＩＶＣが種々の値に固定されている場合、吸気弁作用角ＶＣＡＭが変化すると前半燃焼期間長ａが増減する方向に平行移動していることを示している。
【０１４０】
同様に、図３４に示されたグラフのそれぞれは、機関回転速度ＮＥが２４００［ｒｐｍ］である場合において、吸気弁閉弁時期ＩＶＣを種々の値に固定し、吸気弁作用角ＶＣＡＭを種々の値に変化させた実験において「実測された熱発生率最大時期θdQpeak及び実測された前半燃焼期間長ａ」の関係をプロットしたグラフである。なお、図３４に示された各グラフにおける凡例も表２に示した凡例と同じである。
【０１４１】
図３４に示されたグラフは、機関回転速度ＮＥが２４００［ｒｐｍ］である場合において吸気弁閉弁時期ＩＶＣが種々の値に固定されている場合、吸気弁作用角ＶＣＡＭが変化すると前半燃焼期間長ａが増減する方向に平行移動していることを示している。更に、図３３に示されたグラフと図３４に示されたグラフとから、機関回転速度ＮＥが相違しても、機関回転速度ＮＥが一定に維持されおり、且つ、吸気弁閉弁時期ＩＶＣが一定に維持されている限り、吸気弁作用角ＶＣＡＭが変化すると前半燃焼期間長ａが増減する方向に平行移動することが理解される。
【０１４２】
以上から、図３３及び図３４に示されたグラフによっても、上記知見２の正当性が確認される。従って、知見１及び知見２から上述したように導き出せる知見３の正当性も確認された。
【０１４３】
以上、説明したように、本実施装置１０は本実施方法を用いて、前半燃焼期間長ａを含む熱発生率波形を推定するために必要なパラメータ（着火遅れ期間長τ、着火時期ＦＡ、前半燃焼期間長ａ、熱発生率最大時期θdQpeak、熱発生率傾きｂ／ａ及び熱発生量Ｑ１、Ｑ２及びＱall、後半燃焼期間長ｃ等）を精度良く推定することができる。従って、内燃機関の開発及び設計等に極めて有用な情報を提供することができる。
【０１４４】
特に、本実施方法では、上記知見１及び知見２に基づいて、吸気弁の開弁特性が「ある基準特性」と異なる特性（特定特性）となった場合であっても、前半燃焼期間長ａを容易に推定することができる。換言すれば、吸気弁の開弁特性毎に前半燃焼期間長ａの予測式（上記（３）式）を独立に作成する必要がない。
【０１４５】
即ち、吸気弁開弁特性が基準特性である場合における「前半燃焼期間長ａと熱発生率最大時期θdQpeakとの関係（基準関係）」を実測するか、或いは、上記（３）式又はＷｉｅｂｅ関数等に基づいて推定して取得するとともに、吸気弁開弁特性が特定特性である場合における「ある熱発生率最大時期θdQpeakにおける前半燃焼期間長ａ」を一点だけ取得することにより（又は、前半燃焼期間長ａを２〜３の互いに異なる熱発生率最大時期θdQpeakに対して取得し、それらの平均を求めることにより）、その異なる特定特性におけるあらゆる熱発生率最大時期に対する前半燃焼期間長ａを上記（１Ｂ）式及び上記（２Ｂ）式を用いて簡単に推定することができる。
【０１４６】
更に、具体的に述べると、本実施方法においては、上記（３）式は「上述した関数ｆ(θdQpeak）又は関数ｇ(θdQpeak）」として採用され得る。
【０１４７】
ここで、上記（３）式は、上述したように、負荷率ＫＬ、ＥＧＲ率Ｇegr、空燃比Ａ／Ｆ及び冷却水温ＴＨＷ（機関暖機状態）に関わらず成立するから、これらの運転状態パラメータ毎に実験を繰り返して（３）式（（３）式に用いられるＣ３、α及びβ）を独立に決定する必要がない。従って、これらの運転状態パラメータの組み合わせに対して前半燃焼期間長ａを求めるのに必要な時間及び労力を削減できる。なお、他のモデル式においても同様であるが、Ｃ３、α及びβ等の各モデル式で使用される値は、代表的なデータ点に基づいて決定することができる。
【０１４８】
更に、上記（３）式が作成されれば、吸気弁開弁特性が基準特性と異なる特定特性にある場合において熱発生率最大時期θdQpeakが任意の基準値（上記θ１又は上記θ３）であるときの前半燃焼期間長ａのみを実測するだけで、吸気弁開弁特性が特定特性にある場合の「任意の熱発生率最大時期θdQpeakに対する前半燃焼期間長ａ」を推定することができる。その結果、前半燃焼期間長ａを求めるのに必要な時間及び労力を大幅に削減することができる。
【０１４９】
（適用例１）
更に、例えば、上記（３）式及び本実施方法の知見１（上記（１Ａ）式又は（１Ｂ）式）は、図３５にフローチャートにより示したルーチンにおいても用いられる。このルーチンは、吸気弁位相角ＩＮＶＴが特定位相角ＩＮＶＴｘであり且つ吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒである場合（即ち、吸気弁の開弁特性が、基準特性に対して、吸気弁位相角ＩＮＶＴのみが相違する特定特性である場合）、図１５に示したルーチンに代わって用いられる。以下、このルーチンについて簡単に説明する。なお、図３５に示したステップであって図１５に示されているステップには同じ符号を付し、それらのステップについての説明は省略される。更に、ＣＰＵはステップ１５０５にて特定位相角ＩＮＶＴｘを入力する。
【０１５０】
ＣＰＵは、ステップ１５１５にて仮想熱発生率最大時期θzdQpeakにおける燃焼室容積Ｖ@θzdQpeakを算出すると、以下に述べるステップ３５１０乃至ステップ３５４０の処理を順に行い、ステップ１５２５に進む。
【０１５１】
ステップ３５１０：ＣＰＵは、上記（３）式と同じ式に燃焼室容積Ｖ@θzdQpeak及び機関回転速度ＮＥを代入することによって、仮想基準前半燃焼期間長ａza［ＣＡ］を算出する。仮想基準前半燃焼期間長ａzaは、吸気弁開弁特性が基準特性に設定されている場合に対して作成された上記（３）式により求められる前半燃焼期間長の基準値である。
【０１５２】
ステップ３５２０：ＣＰＵは、基準作用角ＶＣＡＭｒと機関回転速度ＮＥとの組み合わせに基づいてルックアップテーブルMapΔａ１(ＩＮＶＴ）を選択する。テーブルMapΔａ１(ＩＮＶＴ）は、「基準作用角ＶＣＡＭｒ及び機関回転速度ＮＥの組み合わせ」毎にメモリに格納されている。
【０１５３】
テーブルMapΔａ１(ＩＮＶＴ）は、上記知見１に基づいて予め作成される。即ち、テーブルMapΔａ１(ＩＮＶＴ）は、機関回転速度ＮＥの下で吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒに維持されている状態において吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒから所定位相角ＩＮＶＴｘに変更された場合に対して、熱発生率最大時期θdQpeakが任意の基準値θ１であるときの前半燃焼期間長ａの増大量（上記（１Ａ）式における（ａd1−ａr1）に相当する値であり、負の値を含む）Δａ１を予め実測及び計算等により求め、その増大量Δａ１と所定位相角ＩＮＶＴｘとの関係を「基準作用角ＶＣＡＭｒ及び機関回転速度ＮＥの組み合わせ」に関連付けて記憶したテーブルである。
【０１５４】
ステップ３５３０：ＣＰＵは、選択されたテーブルMapΔａ１(ＩＮＶＴ）に吸気弁位相角ＩＮＶＴｘを適用することにより（即ち、変数（引数）ＩＮＶＴとして吸気弁位相角ＩＮＶＴｘを用いることにより）補正量Δａ１を取得する。
【０１５５】
ステップ３５４０：ＣＰＵは、取得された補正量Δａ１を仮想基準前半燃焼期間長ａzaに加えることにより、仮想前半燃焼期間長ａｚを推定する。
【０１５６】
このように、知見１に基づく本実施方法によれば、吸気弁位相角ＩＮＶＴが変化した場合であっても、簡単に仮想前半燃焼期間長ａｚを推定することができる。更に、テーブルMapΔａ１(ＩＮＶＴ）は、吸気弁位相角ＩＮＶＴが種々の値に変化した場合に上記（３）式を独立して作成する場合に比べ、極めて少ない時間及び労力で作成することができる。
【０１５７】
（適用例２）
更に、例えば、上記（３）式及び本実施方法の知見２（上記（２Ａ）式又は（２Ｂ）式）は、図３６にフローチャートにより示したルーチンにおいても用いられる。このルーチンは、吸気弁作用角ＶＣＡＭが特定作用角ＶＣＡＭｘであり且つ吸気弁閉弁時期ＩＶＣが基準閉弁時期ＩＶＣｒである場合（即ち、吸気弁の開弁特性が、基準特性に対して、吸気弁作用角ＶＣＡＭのみが相違する特定特性である場合）、図１５に示したルーチンに代わって用いられる。以下、このルーチンについて簡単に説明する。なお、図３６に示したステップであって図１５に示されているステップには同じ符号を付し、それらのステップについての説明は省略される。更に、ＣＰＵはステップ１５０５にて特定作用角ＶＣＡＭｘを入力する。
【０１５８】
ＣＰＵは、ステップ１５１５にて仮想熱発生率最大時期θzdQpeakにおける燃焼室容積Ｖ@θzdQpeakを算出すると、以下に述べるステップ３６１０乃至ステップ３６４０の処理を順に行い、ステップ１５２５に進む。
【０１５９】
ステップ３６１０：ＣＰＵは、上記（３）式と同じ式に燃焼室容積Ｖ@θzdQpeak及び機関回転速度ＮＥを代入することによって、仮想基準前半燃焼期間長ａza［ＣＡ］を算出する。仮想基準前半燃焼期間長ａzaは、「吸気弁開弁特性が基準特性に設定されている場合（吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒに設定され且つ吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒに設定されている場合）に対して作成された上記（３）式」により求められる前半燃焼期間長の基準値である。
【０１６０】
ステップ３６２０：ＣＰＵは、基準閉弁時期ＩＶＣｒと機関回転速度ＮＥとの組み合わせに基づいてルックアップテーブルMapΔａ3(ＶＣＡＭ）を選択する。テーブルMapΔａ3(ＶＣＡＭ）は、「基準閉弁時期ＩＶＣｒ及び機関回転速度ＮＥの組み合わせ」毎にメモリに格納されている。
【０１６１】
テーブルMapΔａ3(ＶＣＡＭ）は、上記知見２に基づいて予め作成される。即ち、テーブルMapΔａ3(ＶＣＡＭ）は、機関回転速度ＮＥの下で吸気弁閉弁時期ＩＶＣが基準閉弁時期ＩＶＣｒに維持されている状態において吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒから所定作用角ＶＣＡＭｘに変更された場合に対して、熱発生率最大時期θdQpeakが任意の基準値θ３であるときの前半燃焼期間長ａの増大量（上記（２Ａ）式における（ａd3−ａr3）に相当する値であり、負の値を含む）Δａ３を予め実測及び計算等により求め、その増大量Δａ３と所定作用角ＶＣＡＭｘとの関係を「基準閉弁時期ＩＶＣｒ及び機関回転速度ＮＥの組み合わせ」に関連付けて記憶したテーブルである。
【０１６２】
ステップ３６３０：ＣＰＵは、選択されたテーブルMapΔａ3(ＶＣＡＭ）に作用角ＶＣＡＭｘを適用することにより（即ち、変数（引数）ＶＣＡＭとして作用角ＶＣＡＭｘを用いることにより）補正量Δａ３を取得する。
【０１６３】
ステップ３６４０：ＣＰＵは、取得された補正量Δａ３を仮想基準前半燃焼期間長ａzaに加えることにより、仮想前半燃焼期間長ａｚを推定する。
【０１６４】
このように、知見２に基づく本実施方法によれば、吸気弁作用角ＶＣＡＭが変化した場合であっても、簡単に仮想前半燃焼期間長ａｚを推定することができる。更に、テーブルMapΔａ3(ＶＣＡＭ）は、吸気弁作用角ＶＣＡＭが種々の値に変化した場合に上記（３）式を独立して作成する場合に比べ、極めて少ない時間及び労力で作成することができる。
【０１６５】
以上、説明したように、本実施装置１０及び本実施方法においては、燃焼状態を表す重要なパラメータの一つである前半燃焼期間長ａが吸気弁の開弁特性に応じて精度良く且つ少ない労力にて推定される。従って、内燃機関の開発工数を低減することができる。更に、仮に、燃焼状態を表す他のパラメータ（着火遅れ期間長τ及び熱発生率傾きｂ／ａ等）が吸気弁開弁特性の影響を受けて真値と相違していても、少なくとも前半燃焼期間長ａは吸気弁開弁特性を精度良く反映した値となるから、全体として熱発生率波形の精度も向上することができる。
【０１６６】
なお、本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記知見１及び知見２は、組み合わせて用いることができる。つまり、吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒであり且つ吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒである場合における「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係」を基準関数ｆ１として求めておき、先ず、吸気弁位相角ＩＮＶＴのみが位相角ＩＮＶＴｘへと変化した場合の「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係」を表す関数ｆ２を上記知見１に基づいて求める。そして、その関数ｆ２を新たな基準関数とし（即ち、位相角ＩＮＶＴｘを基準位相角と考え）、吸気弁作用角ＶＣＡＭのみが作用角ＶＣＡＭｘへと変化した場合の「熱発生率最大時期θdQpeakと前半燃焼期間長ａとの関係」を表す関数ｆ３を上記知見２に基づいて求める。
【０１６７】
或いは、上記ステップ３５４０又は上記３６４０にて仮想基準前半燃焼期間長ａzaに加えられる補正量Δａを、テーブルMapΔａ(ＮＥ，ＩＮＶＴ，ＶＣＡＭ）として求めておき、このテーブルMapΔａ(ＮＥ，ＩＮＶＴ，ＶＣＡＭ）に、実際の機関回転速度ＮＥ、吸気弁位相角ＩＮＶＴｘ及び吸気弁作用角ＶＣＡＭｘを適用して補正量Δａを取得し、その補正量Δａ（Δａ１又はΔａ３）を仮想基準前半燃焼期間長ａzaに加えることにより、吸気弁位相角ＩＮＶＴ及び吸気弁作用角ＶＣＡＭが基準特性のそれら（即ち、基準位相角ＩＮＶＴｒ及び基準作用角ＶＣＡＭｒのそれぞれ）から変化したときの前半燃焼期間長ａ(仮想前半燃焼期間長ａｚ）を推定してもよい。
【０１６８】
このテーブルMapΔａ(ＮＥ，ＩＮＶＴ，ＶＣＡＭ）は、上記知見１及び２に基づいて予め作成される。即ち、テーブルMapΔａ(ＮＥ，ＩＮＶＴ，ＶＣＡＭ）は、機関回転速度ＮＥの下で吸気弁位相角ＩＮＶＴが基準位相角ＩＮＶＴｒに設定され且つ吸気弁作用角ＶＣＡＭが基準作用角ＶＣＡＭｒに設定されている場合から、機関回転速度ＮＥの下で吸気弁位相角ＩＮＶＴが任意位相角ＩＮＶＴｘに変更され且つ吸気弁作用角ＶＣＡＭが任意作用角ＶＣＡＭｘに変更された場合の「熱発生率最大時期θdQpeakが任意の基準値θ１であるときの前半燃焼期間長ａの増大量（補正量、差分値）Δａ」を予め実測及び計算等により求め、その増大量Δａを、機関回転速度ＮＥ、任意位相角ＩＮＶＴｘ及び任意作用角ＶＣＡＭｘに関連付けて記憶したテーブルである。
【０１６９】
従って、この補正量（差分値）Δａを求めるステップ（ステップ３５３０又はステップ３６３０に対応するステップ）は、
機関回転速度が所定の基準回転速度であり且つ機関の吸気弁の開弁特性が所定の基準開弁特性に設定された基準状態において熱発生率最大時期が所定の第１時期（任意の基準値θ１）であるときの前半燃焼期間長を基準関係（基準関数ｆ１）に基づいて第１基準期間長として取得するステップ（第１ステップ）；
機関回転速度が前記基準回転速度であり且つ前記吸気弁の開弁特性が「前記基準開弁特性と異なる特定開弁特性」に設定された特定状態において前記熱発生率最大時期が前記第１時期であるときの前記前半燃焼期間長を第１特定期間長として実測により取得するステップ（第２ステップ）；及び
前記第１特定期間長から前記第１基準期間長を減じることにより第１差分値（例えば、Δａ１＝ａd1−ａr1、又は、Δａ３＝ａd3−ａr3）を取得するか、又は、前記第１基準期間長から前記第１特定期間長を減じることにより第２差分値（例えば、−Δａ１＝ａr1−ａd1、又は、−Δａ３＝ａr3−ａd3）を取得するステップ（第３ステップ）；
の各ステップを実行することに相当する。
なお、第１ステップでは、第１基準期間長を実測に求めても良い。
【０１７０】
また、仮想基準前半燃焼期間長ａzaを求めるステップ（ステップ３５１０又はステップ３６１０に対応するステップ）は、前記基準状態において熱発生率最大時期が前記第１時期と異なる第２時期（仮想熱発生率最大時期θzdQpeak）であるときの前半燃焼期間長を前記基準関係に基づいて第２基準期間長として取得するステップ（第４ステップ）に相当し、
補正量Δａ（Δａ１又はΔａ３）を仮想基準前半燃焼期間長ａzaに加えるステップ（ステップ３５４０又はステップ３６４０に対応するステップ）は、前記特定状態において前記熱発生率最大時期が前記第２時期であるときの前半燃焼期間長（仮想前半燃焼期間長ａｚ）を、前記第２基準期間長（仮想基準前半燃焼期間長ａza）に前記第１差分値（Δａ１、Δａ３）を加えることにより推定するステップ（第５ステップ）、又は、前記第２基準期間長（仮想基準前半燃焼期間長ａza）から前記第２差分値（−Δａ１、−Δａ３）を減じることにより推定するステップに相当する、と言うことができる。
【０１７１】
更に、各モデル式の前記前提条件が成立しない場合、その成立していない条件に対する補正係数を別途求めておき、その補正係数によって「着火遅れ期間の長さτ、前半燃焼期間長ａ及び熱発生率傾きｂ／ａ」のそれぞれを補正してもよい。
【符号の説明】
【０１７２】
１０…シミュレーション装置、１１…着火遅れ期間推定部（着火遅れ期間モデル）、１２…前半燃焼期間推定部（前半燃焼期間推定モデル）、１３…熱発生率傾き推定部（熱発生率傾きモデル）、１４…熱発生量推定部（熱発生量モデル）。
【図１】
2016160761.tif 000014
【図２】
2016160761.tif 000015
【図３】
2016160761.tif 000016
【図４】
2016160761.tif 000017
【図５】
2016160761.tif 000018
【図６】
2016160761.tif 000019
【図７】
2016160761.tif 000020
【図８】
2016160761.tif 000021
【図９】
2016160761.tif 000022
【図１０】
2016160761.tif 000023
【図１１】
2016160761.tif 000024
【図１２】
2016160761.tif 000025
【図１３】
2016160761.tif 000026
【図１４】
2016160761.tif 000027
【図１５】
2016160761.tif 000028
【図１６】
2016160761.tif 000029
【図１７】
2016160761.tif 000030
【図１８】
2016160761.tif 000031
【図１９】
2016160761.tif 000032
【図２０】
2016160761.tif 000033
【図２１】
2016160761.tif 000034
【図２２】
2016160761.tif 000035
【図２３】
2016160761.tif 000036
【図２４】
2016160761.tif 000037
【図２５】
2016160761.tif 000038
【図２６】
2016160761.tif 000039
【図２７】
2016160761.tif 000040
【図２８】
2016160761.tif 000041
【図２９】
2016160761.tif 000042
【図３０】
2016160761.tif 000043
【図３１】
2016160761.tif 000044
【図３２】
2016160761.tif 000045
【図３３】
2016160761.tif 000046
【図３４】
2016160761.tif 000047
【図３５】
2016160761.tif 000048
【図３６】
2016160761.tif 000049