JP2005226481A - エンジンのデポジット量検出装置及びノック制御装置 - Google Patents

Abstract

【課題】燃焼室内のデポジット量を直接的に算出し得る装置を提供する。
【解決手段】熱発生率を算出する熱発生率算出手段（３１）と、ノック発生時期を検出または推定するノック発生時期検出・推定手段（３１）と、実ノック強度を検出する実ノック強度検出手段（４７）と、これら熱発生率、ノック発生時期及び実ノック強度に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段（３１）とを備える。
【選択図】図１

Description

本発明は、エンジン（内燃機関）のデポジット量検出装置、特に燃焼室内のデポジット量を検出（推定）するもの及びその検出したデポジット量を用いたノック制御装置に関する。

可変バルブタイミング機構を備えたエンジンにおいて、加速時の空燃比に基づいて吸気系におけるデポジット量を算出し、デポジット量に応じてバルブタイミングを遅角させてノッキングを回避するものがある（特許文献１参照）。
特開平９−３０３１６５号公報

ところで、上記の特許文献１の技術は、吸気系に堆積するデポジットが増えれば、燃焼室内に堆積するデポジットも増えるであろうとみなすものであり、燃焼室内のデポジット量を直接的に算出してはいない。

このため、吸気系にデポジットが堆積していても燃焼室内にデポジットが堆積していない場合や吸気系にデポジットが堆積していないのに燃焼室内にはデポジットが堆積している場合に燃焼室内のデポジット量の算出に誤差が生じてしまう。

そこで本発明は、燃焼室内のデポジット量を直接的に算出し得る装置及びその算出した燃焼室内のデポジット量を用いたノック制御装置を提供することを目的とする。

本発明は、熱発生率と燃焼質量割合（ＢＲ）とのいずれかを算出し、ノック発生時期（θｋｎｋ）を検出または推定し、実ノック強度（ＫＩＤ）を検出し、これら熱発生率と燃焼質量割合（ＢＲ）とのいずれか、ノック発生時期（θｋｎｋ）及び実ノック強度（ＫＩＤ）に基づいて燃焼室内のデポジット量（ＶＤＥＰ）を算出するように構成する。

また、本発明は、こうした構成にさらに前記算出したデポジット量に基づいてノックが生じないように点火時期を遅角側に制御する。この場合に、吸気弁開閉時期を可変に調整し得る可変バルブタイミング機構を備えるエンジンでは、前記算出したデポジット量に基づいてノックが生じないように吸気弁開時期を遅角側に制御する。

ノック発生時期に基づいてノック発生時期の燃焼室容積（Ｖｋｎｋ）を算出することができ、ノック発生時期の熱発生率やノック発生時期の燃焼質量割合（ＢＲｋｎｋ）をこのノック発生時期の燃焼室容積（Ｖｋｎｋ）で割った値はノック強度と相関がある。この値を用いてノック強度推定値を算出し、このノック強度推定値と実ノック強度を比べることで燃焼室内のデポジット量を算出できる。本発明はこうした考え方に基づき、熱発生率やこの微分値に相当する燃焼質量割合を算出し、ノック発生時期（θｋｎｋ）を検出または推定し、実ノック強度（ＫＩＤ）を検出し、これら熱発生率と燃焼質量割合のいずれか、ノック発生時期（θｋｎｋ）及び実ノック強度（ＫＩＤ）に基づいて燃焼室内のデポジット量（ＶＤＥＰ）を算出するので、吸気系にデポジットが堆積していても燃焼室内にデポジットが堆積していない場合や吸気系にデポジットが堆積していないのに燃焼室内にはデポジットが堆積している場合においても、燃焼室内のデポジット量を精度よく算出できる。

また、デポジット量の算出結果を情報として知らせるようにすることで、サービス工場における清掃等の対応が可能となり、エンジンの効率低下や排気悪化等を事前に防ぐことができる。

一方、燃焼室内にデポジットが堆積すると圧縮比が上がりノッキングが発生しやすくなる。この場合に本発明によれば算出したデポジット量に基づいてノックが生じないように点火時期や吸気弁開時期を遅角側に制御するので、燃焼室内のデポジットの堆積に伴うノッキングの発生を防止できる。

以下、図面に基づき本発明の実施形態について説明する。

図１は、本発明のシステムを説明するための概略図である。

空気は吸気コレクタ２に蓄えられた後、吸気マニホールド３を介して各気筒の燃焼室５に導入される。燃料は各気筒の吸気ポート４に配置された燃料インジェクタ２１より噴射供給される。空気中に噴射された燃料は気化しつつ空気と混合してガス（混合気）を作り、燃焼室５に流入する。この混合気は吸気弁１５が閉じることで燃焼室５内に閉じこめられ、ピストン６の上昇によって圧縮される。

この圧縮混合気に対して高圧火花により点火を行うため、パワートランジスタ内蔵の点火コイルを各気筒に配した電子配電システムの点火装置１１を備える。すなわち、点火装置１１は、バッテリからの電気エネルギーを蓄える点火コイル１３と、点火コイル１３の一次側への通電、遮断を行うパワートランジスタと、燃焼室５の天井に設けられ点火コイル１３の一次電流の遮断によって点火コイル１３の二次側に発生する高電圧を受けて、火花放電を行う点火プラグ１４とからなっている。

圧縮上死点より少し手前で点火プラグ１４により火花が飛ばされ圧縮混合気に着火されると、火炎が広がりやがて爆発的に燃焼し、この燃焼によるガス圧がピストン６を押し下げる仕事を行う。この仕事はクランクシャフト７の回転力として取り出される。燃焼後のガス（排気）は排気弁１６が開いたとき排気通路８へと排出される。

排気通路８には三元触媒９を備える。三元触媒９は排気の空燃比が理論空燃比を中心とした狭い範囲（ウインドウ）にあるとき、排気に含まれるＨＣ、ＣＯ、ＮＯｘといった有害三成分を同時に効率よく除去できる。空燃比は吸入空気量と燃料量の比であるので、エンジンの１サイクル（４サイクルエンジンではクランク角で７２０°区間）当たりに燃焼室５に導入される吸入空気量と、燃料インジェクタ２１からの燃料噴射量との比が理論空燃比となるように、エンジンコントローラ３１ではエアフローメータ３２からの吸入空気流量の信号とクランク角センサ（３３、３４）からの信号に基づいて燃料インジェクタ２１からの燃料噴射量を定めると共に、三元触媒９の上流に設けたＯ₂センサ３５からの信号に基づいて空燃比をフィードバック制御している。

吸気コレクタ２の上流には絞り弁２３がスロットルモータ２４により駆動される、いわゆる電子制御スロットル２２を備える。運転者が要求するトルクはアクセルペダル４１の踏み込み量（アクセル開度）に現れるので、エンジンコントローラ３１ではアクセルセンサ４２からの信号に基づいて目標トルクを定め、この目標トルクを実現するための目標空気量を定め、この目標空気量が得られるようにスロットルモータ２４を介して絞り弁２３の開度を制御する。

吸気弁用カムシャフト２５、排気弁用カムシャフト２６及びクランクシャフト７の各前部にはそれぞれカムスプロケット、クランクスプロケットが取り付けられ、これらスプロケットにタイミングチェーン（図示しない）を掛け回すことで、カムシャフト２５、２６がエンジンのクランクシャフト７により駆動されるのであるが、このカムスプロケットと吸気弁用カムシャフト２５との間に介在して、作動角一定のまま吸気弁用カムの位相を連続的に制御し得る可変吸気バルブタイミングコントロール機構（以下、「吸気ＶＴＣ機構」という。）２７と、カムスプロケットと排気弁用カムシャフト２６との間に介在して、作動角一定のまま排気弁用カムの位相を連続的に制御し得る可変排気バルブタイミングコントロール機構（以下、「排気ＶＴＣ機構」という。）２８とを備える。吸気弁１５の開閉時期や排気弁１６の開閉時期を変えると燃焼室５に残留する不活性ガスの量が変化する。燃焼室５内の不活性ガスの量が増えるほどポンピングロスが減って燃費がよくなるので、運転条件によりどのくらいの不活性ガスが燃焼室５内に残留したらよいかを目標吸気弁閉時期や目標排気弁閉時期にして予め定めており、エンジンコントローラ３１ではそのときの運転条件（エンジンの負荷と回転速度）より目標吸気弁閉時期と目標排気弁閉時期を定め、それら目標値が得られるように吸気ＶＴＣ機構２７、排気ＶＴＣ機構２８の各アクチュエータを介して吸気弁閉時期と排気弁閉時期を制御する。

吸気温度センサ４３からの吸気温度の信号、吸気圧力センサ４４からの吸気圧力の信号、排気温度センサ４５からの排気温度の信号、排気圧力センサ４６からの排気圧力の信号が、水温センサ３７からの冷却水温の信号と共に入力されるエンジンコントローラ３１では、パワートランジスタ１３を介して点火プラグ１４の一次側電流の遮断時期である点火時期を制御する。

図２はエンジンコントローラ３１内で行われる点火時期制御のブロック図で、大きくは点火時期演算部５１と点火時期制御部６１とからなる。点火時期演算部５１はさらに初期燃焼期間算出部５２、主燃焼期間算出部５３、燃焼期間算出部５４、基本点火時期算出部５５、前回燃焼開始時期算出部５６及び点火時期指令値算出部５７からなる。

初期燃焼期間算出部５２では、混合気が着火してから火炎核が形成されるまでの期間を初期燃焼期間ＢＵＲＮ１として算出する。主燃焼期間算出部５３では、火炎核が形成されてから燃焼圧力が最大値Ｐｍａｘに達するまでの期間を主燃焼期間ＢＵＲＮ２として算出する。燃焼期間算出部５４では、これら初期燃焼期間ＢＵＲＮ１と主燃焼期間ＢＵＲＮ２との合計を、燃焼開始より最大燃焼圧力Ｐｍａｘに至るまでの燃焼期間ＢＵＲＮとして算出する。基本点火時期算出部５５では、この燃焼期間ＢＵＲＮに基づいてＭＢＴの得られる点火時期（この点火時期を「基本点火時期」という。）ＭＢＴＣＡＬを算出する。

点火時期指令値算出部５７は図４０に示したように、ノックによる圧力上昇量推定部８０、ノック強度推定値算出部８６、ノック判定部８７、点火時期遅角補正部８８、デポジット量算出部８９からなる。

ノックによる圧力上昇量推定部８０では運転条件に基づいて燃焼室５内のノックによる圧力上昇分ＤＰ（ノックによる圧力上昇量）を推定する。このノックによる圧力上昇量推定部８０は、さらに自着火時期算出部８１、未燃燃料量算出部８２、自着火時容積算出部８３、既燃ガス比熱算出部８４、ノックによる圧力上昇量算出部８５からなる。

自着火時期算出部８１では燃焼室５内の燃料が自着火にいたるまでの時間τの逆数の分布を表す特性に基づいて燃焼室５内の燃料が自着火する時期（以下「自着火時期」という。）θｋｎｋを算出（推定）する。未燃燃料量算出部８２では、燃料量ＱＩＮＪとこの自着火時期θｋｎｋとに基づいて計算式により自着火する燃料量である未燃燃料量ＭＵＢを算出する。自着火時容積算出部８３では自着火時期θｋｎｋに基づいて燃焼室５の自着火時期における容積Ｖｋｎｋを算出する。既燃ガス比熱算出部８４では既燃ガス（自着火する燃料ガス）のエンタルピＥと自着火時期θｋｎｋとに基づいて計算式により既燃ガスの定容比熱Ｃｖを算出する。ノックによる圧力上昇分算出部８５では、これら未燃燃料量ＭＵＢ、自着火時期容積Ｖｋｎｋ、既燃ガスの定容比熱Ｃｖとに基づいてノックによる圧力上昇分ＤＰを計算式により算出する。

ノック強度推定値算出部８６ではこのノックによる圧力上昇分ＤＰに基づいてノック強度推定値ＫＩＣを算出する。ノック判定部８７ではこのノック強度推定値ＫＩＣと所定のスライスレベルＳＬとの比較によりノック強度が大きいか否かを判定する。点火時期遅角補正部８８ではこの判定結果に基づきノック強度が大きい場合に基本点火時期ＭＢＴＣＡＬを遅角側に補正した値を点火時期指令値ＱＡＤＶとして求め、これに対してノック強度が大きくない場合に基本点火時期ＭＢＴＣＡＬをそのまま点火時期指令値ＱＡＤＶとする。

燃焼室内デポジット量算出部８９では実ノック強度ＫＩＤとノック強度推定値ＫＩＣに基づいて燃焼室５内のデポジット量ＶＤＥＰを算出する。

図２に戻り点火時期制御部６１ではこのようにして算出された点火時期指令値ＱＡＤＶで点火プラグ１４が燃焼室５内の混合気に対して着火するように、イグニッションコイル１３への通電角と非通電角を制御する。

上記のように燃焼期間ＢＵＲＮを初期燃焼期間ＢＵＲＮ１と主燃焼期間ＢＵＲＮ２に分けて算出し、燃焼期間ＢＵＲＮに応じて基本点火時期ＭＢＴＣＡＬを求めるようにしたのは、燃焼解析より得られた結果に基づくものである。以下、燃焼解析に基づくこの点火時期制御をさらに説明する。

図３に示すようにＭＢＴ（最大トルクの得られる最小進角値）で混合気に点火した場合に混合気の燃焼圧力が最大値Ｐｍａｘとなるクランク角を基準クランク角θＰＭＡＸ［ｄｅｇＡＴＤＣ］とする。基準クランク角θＰＭＡＸは燃焼方式によらずほぼ一定であり、一般に圧縮上死点後１２〜１５度、最大で圧縮上死点後１０〜２０度の範囲にある。

図４に火花点火エンジンにおける燃焼室内の燃焼解析により得られた燃焼質量割合ＢＲ（燃焼ガス質量割合）の変化を示す。燃焼室に供給された燃料に対する燃焼質量の比率を表す燃焼質量割合ＢＲは、点火時に０％であり、完全燃焼によって１００％に達する。基準クランク角θＰＭＡＸにおける燃焼質量割合は一定で約６０％であることが実験により確かめられている。

燃焼質量割合ＢＲが０％から基準クランク角θＰＭＡＸ相当の約６０％に達するまでの変化代に相当する燃焼期間は、燃焼開始直後で燃焼質量割合にも燃焼圧力にもほとんど変化のない期間である初期燃焼期間と、燃焼質量割合と燃焼圧力が急激に増加する主燃焼期間とに分けられる。初期燃焼期間は、燃焼開始から火炎核が形成されるまでの段階であり、火炎核が形成されるのは燃焼質量割合が０％から２％〜１０％まで変化したときである。この初期燃焼期間中は、燃焼圧力や燃焼温度の上昇速度が小さく、燃焼質量割合の変化に対して初期燃焼期間は長い。初期燃焼期間の長さは燃焼室内の温度や圧力の変化の影響を受けやすい。

一方、主燃焼期間においては、火炎核から外側へと火炎が伝播するのであり、その火炎速度（つまり燃焼速度）が急上昇する。そのため、主燃焼期間の燃焼質量割合の変化は初期燃焼期間の燃焼質量割合の変化に比べて大きい。

エンジンコントローラ３１では、燃焼質量割合が２％に達する（変化する）までを初期燃焼期間ＢＵＲＮ１［ｄｅｇ］とし、初期燃焼期間ＢＵＲＮ１の終了後、基準クランク角θＰＭＡＸに至るまでの区間（燃焼室量割合でいえば２％より約６０％に達するまでの間）を主燃焼期間ＢＵＲＮ２［ｄｅｇ］として区別する。そして、初期燃焼期間ＢＵＲＮ１に主燃焼期間ＢＵＲＮ２を加えた合計である燃焼期間ＢＵＲＮ［ｄｅｇ］を算出し、この燃焼期間ＢＵＲＮから基準クランク角θＰＭＡＸ［ｄｅｇＡＴＤＣ］を差し引き、さらに後述する点火無駄時間相当クランク角ＩＧＮＤＥＡＤ［ｄｅｇ］を加えたクランク角位置を、ＭＢＴの得られる点火時期である基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］として設定する。

火炎核の形成される初期燃焼期間での燃焼室５内の圧力、温度は、点火時の圧力、温度とほぼ等価になるが、これから点火時期を算出しようとしているのに、最初から正確な点火時期を設定することはできない。そこで、図２に示したように前回燃焼開始時期算出部５６で基本点火時期の前回値を前回燃焼開始時期ＭＢＴＣＹＣＬ［ｄｅｇＢＴＤＣ］として算出し、この値を初期燃焼期間算出部５２に対して与えるようにし、初期燃焼期間算出部５２において初期燃焼期間の算出をサイクリックに繰り返すことで、精度の高い結果を時間遅れなしに出すようにしている。

次に、エンジンコントローラ３１で実行される上記の基本点火時期ＭＢＴＣＡＬの算出を以下のフローチャートを参照しながら詳述する。

図５は点火時期の算出に必要な各種の物理量を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。

まずステップ１１では、吸気弁閉時期ＩＶＣ［ｄｅｇＢＴＤＣ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、内部不活性ガス率ＭＲＥＳＦＲ［％］、温度センサ３７により検出される冷却水温ＴＷＫ［Ｋ］、目標当量比ＴＦＢＹＡ、クランク角センサにより検出されるエンジン回転速度ＮＲＰＭ［ｒｐｍ］、点火無駄時間ＤＥＡＤＴＩＭＥ［μｓｅｃ］を読み込む。

ここで、クランク角センサはクランクシャフト７のポジションを検出するポジションセンサ３３と、吸気用カムシャフト２５ポジションを検出するフェーズセンサ３４とからなり、これら２つのセンサ３３、３４からの信号に基づいてエンジン回転速度ＮＲＰＭ［ｒｐｍ］が算出されている。

吸気弁閉時期ＩＶＣは吸気ＶＴＣ機構２７に与える指令値から既知である。あるいはフェーズセンサ３４により実際の吸気弁閉時期を検出してもかまわない。

内部不活性ガス率ＭＲＥＳＦＲは燃焼室内に残留する不活性ガス量を燃焼室内の総ガス量で除した値で、その算出については後述する。点火無駄時間ＤＥＡＤＴＩＭＥは一定値である。

目標当量比ＴＦＢＹＡは図示しない燃料噴射量の算出フローにおいて算出されている。目標当量比ＴＦＢＹＡは無名数であり、理論空燃比を１４．７とすると、次式により表される値である。

ＴＦＢＹＡ＝１４．７／目標空燃比 …（１）
例えば（１）式より目標空燃比が理論空燃比のときＴＦＢＹＡ＝１．０となり、目標空燃比が例えば２２．０といったリーン側の値であるとき、ＴＦＢＹＡは１．０未満の正の値である。

ステップ１２では、燃焼室５の吸気弁閉時期ＩＶＣにおける容積（つまり圧縮開始時期での容積）ＶＩＶＣ［ｍ³］を算出する。燃焼室５の吸気弁閉時期における容積ＶＩＶＣは、ピストン６のストローク位置によって決まる。ピストン６のストローク位置はエンジンのクランク角位置によって決まる。

図６を参照して、エンジンのクランクシャフト７１の回転中心７２がシリンダの中心軸７３からオフセットしている場合を考える。コネクティングロッド７４、コネクティングロッド７４とクランクシャフト７１との結節点７５、コネクティングロッド７４とピストンをつなぐピストンピン７６が図に示す関係にあるとする。このときの、燃焼室５の吸気弁閉時期における容積ＶＩＶＣは次式（２）〜（６）で表すことができる。

ＶＩＶＣ＝ｆ１（θivc）＝Ｖｃ＋（π／４）Ｄ²・Ｈivc …（２）
Ｖｃ＝（π／４）Ｄ²・Ｈx／（ε−１） …（３）
Ｈivc＝｛（ＣＮＤ＋ＳＴ²／２）−（ＣＲoff−ＰＩＳoff）²｝^1/2
−｛（ＳＴ／２）・cos（θivc＋θoff）｝
＋（ＣＮＤ²−Ｘ²）^1/2 …（４）
Ｘ ＝（ＳＴ／２）・sin（θivc＋θoff）−ＣＲoff＋ＰＩＳoff …（５）
θoff＝arcsin｛（ＣＲoff−ＰＩＳoff）／（ＣＮＤ・（ＳＴ／２））｝…（６）
ただし、Ｖｃ：隙間容積［ｍ³］、
ε ：圧縮比、
Ｄ ：シリンダボア径［ｍ］、
ＳＴ ：ピストンの全ストローク［ｍ］、
Ｈivc ：吸気弁閉時期におけるピストンピン７６の
ＴＤＣからの距離［ｍ］、
Ｈx ：ピストンピン７６のＴＤＣからの距離の最大値と最小値の
差［ｍ］、
ＣＮＤ ：コネクティングロッド７４の長さ［ｍ］、
ＣＲoff ：結節点７５のシリンダ中心軸７３からのオフセット距離
［ｍ］、
ＰＩＳoff：クランクシャフト回転中心７２のシリンダ中心軸７３から
のオフセット距離［ｍ］、
θivc ：吸気弁閉時期のクランク角［ｄｅｇＡＴＤＣ］、
θoff ：ピストンピン７６とクランクシャフト回転中心７２とを結
ぶ線がＴＤＣにおいて垂直線となす角度［ｄｅｇ］、
Ｘ ：結節点７５とピストンピン７６との水平距離［ｍ］、
吸気弁閉時期のクランク角θivcは前述のように、エンジンコントローラ３１から吸気ＶＴＣ機構２７への指令信号によって決まるので、既知である。式（２）〜（６）にこのときのクランク角θivc（＝ＩＶＣ）を代入すれば、燃焼室５の吸気弁閉時期における容積ＶＩＶＣを算出することができる。したがって、実用上は燃焼室５の吸気弁閉時期における容積ＶＩＶＣは吸気弁閉時期ＩＶＣをパラメータとするテーブルで設定したものを用いる。吸気ＶＴＣ機構２７を備えないときには定数で与えることができる。

ステップ１３では、燃焼室５の吸気弁閉時期ＩＶＣにおける温度（つまり圧縮開始時期温度）ＴＩＮＩ［Ｋ］を算出する。燃焼室５に流入するガスの温度は、燃焼室５に流入する新気と燃焼室５に残留する不活性ガスとが混じったガスの温度であり、燃焼室５に流入する新気の温度は吸気コレクタ２内の新気温度ＴＣＯＬに等しく、また燃焼室５内に残留する不活性ガスの温度は排気ポート部近傍の排気温度ＴＥＸＨで近似できるので、燃焼室５の吸気弁閉時期ＩＶＣにおける温度ＴＩＮＩは吸気弁閉時期ＩＶＣになったタイミングでの、吸気コレクタ２内の新気温度ＴＣＯＬ、排気温度ＴＥＸＨ、燃焼室５内に残留する不活性ガスの割合である内部不活性ガス率ＭＲＥＳＦＲから次式により求めることができる。

ＴＩＮＩ＝ＴＥＸＨ×ＭＲＥＳＦＲ＋ＴＣＯＬ×（１−ＭＲＥＳＦＲ）…（７）
ステップ１４では燃焼室５の吸気弁閉時期ＩＶＣにおける圧力（つまり圧縮開始時期圧力）ＰＩＮＩ［Ｐａ］を算出する。すなわち、吸気弁閉時期ＩＶＣになったタイミングでのコレクタ内圧力ＰＣＯＬを吸気弁閉時期ＩＶＣにおける圧力ＰＩＮＩとして取り込む。

ステップ１５では、燃焼室５内の混合気の燃えやすさを表す反応確率ＲＰＲＯＢＡ［％］を算出する。反応確率ＲＰＲＯＢＡは無次元の値であり、残留不活性ガス率ＭＲＥＳＦＲ、冷却水温ＴＷＫ［Ｋ］、目標当量比ＴＦＢＹＡの３つのパラメータに依存するので、次式により表すことができる。

ＲＰＲＯＢＡ＝ｆ３（ＭＲＥＳＦＲ、ＴＷＫ、ＴＦＢＹＡ） …（８）
具体的に説明すると、ＭＲＥＳＦＲ、ＴＷＫ、ＴＦＢＹＡの３つのパラメータの組み合わせによって得られる反応確率の最大値を１００％とし、これらのパラメータと反応確率ＲＰＲＯＢＡの関係を実験的に求め、求めた反応確率ＲＰＲＯＢＡをパラメータに応じたテーブルとしてエンジンコントローラ３１のメモリに予め格納しておく。ステップ１４ではパラメータに応じてこのテーブルを検索することにより反応確率ＲＰＲＯＢＡを求める。

具体的には、冷却水温ＴＷＫに応じて図７に示すような特性を有する水温補正係数のテーブルと、同様に設定された内部不活性ガス率補正係数のテーブル（図示しない）と、目標当量比ＴＦＢＹＡに応じて図８に示すような特性を有する当量比補正係数のテーブルを予めメモリに格納しておく。各補正係数の最大値はそれぞれ１．０であり、３種類の補正係数の積に反応確率の最大値１００％を掛け合わせることで、反応確率ＲＰＲＯＢＡを算出する。

各テーブルを説明すると、図７に示す水温補正係数は冷却水温ＴＷＫが高いほど大きく、冷却水温ＴＷＫが８０℃以上では１．０になる。図８に示す当量比補正係数は目標当量比ＴＦＢＹＡが１．０のとき、つまり理論空燃比のときに最大値の１．０となり、目標当量比が１．０より大きくても小さくても当量比補正係数は減少する。内部不活性ガス率補正係数は図示しないが、内部不活性ガス率ＭＲＥＳＦＲがゼロの場合に１．０となる。

ステップ１６では、基準クランク角θＰＭＡＸ［ｄｅｇＡＴＤＣ］を算出する。前述のように基準クランク角θＰＭＡＸはあまり変動しないが、それでもエンジン回転速度ＮＲＰＭの上昇に応じて進角する傾向があるため、基準クランク角θＰＭＡＸはエンジン回転速度ＮＲＰＭの関数として次式で表すことができる。

θＰＭＡＸ＝ｆ４（ＮＲＰＭ） …（９）
具体的にはエンジン回転速度ＮＲＰＭから、エンジンコントローラ３１のメモリに予め格納された図９に示す特性のテーブルを検索することにより基準クランク角θＰＭＡＸを求める。算出を容易にするために、基準クランク角θＰＭＡＸを一定とみなすことも可能である。

最後にステップ１７では、点火無駄時間相当クランク角ＩＧＮＤＥＡＤ［ｄｅｇ］を算出する。点火無駄時間相当クランク角ＩＧＮＤＥＡＤは、エンジンコントローラ３１から点火コイル１３の一次電流を遮断する信号を出力したタイミングから点火プラグ１４が実際に点火するまでのクランク角区間で、次式により表すことができる。

ＩＧＮＤＥＡＤ＝ｆ５（ＤＥＡＤＴＩＭＥ、ＮＲＰＭ） …（１０）
ここでは、点火無駄時間ＤＥＡＤＴＩＭＥを２００μｓｅｃとする。（１０）式は、エンジン回転速度ＮＲＰＭから点火無駄時間ＤＥＡＤＴＩＭＥに相当するクランク角である点火無駄時間相当クランク角ＩＧＮＤＥＡＤを算出するためのものである。

図１０は初期燃焼期間ＢＵＲＮ１［ｄｅｇ］を算出するためのもの、また図１２は主燃焼期間ＢＵＲＮ２［ｄｅｇ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１０、図１２は図５に続けて実行する。
図１０、図１２はどちらを先に実行してもかまわない。

まず図１０から説明すると、ステップ１６１では、前回燃焼開始時期ＭＢＴＣＹＣＬ［ｄｅｇＢＴＤＣ］、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、図５のステップ１３で算出されている燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出されている燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、図５のステップ１５で算出されている反応確率ＲＰＲＯＢＡ［％］を読み込む。

ここで、前回燃焼開始時期ＭＢＴＣＹＣＬは、基本点火時期ＭＢＴＣＡＬの［ｄｅｇＢＴＤＣ］の１サイクル前の値であり、その算出については図１３により後述する。

ステップ１６２では燃焼室５の燃焼開始時期における容積Ｖ０［ｍ³］を算出する。前述したように、ここでの点火時期（燃焼開始時期）は今回のサイクルで演算する基本点火時期ＭＢＴＣＡＬではなく基本点火時期の１サイクル前の値である。すなわち、基本点火時期の１サイクル前の値であるＭＢＴＣＹＣＬから次式により燃焼室５の燃焼開始時期における容積Ｖ０を算出する。

Ｖ０＝ｆ６（ＭＢＴＣＹＣＬ） …（１１）
具体的には前回燃焼開始時期ＭＢＴＣＹＣＬにおけるピストン６のストローク位置と、燃焼室５のボア径から、燃焼室５のＭＢＴＣＹＣＬにおける容積Ｖ０を算出する。図５のステップ１２では、燃焼室５の吸気弁閉時期ＩＶＣにおける容積ＶＩＶＣを、吸気弁閉時期をパラメータとする吸気弁閉時期容積のテーブルを検索することにより求めたが、ここではＭＢＴＣＹＣＬをパラメータとする前回燃焼開始時期容積のテーブルを検索することにより、燃焼室５の前回燃焼開始時期ＭＢＴＣＹＣＬにおける容積Ｖ０を求めればよい。

ステップ１６３では燃焼開始時期における有効圧縮比Ｅｃを算出する。有効圧縮比Ｅｃは無次元の値であり、次式に示すように燃焼室５の燃焼開始時期における容積Ｖ０を燃焼室５の吸気弁閉時期における容積ＶＩＶＣで除した値である。

Ｅｃ＝ｆ７（Ｖ０−ＶＤＥＰ、ＶＩＶＣ）
＝（Ｖ０−ＶＤＥＰ）／ＶＩＶＣ …（１２）
ただし、ＶＤＥＰ：燃焼室５内のデポジット量、
（１２）式では燃焼室５の燃焼開始時期における容積Ｖ０から、燃焼室５内におけるデポジットの堆積による燃焼室５の容積減少分を表すデポジット量ＶＤＥＰを減算している。これは次の理由による。燃焼室５内にデポジットが堆積して燃焼室５の容積が減少している状態になるとそのぶん圧縮比が上昇する。従って、燃焼室５内にデポジットが堆積して燃焼室５の容積が減少している状態になっても、燃焼室５内にデポジットが堆積していないときの燃焼室容積（ステップ１６２のＶ０）を用いて有効圧縮比Ｅｃを算出したのでは、燃焼室容積の減少に伴う誤差が生じるので、これを避けるためである。燃焼室５内のデポジット量ＶＤＥＰの算出については後述する（図３７参照）。

ステップ１６４では吸気弁閉時期ＩＶＣから燃焼開始時期に至る間の燃焼室５内の温度上昇率ＴＣＯＭＰを次式に示すように有効圧縮比Ｅｃに基づいて算出する。

ＴＣＯＭＰ＝ｆ８（Ｅｃ）＝Ｅｃ＾（κ−１） …（１３）
ただし、κ：比熱比、
（１３）式は断熱圧縮されるガスの温度上昇率の式である。なお、（１３）式右辺の「＾」は累乗計算を表している。この記号は後述する式でも使用する。

κは断熱圧縮されるガスの定圧比熱を定容比熱で除した値で、断熱圧縮されるガスが空気であればκ＝１．４であり、簡単にはこの値を用いればよい。ただし、混合気に対してκの値を実験的に求めることで、一層の算出精度の向上が可能である。

図１１は（１３）式を図示したものである。従って、このような特性のテーブルを予めエンジンコントローラ３１のメモリに格納しておき、有効圧縮比Ｅｃに基づき当該テーブルを検索することにより温度上昇率ＴＣＯＭＰを求めることも可能である。

ステップ１６５では、燃焼室５の燃焼開始時期における温度Ｔ０［Ｋ］を、燃焼室５の吸気弁閉時期における温度ＴＩＮＩに温度上昇率ＴＣＯＭＰを乗じることで、つまり
Ｔ０＝ＴＩＮＩ×ＴＣＯＭＰ …（１４）
の式により算出する。

ステップ１６６、１６７はステップ１６４、１６５と同様である。すなわち、ステップ１６６では吸気弁閉時期ＩＶＣから燃焼開始時期に至る間の燃焼室５内の圧力上昇率ＰＣＯＭＰを次式に示すように有効圧縮比Ｅｃに基づいて算出する。

ＰＣＯＭＰ＝ｆ９（Ｅｃ）＝Ｅｃ＾κ…（４１）
ただし、κ：比熱比、
（４１）式も（１３）式と同じに断熱圧縮されるガスの圧力上昇率の式である。（４１）式右辺の「＾」も（１３）式と同じに累乗計算を表している。

κは上記（１３）式で用いている値と同じで、断熱圧縮されるガスが空気であればκ＝１．４であり、簡単にはこの値を用いればよい。ただし、混合気に対してその組成、温度からκの値を求めることで、一層の算出精度の向上が可能である。

図１１と同様の特性のテーブルを予めエンジンコントローラ３１のメモリに格納しておき、有効圧縮比Ｅｃに基づき当該テーブルを検索することにより圧力上昇率ＰＣＯＭＰを求めることも可能である。

ステップ１６７では、燃焼室５の燃焼開始時期における圧力Ｐ０［Ｐａ］を、燃焼室５の吸気弁閉時期における圧力ＰＩＮＩに圧力上昇率ＰＣＯＭＰを乗じることで、つまり
Ｐ０＝ＰＩＮＩ×ＰＣＯＭＰ …（４２）
の式により算出する。

ステップ１６８では、初期燃焼期間における層流燃焼速度ＳＬ１［ｍ／ｓｅｃ］を次式（公知）により算出する。

ＳＬ１＝ｆ１０（Ｔ０、Ｐ０）
＝ＳＬstd×（Ｔ０／Ｔstd）^2.18×（Ｐ０／Ｐstd）^-0.16 …（１５）
ただし、Ｔstd ：基準温度［Ｋ］、
Ｐstd ：基準圧力［Ｐａ］、
ＳＬstd：基準温度Ｔstdと基準圧力Ｐstdにおける基準層流燃焼
速度［ｍ／ｓｅｃ］、
Ｔ０ ：燃焼室５の燃焼開始時期における温度［Ｋ］、
Ｐ０ ：燃焼室５の燃焼開始時期における圧力［Ｐａ］、
層流燃焼速度（層流火炎速度）は気体の流れがない状態での火炎の伝播速度のことであり、燃焼室５内の圧縮速度、燃焼室５内の吸気流速に因らず、燃焼室５の温度及び圧力の関数となることが知られていることから、初期燃焼期間における層流燃焼速度を燃焼開始時温度Ｔ０と燃焼開始時圧力Ｐ０の関数として、また後述するように主燃焼期における層流燃焼速度を圧縮上死点時温度ＴＴＤＣと圧縮上死点圧力ＰＴＤＣの関数としている。これは、層流燃焼速度は一般的に、エンジン負荷、燃焼室５内の不活性ガス率、吸気弁閉時期、比熱比、吸気温度により変化するのであるが、これらは燃焼室５内の温度Ｔと圧力Ｐに影響する因子であるので、層流燃焼速度は最終的に燃焼室５内の温度Ｔと圧力Ｐにより規定できるとするものである。

上記の（１５）式において基準温度Ｔstdと基準圧力Ｐstdと基準層流燃焼速度ＳＬstdは実験により予め定められる値である。

燃焼室５の通常の圧力である２ｂａｒ以上の圧力下では、（１５）式の圧力項（Ｐ０／Ｐstd）^-0.16は小さな値となる。従って、圧力項（Ｐ０／Ｐstd）^-0.16を一定値として、基準層流燃焼速度ＳＬstdを基準温度Ｔstdのみで規定することも可能である。

従って、基準温度Ｔstdが５５０［Ｋ］で、基準層流燃焼速度ＳＬstdが１．０［ｍ／ｓｅｃ］で、圧力項が０．７である場合の燃焼開始時期における温度Ｔ０と層流燃焼速度ＳＬ１との関係は近似的に次式で定義することができる。

ＳＬ１＝ｆ１１（Ｔ０）
＝１．０×０．７×（Ｔ０／５５０）^2.18 …（１６）
ステップ１６９では、初期燃焼期間におけるガス流動の乱れ強さＳＴ１を算出する。このガス流動の乱れ強さＳＴ１は無次元の値であり、燃焼室５に流入する新気の流速と燃料インジェクタ２１の噴射燃料のペネトレーションとに依存する。

燃焼室５に流入する新気の流速は、吸気通路の形状と、吸気弁１５の作動状態と、吸気弁１５を設ける吸気ポート４の形状に依存する。噴射燃料のペネトレーションは燃料インジェクタ２１の噴射圧力と、燃料噴射期間と、燃焼噴射タイミングに依存する。

最終的に、初期燃焼期間におけるガス流動の乱れ強さＳＴ１は、エンジン回転速度ＮＲＰＭの関数として次式で表すことができる。

ＳＴ１＝ｆ１２（ＮＲＰＭ）＝Ｃ１×ＮＲＰＭ …（１７）
ただし、Ｃ１：定数、
乱れ強さＳＴ１を回転速度ＮＲＰＭをパラメータとするテーブルから求めることも可能である。

ステップ１７０では層流燃焼速度Ｓ１と乱れ強さＳＴ１から、初期燃焼期間におけるガスの燃焼速度ＦＬＡＭＥ１［ｍ／ｓｅｃ］を次式により算出する。

ＦＬＡＭＥ１＝ＳＬ１×ＳＴ１ …（１８）
燃焼室５内にガス乱れがあるとガスの燃焼速度が変化する。（１８）式はこのガス乱れに伴う燃焼速度への寄与（影響）を考慮したものである。

ステップ１７１では、次式により初期燃焼期間ＢＵＲＮ１［ｄｅｇ］を算出する。

ＢＵＲＮ１＝｛（ＮＲＰＭ×６）×ＢＲ１×（Ｖ０−ＶＤＥＰ）｝
／（ＲＰＲＯＢＡ×ＡＦ１×ＦＬＡＭＥ１） …（１９）
ただし、ＶＤＥＰ：燃焼室５内のデポジット量、
ＡＦ１ ：火炎核の反応面積（固定値）［ｍ²］、
この（１９）式および後述する（２２）式は、燃焼ガス質量を燃焼速度で割ると燃焼期間が得られるとする次の基本式より導いたものであるが、（１９）、（２２）式右辺の分子、分母ががただちに燃焼ガス質量、燃焼速度を表すものではない。

燃焼期間［ｓｅｃ］＝シリンダ内総質量［ｇ］／（未燃ガス密度［ｇ／ｍ³］
×火炎表面積［ｍ²］×火炎速度［ｍ／ｓｅｃ］）
…（補１）
（補１）式右辺分母の未燃ガス密度は、未燃ガス質量［ｇ］を未燃ガス体積［ｍ³］で割った値であるので、従来装置のように質量に相当する充填効率ＩＴＡＣのみの関数では未燃ガス密度を正確に計算できているとはいえない。そこで、（補１）式に対して実験結果とを照らし合わせつつ所定の近似を導入して初めて得られたのが上記（１９）式及び後述する（２２）式に示す実験式である。

ここで、（１９）式右辺のＢＲ１は燃焼開始時期より初期燃焼期間ＢＵＲＮ１の終了時期までの燃焼質量割合の変化代であり、ここではＢＲ１＝２％に設定している。（１９）式右辺の（ＮＲＰＭ×６）は単位をｒｐｍからクランク角（ｄｅｇ）に変換するための処理である。火炎核の反応面積ＡＦ１は実験的に設定される。

また、初期燃焼期間中はほぼ燃焼室容積は変わらないとみなすことができる。従って、初期燃焼期間ＢＵＲＮ１を算出するに際して最初の燃焼室容積である燃焼開始時の燃焼室容積Ｖ０を採用している。

さらに（１９）式では燃焼室５の燃焼開始時期における容積Ｖ０から、燃焼室５内のデポジット量ＶＤＥＰを減算している。これは次の理由による。燃焼室５内にデポジットが堆積して燃焼室５の容積が減少している状態になるとそのぶん圧縮比が上昇して燃焼の速度が上昇し燃焼期間が短くなる。従って、燃焼室５内にデポジットが堆積して燃焼室５の容積が減少している状態になっても、燃焼室５内にデポジットが堆積していないときの燃焼室容積（ステップ１６２のＶ０）を用いて初期燃焼期間ＢＵＲＮ１を算出したのでは、圧縮比の上昇に伴う誤差が生じるので、これを避けるためである。燃焼室５内のデポジット量ＶＤＥＰの算出については後述する（図３７参照）。

次に図１２のフローに移ると、ステップ１８１では図１０のステップ１６１と同様に、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、図５のステップ１３で算出されている燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出されている燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、図５のステップ１５で算出されている反応確率ＲＰＲＯＢＡ［％］を読み込み、さらにシリンダ新気量ＭＡＣＹＬ［ｇ]、目標当量比ＴＦＢＹＡ、内部不活性ガス量ＭＲＥＳ［ｇ］、外部不活性ガス量ＭＥＧＲ［ｇ］を読み込む。

ここで、図１には外部ＥＧＲ装置は示していないが、図１２に関する限り外部ＥＧＲ装置を備えているエンジンを前提として説明する。この場合に、外部不活性ガス量ＭＥＧＲは例えば公知の手法（特開平１０−１４１１５０号公報参照）を用いて算出すればよい。なお、図１に示す本実施形態のように外部ＥＧＲ装置を備えていないエンジンを対象とするときには外部不活性ガス量ＭＥＧＲ＝０で扱えば足りる。シリンダ新気量ＭＡＣＹＬ、内部不活性ガス量ＭＲＥＳの算出については図１４以降で後述する。

ステップ１８２、１８３は図１０のステップ１６３、１６４と同様である。すなわち、ステップ１８２で圧縮上死点時期における有効圧縮比Ｅｃ ２を算出する。有効圧縮比Ｅｃ ２も上記（１２）式の有効圧縮比Ｅｃと同様に無次元の値であり、次式に示すように燃焼室５の圧縮上死点時における容積ＶＴＤＣを燃焼室５の吸気弁閉時期における容積ＶＩＶＣで除した値である。

Ｅｃ ２＝ｆ１３（ＶＴＤＣ、ＶＩＶＣ）＝ＶＴＤＣ／ＶＩＶＣ
…（４３）
（４３）式において燃焼室５の圧縮上死点時における容積ＶＴＤＣは運転条件によらず一定であり、予めエンジンコントローラ３１のメモリに格納しておけばよい。

ステップ１８３では吸気弁閉時期ＩＶＣから圧縮上死点に至る間の燃焼室５内の断熱圧縮による温度上昇率ＴＣＯＭＰ ２を次式に示すように有効圧縮比Ｅｃ ２に基づいて算出する。

ＴＣＯＭＰ ２＝ｆ１４（Ｅｃ ２）
＝Ｅｃ ２＾（κ−１）…（４４）
ただし、κ：比熱比、
図１１と同様の特性のテーブルを予めエンジンコントローラ３１のメモリに格納しておき、有効圧縮比Ｅｃ ２から当該テーブルを検索することにより温度上昇率ＴＣＯＭＰ ２を求めることも可能である。

ステップ１８４ではシリンダ新気量ＭＡＣＹＬ、目標当量比ＴＦＢＹＡ、内部不活性ガス量ＭＲＥＳ、外部不活性ガス量ＭＥＧＲから次式により燃焼室５の総ガス質量ＭＧＡＳ［ｇ］を算出する。

ＭＧＡＳ＝ＭＡＣＹＬ×(１＋ＴＦＢＹＡ／１４．７)＋ＭＲＥＳ＋ＭＥＧＲ
…（４５）
（４５）式右辺の括弧内の「１」は新気分、「ＴＦＢＹＡ／１４．７」は燃料分である。

ステップ１８５ではこの燃焼室５の総ガス質量ＭＧＡＳと、シリンダ新気量ＭＡＣＹＬ、目標当量比ＴＦＢＹＡを用い、次式により混合気の燃焼による温度上昇量（燃焼上昇温度）ＴＢＵＲＮ［Ｋ］を算出する。

ＴＢＵＲＮ＝｛ＭＡＣＹＬ×ＴＦＢＹＡ／１４．７×ＢＲｋ×Ｑ｝
／（Ｃｖ×ＭＧＡＳ）…（４６）
ただし、Ｑ ：燃料の定発熱量、
ＢＲｋ：シリンダ内燃料の燃焼質量割合、
Ｃｖ ：定積比熱、
（４６）式右辺の分子はシリンダ内燃料による発生総熱量［Ｊ］、分母は単位発生熱量当たりの温度上昇率［Ｊ／Ｋ］を意味している。すなわち、（４６）式は熱力学の公式に当てはめた近似式である。

ここで、シリンダ内燃料の燃焼質量割合ＢＲｋとしては予め実験等で適合しておく。簡易的には例えば６０％／２＝３０％を設定する。これは、本実施形態では燃焼質量割合ＢＲが約６０％に達するまでを燃焼期間として扱うので、そのちょうど中間の３０％をＢＲｋとして設定するものである。

燃料の定発熱量Ｑは燃料の種類により異なる値であるので、燃料の種類に応じ予め実験等で求めておく。定積比熱Ｃｖは２〜３の値であり予め実験等で代表値を適合しておく。ただし、混合気に対してその組成、温度から定積比熱Ｃｖの値を求めることで、一層の算出精度の向上が可能である。

ステップ１８６では、燃焼室５の圧縮上死点における温度ＴＴＤＣ［Ｋ］を、燃焼室５の吸気弁閉時期における温度ＴＩＮＩに圧縮上死点までの温度上昇率ＴＣＯＭＰ ２を乗じその乗算値に上記の燃焼上昇温度ＴＢＵＲＮを加算することで、つまり次式により算出する。

ＴＴＤＣ＝ＴＩＮＩ×ＴＣＯＭＰ ２＋ＴＢＵＲＮ
…（４７）
ステップ１８７では、この燃焼室５の圧縮上死点における温度ＴＴＤＣと容積ＶＴＤＣ及び燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ、容積ＶＩＶＣ及び温度ＴＩＮＩから次式により燃焼室５の圧縮上死点における圧力ＰＴＤＣ［Ｋ］を算出する。

ＰＴＤＣ＝ＰＩＮＩ×ＶＩＶＣ×ＴＴＤＣ／（ＶＴＤＣ×ＴＩＮＩ）
…（４８）
（４８）式は状態方程式を用いて得たものである。すなわち、吸気弁閉時期における圧力、容積及び温度（ＰＩＮＩ、ＶＩＶＣ、ＴＩＮＩ）を用いて次の状態方程式が成立する。

ＰＩＮＩ×ＶＩＶＣ＝ｎ・Ｒ・ＴＩＮＩ…（補２）
ただし、ｎ：モル数、
Ｒ：ガス定数、
圧縮上死点近傍では容積はほぼ等しいので、圧縮上死点での圧力、容積及び温度（ＰＴＤＣ、ＶＴＤＣ、ＴＴＤＣ）を用いて次の状態方程式が成立する。

ＰＴＤＣ×ＶＴＤＣ＝ｎ・Ｒ・ＴＴＤＣ…（補３）
この（補３）式と上記（補２）との両式からｎ・Ｒを消去しＰＴＤＣについて解くと、上記（４８）式が得られる。

ステップ１８８では図１０のステップ１６８と同様にして、次式（公知）により、主燃焼期間における層流燃焼速度ＳＬ２［ｍ／ｓｅｃ］を算出する。

ＳＬ２＝ｆ１５（ＴＴＤＣ、ＰＴＤＣ）
＝ＳＬstd×（ＴＴＤＣ／Ｔstd）^2.18×（ＰＴＤＣ／Ｐstd）^-0.16
…（４９）
ただし、Ｔstd ：基準温度［Ｋ］、
Ｐstd ：基準圧力［Ｐａ］、
ＳＬstd：基準温度Ｔstdと基準圧力Ｐstdにおける基準層流燃焼速度
［ｍ／ｓｅｃ］、
ＴＴＤＣ：燃焼室５の圧縮上死点における温度［Ｋ］、
ＰＴＤＣ：燃焼室５の圧縮上死点における圧力［Ｐａ］、
（４９）式の解説は上記（１６）式と同様ある。すなわち、（４９）式の基準温度Ｔstdと基準圧力Ｐstdと基準層流燃焼速度ＳＬstdは実験により予め定められる値である。燃焼室５の通常の圧力である２ｂａｒ以上の圧力下では、（４９）式の圧力項（ＰＴＤＣ／Ｐstd）^-0.16は小さな値となる。従って、圧力項（ＰＴＤＣ／Ｐstd）^-0.16を一定値として、基準層流燃焼速度ＳＬstdを基準温度Ｔstdのみで規定することも可能である。よって、基準温度Ｔstdが５５０［Ｋ］で、基準層流燃焼速度ＳＬstdが１．０［ｍ／ｓｅｃ］で、圧力項が０．７である場合の圧縮上死点における温度ＴＴＤＣと層流燃焼速度ＳＬ２との関係は近似的に次式で定義することができる。

ＳＬ２＝ｆ１６（ＴＴＤＣ）
＝１．０×０．７×（ＴＴＤＣ／５５０）^2.18
…（５０）
ステップ１８９期間では主燃焼期間におけるガス流動の乱れ強さＳＴ２を算出する。このガス流動の乱れ強さＳＴ２も初期燃焼期間におけるガス流動の乱れ強さＳＴ１と同様に、エンジン回転速度ＮＲＰＭの関数として次式で表すことができる。

ＳＴ２＝ｆ１７（ＮＲＰＭ）＝Ｃ２×ＮＲＰＭ …（２０）
ただし、Ｃ２：定数、
乱れ強さＳＴ２を回転速度をパラメータとするテーブルから求めることも可能である。

ステップ１９０では、層流燃焼速度ＳＬ２［ｍ／ｓｅｃ］と主燃焼期間におけるガス流動の乱れ強さＳＴ２とから、主燃焼期間における燃焼速度ＦＬＡＭＥ２［ｍ／ｓｅｃ］を次式により算出する。

ＦＬＡＭＥ２＝ＳＬ２×ＳＴ２ …（２１）
ただし、ＳＬ２：層流燃焼速度［ｍ／ｓｅｃ］、
（２１）式は（１８）式と同様、ガス乱れに伴う燃焼速度への寄与を考慮したものである。

ステップ１９１では、主燃焼期間ＢＵＲＮ２［ｄｅｇ］を（１９）式に類似した次式で算出する。

ＢＵＲＮ２＝｛（ＮＲＰＭ×６）×（ＢＲ２×ＶＴＤＣ）｝
／（ＲＰＲＯＢＡ×ＡＦ２×ＦＬＡＭＥ２） …（２２）
ただし、ＡＦ２：火炎核の反応面積［ｍ²］、
ここで、（２２）式右辺のＢＲ２は主燃焼期間の開始時期より終了時期までの燃焼質量割合の変化代である。初期燃焼期間の終了時期に燃焼質量割合ＢＲが２％になり、その後、主燃焼期間が開始し、燃焼質量割合ＢＲが６０％に達して主燃焼期間が終了すると考えているので、ＢＲ２＝６０％−２％＝５８％を設定している。ＡＦ２は火炎核の成長行程における平均の反応面積であり、（１９）式のＡＦ１と同様に、予め実験的に定めた固定値とする。

主燃焼期間では圧縮上死点を挟んで燃焼室容積が変化する。つまり、主燃焼期間の開始時期と、主燃焼期間の終了時期のほぼ中央に圧縮上死点位置が存在するとみなすことができる。また、圧縮上死点付近ではクランク角が変化しても燃焼室容積があまり変化しない。そこで主燃焼期間での燃焼室容積としてはこの圧縮上死点での燃焼室容積ＶＴＤＣで代表させることとしている。

図１３は基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１０、図１２のうち遅く実行されるフローに続けて実行する。

ステップ４１では、図１０のステップ１７１で算出されている初期燃焼期間ＢＵＲＮ１、図１２のステップ１９１で算出されている主燃焼期間ＢＵＲＮ２、図５のステップ１７で算出されている点火時期無駄時間相当クランク角ＩＧＮＤＥＡＤ、図５のステップ１６で算出されている基準クランク角θＰＭＡＸを読み込む。

ステップ４２では、初期燃焼期間ＢＵＲＮ１と主燃焼期間ＢＵＲＮ２の合計を燃焼期間ＢＵＲＮ［ｄｅｇ］として算出する。

ステップ４３では次式により基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出する。

ＭＢＴＣＡＬ＝ＢＵＲＮ−θＰＭＡＸ＋ＩＧＮＤＥＡＤ …（２３）
ステップ４４では、この基本点火時期ＭＢＴＣＡＬから点火無駄時間相当クランク角ＩＧＮＤＥＡＤを差し引いた値を前回燃焼開始時期ＭＢＴＣＹＣＬ［ｄｅｇＢＴＤＣ］として算出する。

このようにして算出した基本点火時期ＭＢＴＣＡＬは、点火時期指令値として点火レジスタに移され、実際のクランク角がこの点火時期指令値と一致したタイミングでエンジンコントローラ３１より一次電流を遮断する点火信号が点火コイル１３に出力される。

また、今サイクルの点火時期指令値としてステップ４３で算出された基本点火時期ＭＢＴＣＡＬが用いられたとすると、次サイクルの点火時期になるまでの間、ステップ４４で算出された前回燃焼開始時期ＭＢＴＣＹＣＬが図１０のステップ１６２において用いられる。

次に、図１４は燃焼室５内の内部不活性ガス率ＭＲＥＳＦＲを算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。このフローは上記図５のフローに先立って実行する。

ステップ５１ではエアフローメータ３２の出力と目標当量比ＴＦＢＹＡを読み込む。ステップ５２ではエアフロメータ３２の出力に基づいて、燃焼室５に流入する新気量（シリンダ新気量）ＭＡＣＹＬを算出する。このシリンダ新気量ＭＡＣＹＬの算出方法については公知の方法を用いればよい（特開２００１−５００９１号公報参照）。

ステップ５３では、燃焼室５内の内部不活性ガス量ＭＲＥＳを算出する。この内部不活性ガス量ＭＲＥＳの算出については、図１５のフローにより説明する。

図１５（図１４ステップ５３のサブルーチン）においてステップ６１では、燃焼室５内の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬを算出する。この不活性ガス量ＭＲＥＳＣＹＬの算出についてはさらに図１６のフローにより説明する。

図１６（図１５ステップ６１のサブルーチン）においてステップ７１では、排気弁閉時期ＥＶＣ［ｄｅｇＢＴＤＣ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、圧力センサ４６により検出される排気圧力ＰＥＸＨ［ｋＰａ］を読み込む。

ここで、吸気弁閉時期ＩＶＣが吸気ＶＴＣ機構２７に与える指令値から既知であったように、排気弁閉時期ＥＶＣも排気ＶＴＣ機構２８に与える指令値から既知である。

ステップ７２では燃焼室５の排気弁閉時期ＥＶＣにおける容積ＶＥＶＣを算出する。これは吸気弁閉時期ＩＶＣにおける容積ＶＩＶＣと同様に、排気弁閉時期をパラメータとするテーブルを検索することにより求めればよい。すなわち、排気弁ＶＴＣ機構２８を備える場合には、排気弁閉時期ＥＶＣから図２３に示すテーブルを検索することにより、燃焼室５の排気弁閉時期ＥＶＣにおける容積ＶＥＶＣを求めればよい。排気ＶＴＣ機構２８を備えないときには定数で与えることができる。

また、図示しないが圧縮比を変化させる機構を有する場合には、圧縮比の変化量に応じた排気弁閉時期における燃焼室容積ＶＥＶＣをテーブルから求める。排気ＶＴＣ機構２８に加えて圧縮比を変化させる機構をも有する場合には、排気弁閉時期と圧縮比変化量とに応じたマップを検索することにより排気弁閉時期における燃焼室容積を求める。

ステップ７３では、目標当量比ＴＦＢＹＡから図２４に示すテーブルを検索することにより、燃焼室５内の不活性ガスのガス定数ＲＥＸを求める。図２４に示すように、不活性ガスのガス定数ＲＥＸは目標当量比ＴＦＢＹＡが１．０のとき、つまり理論空燃比のとき最も小さく、これより大きくても小さくても大きくなる。

ステップ７４では、排気温度ＴＥＸＨに基づいて燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣを推定する。簡単には排気温度ＴＥＸＨをそのままＴＥＶＣとおけばよい。なお、燃焼室５の排気弁閉時期における温度ＴＥＶＣは、インジェクタ２１の燃料噴射量に応じた熱量により変化するため、このような特性をも加味すれば、ＴＥＶＣの算出精度が向上する。

ステップ７５では、排気圧力ＰＥＸＨに基づいて燃焼室５の排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣを算出する。簡単には排気圧力ＰＥＸＨをＰＥＶＣと置けばよい。

ステップ７６では、燃焼室５の排気弁閉時期ＥＶＣにおける容積ＶＥＶＣ、排気弁閉時期ＥＶＣにおける温度ＴＥＶＣ、排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣ及び不活性ガスのガス定数ＲＥＸから、燃焼室５の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬを次式により算出する。

ＭＲＥＳＣＹＬ＝（ＰＥＶＣ×ＶＥＶＣ）／（ＲＥＸ×ＴＥＶＣ） …（２４）
このようにして燃焼室５の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬの算出を終了したら図１５に戻り、ステップ６２で吸排気弁１５、１６のオーバーラップ（図では「Ｏ／Ｌ」と略記する）中に排気側から吸気側へ吹き返す不活性ガス量であるオーバーラップ中吹き返し不活性ガス量ＭＲＥＳＯＬを算出する。

この不活性ガス量ＭＲＥＳＯＬの算出については図１７のフローにより説明する。

図１７（図１５ステップ６２のサブルーチン）においてステップ８１では、吸気弁開時期ＩＶＯ［ｄｅｇＢＴＤＣ］と、排気弁閉時期ＥＶＣ［ｄｅｇＢＴＤＣ］、図１６のステップ７４で算出されている燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣを読み込む。

ここで、吸気弁開時期ＩＶＯは、吸気弁閉時期ＩＶＣより吸気弁１５の開き角だけ前の時期となるので、吸気弁閉時期ＩＶＣより吸気弁１５の開き角（予め分かっている）とから求めることができる。

ステップ８２では吸気弁開時期ＩＶＯと排気弁閉時期ＥＶＣとから、吸排気弁のオーバーラップ量ＶＴＣＯＬ［ｄｅｇ］を次式により算出する。

ＶＴＣＯＬ＝ＩＶＯ＋ＥＶＣ …（２５）
例えば、吸気ＶＴＣ機構２７用アクチュエータへの非通電時に吸気弁開時期ＩＶＯが吸気上死点位置にあり、吸気ＶＴＣ機構２７用アクチュエータへの通電時に吸気弁開時期が吸気上死点より進角する特性であり、かつ排気ＶＴＣ機構２８用アクチュエータへの非通電時に排気弁閉時期ＥＶＣが排気上死点にあり、排気弁ＶＴＣ機構２８用アクチュエータへの通電時に排気弁閉時期ＥＶＣが排気上死点より進角する特性である場合には、ＩＶＯとＥＶＣの合計が吸排気弁のオーバーラップ量ＶＴＣＯＬとなる。

ステップ８３では、吸排気弁のオーバーラップ量ＶＴＣＯＬから、図２５に示すテーブルを検索することによりオーバーラップ中の積算有効面積ＡＳＵＭＯＬを算出する。図２５に示すようにオーバーラップ中の積算有効面積ＡＳＵＭＯＬは吸排気弁のオーバーラップ量ＶＴＣＯＬが大きくなるほど大きくなる値である。

ここで、図２６は、吸排気弁のオーバーラップ中の積算有効面積ＡＳＵＭＯＬの説明図であり、横軸はクランク角、縦軸は吸気弁１２と排気弁１５とのそれぞれの開口面積を示している。オーバーラップ中の任意の時点における有効開口面積は、排気弁開口面積と吸気弁開口面積とのうち小さい方とする。オーバーラップ中の全期間における積算有効面積ＡＳＵＭＯＬは、吸気弁１５及び排気弁１６が開いている期間の積分値（図中の斜線部）である。

このようにオーバーラップ中積算有効面積ＡＳＵＭＯＬを算出することで、吸気弁１５と排気弁１６とのオーバーラップ量を１つのオリフィス（流出孔）であると近似することができ、排気系の状態と吸気系の状態とからこの仮想オリフィスを通過するガス流量を簡略的に算出し得る。

ステップ８４では、目標当量比ＴＦＢＹＡと、燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣとから、図２７に示すマップを検索することにより、燃焼室５に残留する不活性ガスの比熱比ＳＨＥＡＴＲを算出する。図２７に示したように、燃焼室に残留する不活性ガスの比熱比ＳＨＥＡＴＲは目標当量比ＴＦＢＹＡが１．０の近傍にあるときが最も小さくなり、それより大きくても小さくても大きくなる。また、目標当量比ＴＦＢＹＡが一定の条件では、燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣが高くなるほど小さくなる。

ステップ８５では過給判定フラグＴＢＣＲＧ及びチョーク判定フラグＣＨＯＫＥを設定する。この過給判定フラグＴＢＣＲＧ及びチョーク判定フラグＣＨＯＫＥの設定については図１８のフローにより説明する。

図１８（図１７ステップ８５のサブルーチン）においてステップ１０１では、吸気圧力センサ４４により検出される吸気圧力ＰＩＮと、図１６のステップ７５で算出されている燃焼室５の排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣを読み込む。

ステップ１０２では、吸気圧力ＰＩＮと、燃焼室５の排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣとから、次式により吸気排気圧力比ＰＩＮＢＹＥＸを算出する。

ＰＩＮＢＹＥＸ＝ＰＩＮ／ＰＥＶＣ …（２６）
この吸気排気圧力比ＰＩＮＢＹＥＸは無名数であり、これと１をステップ１０３で比較する。吸気排気圧力比ＰＩＮＢＹＥＸが１以下の場合には過給無しと判断し、ステップ１０４に進んで過給判定フラグＴＢＣＲＧ（ゼロに初期設定）＝０とする。

吸気排気圧力比ＰＩＮＢＹＥＸが１より大きい場合には過給有りと判断し、ステップ１０５へ進んで過給判定フラグＴＢＣＲＧ＝１とする。

ステップ１０６では、図１４のステップ５１で読み込まれている目標当量比ＴＦＢＹＡから図２８に示すテーブルを検索することにより、混合気の比熱比ＭＩＸＡＩＲＳＨＲを求め、これをステップ１０７で不活性ガスの比熱比ＳＨＥＡＴＲと入れ換える。図２８に示したように、混合気の比熱比ＭＩＸＡＩＲＳＨＲは、目標当量比ＴＦＢＹＡが小さくなるほど大きくなる値である。

ステップ１０６、１０７において、不活性ガスの比熱比ＳＨＥＡＴＲを混合気の比熱比ＭＩＸＡＩＲＳＨＲに置き換えるのは、ターボ過給や慣性過給等の過給時を考慮したものである。すなわち、過給時には吸排気弁のオーバーラップ中のガス流れが吸気系から排気系へ向かう（吹き抜ける）ので、この場合においては、上記の仮想オリフィスを通過するガスの比熱比を不活性ガスの比熱比から混合気の比熱比に変更することで、吹き抜けるガス量を精度良く推定し、内部不活性ガス量を精度良く算出するためである。

ステップ１０８では、図１７のステップ８４または図１８のステップ１０６、１０７で算出している不活性ガスの比熱比ＳＨＥＡＴＲに基づき、最小と最大とのチョーク判定しきい値ＳＬＣＨＯＫＥＬ、ＳＬＣＨＯＫＥＨを次式により算出する。

ＳＬＣＨＯＫＥＬ＝｛２／（ＳＨＥＡＴＲ＋１）｝
＾｛ＳＨＥＡＴＲ／（ＳＨＥＡＴＲ−１）｝ …（２７ａ）
ＳＬＣＨＯＫＥＨ＝｛−２／（ＳＨＥＡＴＲ+１）｝
＾｛−ＳＨＥＡＴＲ／（ＳＨＥＡＴＲ−１）｝…（２７ｂ）
これらのチョーク判定しきい値ＳＬＣＨＯＫＥＬ、ＳＬＣＨＯＫＥＨは、チョークする限界値を算出している。

ステップ１０８において、（２７ａ）右辺、（２７ｂ）右辺の各累乗計算が困難な場合には、（２７ａ）、（２７ｂ）式の算出結果を、最小チョーク判定しきい値ＳＬＣＨＯＫＥＬのテーブルと最大チョーク判定しきい値ＳＬＣＨＯＫＥＨのテーブルとしてそれぞれエンジンコントローラ３１のメモリに予め記憶しておき、不活性ガスの比熱比ＳＨＥＡＴＲから当該テーブルを検索することにより求めてもよい。

テップ１０９では、吸気排気圧力比ＰＩＮＢＹＥＸが、最小チョーク判定しきい値ＳＬＣＨＯＫＥＬ以上でかつ最大チョーク判定しきい値ＳＬＣＨＯＫＥＨ以下の範囲内にあるか否か、すなわちチョーク状態にないか否かを判定する。吸気排気圧力比ＰＩＮＢＹＥＸが範囲内にある場合にはチョーク無しと判断し、ステップ１１０に進んでチョーク判定フラグＣＨＯＫＥ（ゼロに初期設定）＝０とする。

吸気排気圧力比Ｐ１ＮＢＹＥＸが範囲内にない場合にはチョーク有りと判断し、ステップ１１１に進んでチョーク判定フラグＣＨＯＫＥ＝１とする。

このようにして過給判定フラグとチョーク判定フラグの設定を終了したら図１７に戻り、ステップ８６〜８８で次の４つの場合分けを行う。

〈１〉過給判定フラグＴＢＣＲＧ＝０かつチョーク判定フラグＣＨＯＫＥ＝０のとき
〈２〉過給判定フラグＴＢＣＲＧ＝０かつチョーク判定フラグＣＨＯＫＥ＝０のとき
〈３〉過給判定フラグＴＢＣＲＧ＝０かつチョーク判定フラグＣＨＯＫＥ＝１のとき
〈４〉過給判定フラグＴＢＣＲＧ＝１かつチョーク判定フラグＣＨＯＫＥ＝０のとき
そして、上記〈１〉のときにはステップ８９に進んで、過給無しかつチョーク無し時のオーバーラップ中の平均吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ１を、上記〈２〉のときにはステップ９０に進んで過給無しかつチョーク有り時のオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ２を、上記〈３〉のときにはステップ９１に進んで過給有りかつチョーク無し時のオーバーラップ中の平均吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ３を、上記〈４〉のときにはステップ９２に進んで過給有りかつチョーク有り時の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ４をそれぞれ算出し、算出結果をオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。

ここで、過給無しかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ１の算出について図１９のフローにより説明する
図１９（図１７ステップ８９のサブルーチン）においてステップ１２１では、図１６のステップ７３、７５で算出されている不活性ガスのガス定数ＲＥＸ、燃焼室５の排気弁閉時期における圧力ＰＥＶＣを読み込む。

ステップ１２２では、不活性ガスのガス定数ＲＥＸと、図１７のステップ８１で読み込まれている燃焼室５の排気弁閉時期における温度ＴＥＶＣとに基づき、後述するガス流量の算出式に用いる密度項ＭＲＳＯＬＤを次式により算出する。

ＭＲＳＯＬＤ＝ＳＱＲＴ｛１／（ＲＥＸ×ＴＥＶＣ）｝ …（２８）
ここで、（２８）式右辺の「ＳＱＲＴ」はすぐ右のカッコ内の値の平方根を計算させる関数である。

なお、密度項ＭＲＳＯＬＤの平方根計算が困難な場合は、（２８）式の算出結果をマップとしてエンジンコントローラ３１のメモリに予め記憶しておき、ガス定数ＲＥＸと燃焼室５の排気弁閉時期における温度ＴＥＶＣとからそのマップを検索することにより求めてもよい。

ステップ１２３では、図１７のステップ８４で算出されている不活性ガスの比熱比ＳＨＥＡＴＲと、図１８のステップ１０２で算出されている吸気排気圧力比ＰＩＮＢＹＥＸとに基づき、後述するガス流量の算出式に用いる圧力差項ＭＲＳＯＬＰを次式により算出する。

ＭＲＳＯＬＰ＝ＳＱＲＴ［ＳＨＥＡＴＲ／（ＳＨＥＡＴＲ−１）
×｛ＰＴＮＢＹＥＸ＾（２／ＳＨＥＡＴＲ）
−ＰＴＮＢＹＥＸ＾（（ＳＨＥＡＴＲ＋１）／ＳＨＥＡＴＲ）｝］…（２９）
ステップ１２４では、これら密度項ＭＲＳＯＬＤ、圧力差項ＭＲＳＯＬＰと、燃焼室５の排気弁閉時期における圧力ＰＥＶＣとから、過給無しかつチョーク無し時のオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ１を次式（ガス流量の算出式）により算出し、その算出値をステップ１２５でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。

ＭＲＥＳＯＬｔｍｐ１＝１．４×ＰＥＶＣ×ＭＲＳＯＬＤ×ＭＲＳＯＬＰ…（３０）
次に、過給無しかつチョーク有り時の吹き返し不活性ガス流量の算出について図２０のフローにより説明する
図２０（図１７ステップ９０のサブルーチン）においてステップ１３１、１３２では、図１９のステップ１２１、１２２と同様にして、不活性ガスのガス定数ＲＥＸ、燃焼室５の排気弁閉時期における圧力ＰＥＶＣを読み込み、これらから前述の（２８）式により密度項ＭＲＳＯＬＤを算出する。

ステップ１３３では、図１７のステップ８４で算出されている不活性ガスの比熱比ＳＨＥＡＴＲに基づき、チョーク時圧力差項ＭＲＳＯＬＰＣを次式により算出する。

ＭＲＳＯＬＰＣ＝ＳＱＲＴ［ＳＨＥＡＴＲ×｛２／（ＳＨＥＡＴＲ＋１）｝ ＾｛（ＳＨＥＡＴＲ＋１）／〔ＳＨＥＡＴＲ−１）｝］…（３１）
なお、（３１）式の累乗計算と平方根計算とが困難な場合には、（３１）式の算出結果を、チョーク時圧力差項ＭＲＳＯＬＰＣのテーブルとしてエンジンコントローラ３１のメモリに予めに記憶しておき、不活性ガスの比熱比ＳＨＥＡＴＲからそのテーブルを検索することにより求めてもよい。

ステップ１３４では、これら密度項ＭＲＳＯＬＤ、チョーク時圧力差項ＭＲＳＯＬＰＣと、燃焼室５の排気弁閉時期における圧力ＰＥＶＣとから、過給無しかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ２を次式により算出し、その算出値をステップ１３５でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。

ＭＲＥＳＯＬｔｍｐ２＝ＰＥＶＣ×ＭＲＳＯＬＤ×ＭＲＳＯＬＰＣ …（３２）
次に、過給有りかつチョーク無し時の吹き返しガス流量の算出について図２１のフローにより説明する
図２１（図１７ステップ９１のサブルーチン）においてステップ１４１では、吸気圧力センサ４４により検出される吸気圧力ＰＩＮを読み込む。

ステップ１４２では、図１８のステップ１０６、１０７で算出されている不活性ガスの比熱比ＳＨＥＡＴＲと、図１８のステップ１０２で算出されている吸気排気圧力比ＰＩＮＢＹＥＸとから、過給時圧力差項ＭＲＳＯＬＰＴを次式により算出する。

ＭＲＳＯＬＰＴ＝ＳＱＲＴ［ＳＨＥＡＴＲ／（ＳＨＥＡＴＲ−１）
×｛ＰＩＮＢＹＥＸ＾（−２／ＳＨＥＡＴＲ）
−ＰＩＮＢＹＥＸ＾（−（ＳＨＥＡＴＲ＋１）／ＳＨＥＡＴＲ）｝］…（３３）
なお、（３３）式の累乗計算と平方根計算とが困難な場合は、（３３）式の算出結果を、過給時圧力差項ＭＲＳＯＬＰＴのマップとしてエンジンコントローラ３１のメモリに予め記憶しておき、不活性ガスの比熱比ＳＨＥＡＴＲと吸気排気圧力比ＰＩＮＢＹＥＸとからそのマップを検索することにより求めてもよい。

ステップ１４３では、この過給時圧力差項ＭＲＳＯＬＰＴと吸気圧力ＰＩＮとに基づいて、過給有りかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ３を次式により算出し、その算出値をステップ１４４でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。

ＭＲＥＳＯＬｔｍｐ３＝−０．１５２×ＰＩＮ×ＭＲＳＯＬＰＴ …（３４）
ここで、（３４）式の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ３は負の値とすることで、オーバーラップ中に吸気系から排気系へ吹き抜ける混合気のガス流量を表すことができる。

次に、過給有りかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出について図２２のフローにより説明する
図２２（図１７ステップ９２のサブルーチン）においてステップ１５１、１５２では、図２１のステップ１４１と同じく吸気圧力センサ４４により検出される吸気圧力ＰＩＮを読み込むと共に、図２０のステップ１３２と同じくチョーク時圧力差項ＭＲＳＯＬＰＣを前述の（３１）式により算出する。

ステップ１５３では、このチョーク時圧力差項ＭＲＳＯＬＰＣと吸気圧力ＰＩＮとに基づいて、過給有りかつチョーク有り時のオーバーラップ中吹き返しガス流量ＭＲＥＳＯＬｔｍｐ４を次式により算出し、その算出値をステップ１５４でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。

ＭＲＥＳＯＬｔｍｐ４＝−０．１０８×ＰＩＮ×ＭＲＳＯＬＰＣ …（３５）
ここで、（３５）式の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ４も、ＭＲＥＳＯＬｔｍｐ３と同様、負の値とすることで、オーバーラップ中に吸気側から排気側へ吹き抜ける混合気のガス流量を表すことができる。

このようにして、過給の有無とチョークの有無との組み合わせにより場合分けした、オーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐの算出を終了したら図１７に戻り、ステップ９３においてこのオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐとオーバーラップ期間中の積算有効面積ＡＳＵＭＯＬとから、次式によりオーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬを算出する。

ＭＲＥＳＯＬ＝（ＭＲＥＳＯＬｔｍＰ×ＡＳＵＭＯＬ×６０）
／（ＮＲＰＭ×３６０） …（３６）
このようにしてオーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬの算出を終了したら図１５に戻り、ステップ６３において燃焼室５内の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬと、このオーバーラップ中吹き返しガス量ＭＲＥＳＯＬとを加算して、つまり次式により内部不活性ガス量ＭＲＥＳを算出する。

ＭＲＥＳ＝ＭＲＥＳＣＹＬ＋ＭＲＥＳＯＬ …（３７）
前述のように、過給有り時にはオーバーラップ中吹き返し不活性ガス流量（ＭＲＥＳＯＬｔｍｐ３、ＭＲＥＳＯＬｔｍｐ４）が負となるため、上記（３６）式のオーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬも負となり、このとき（３７）式によれば、オーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬの分だけ内部不活性ガス量が減じられる。

このようにして内部不活性ガス量ＭＲＥＳの算出を終了したら図１４に戻り、ステップ５４においてこの内部不活性ガス量ＭＲＥＳと、目標当量比ＴＦＢＹＡとを用いて、次式により内部不活性ガス率ＭＲＥＳＦＲ（燃焼室５内の総ガス量に対する内部不活性ガス量の割合）を算出する。

ＭＲＥＳＦＲ＝ＭＲＥＳ
／｛ＭＲＥＳ＋ＭＡＣＹＬ×（１＋ＴＦＢＹＡ／１４．７）｝…（３８）
これで内部不活性ガス率ＭＲＥＳＦＲの算出を総て終了する。

このように本実施形態によれば、内部不活性ガス量ＭＲＥＳを、燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬと、吸排気弁のオーバーラップ中の吹き返しガス量ＭＲＥＳＯＬとで構成し（図１５のステップ６３参照）、この場合に、燃焼室５の排気弁閉時期における温度ＴＥＶ及び圧力ＰＥＶＣを算出し（図１６のステップ７４、７５）、これら温度ＴＥＶＣ、圧力ＰＥＶＣと不活性ガスのガス定数ＲＥＸとに基づいて状態方程式（上記（２４）式）により燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬを算出する（図１６のステップ７６参照）ようにしたので、特に、燃焼室５内部の状態量（ＰＥＶＣ、ＶＥＶＣ、ＴＥＶＣ）が刻々と変化する過渡運転時においても、運転条件に関わらず精度良く燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬを算出（推定）できる。

また、燃焼室５の排気弁閉時期における温度ＴＥＶＣ及び圧力ＰＥＶＣ、不活性ガスのガス定数ＲＥＸ及び比熱比ＳＨＥＡＴＲ、吸気圧力ＰＩＮに基づいてオーバーラップ中の吹き返し不活性ガス流量（ＭＲＥＳＯＬｔｍｐ１、ＭＲＥＳＯＬｔｍｐ２）を算出し（図１９、図２０参照）、このガス流量にオーバーラップ中の積算有効面積ＡＳＵＭＯＬを乗算して、オーバーラップ中の吹き返しガス量ＭＲＥＳＯＬを算出する（図１７のステップ９３参照）ようにしたので、精度良くオーバーラップ中吹き返しガス量ＭＲＥＳＯＬを算出（推定）できる。

このように、燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬ、オーバーラップ中吹き返しガス量ＭＲＥＳＯＬとも精度良く算出（推定）できると、これらの和である内部不活性ガス量ＭＲＥＳも精度良く算出（推定）できることになり、この精度良く推定することが可能となった内部不活性ガス量ＭＲＥＳに基づいて算出される内部不活性ガス率ＭＲＥＳＦＲを、点火時期の算出に用いる燃焼室５内の吸気弁閉時期ＩＶＣにおける温度ＴＩＮＩに活かすことで（図５のステップ１３参照）、燃焼室５内の吸気弁閉時期ＩＶＣにおける温度ＴＩＮＩを精度良く算出できる。また、精度良く推定することが可能となった内部不活性ガス量ＭＲＥＳを、燃料噴射量、バルブ開閉タイミング（オーバーラップ量）などにも活かすことで、エンジンを適切に制御することが可能である。

また、不活性ガスのガス定数ＲＥＸや不活性ガスの比熱比ＳＨＥＡＴＲは、目標当量比ＴＦＢＹＡに応じた値としているので（図２４、図２７参照）、理論空燃比を外れた空燃比での運転時（例えば理論空燃比よりもリーンな空燃比で運転を行うリーン運転時、冷間始動時のようにエンジンが元々不安定な状態を安定させるために理論空燃比の空燃比よりもリッチ側の空燃比で運転するエンジン始動直後、同じく大きな出力が要求されるために理論空燃比の空燃比よりもリッチ側の空燃比で運転する全負荷運転時）にも、燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬ、オーバーラップ中吹き返しガス量ＭＲＥＳＯＬ、これらの合計である内部不活性ガス量ＭＲＥＳ、これに基づく内部不活性ガス率ＭＲＥＳＦＲを精度良く算出できる。

また、オーバーラップ期間の積算有効面積ＡＳＵＭＯＬを仮想オリフィスの面積とし、この仮想オリフィスを排気が燃焼室５から吸気系へと吹き抜けると仮定しているので、オーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬの算出が簡略化されている。

次に、ノッキングについて説明する。図２９はノック発生時の燃焼室５内の圧力履歴を示している。高周波分を取り除いた平均圧力を改めて描いてみると、自着火時期θｋｎｋに燃焼室５内の圧力が一気に上昇していることがわかる。このノックによる圧力上昇は燃焼室５内の未燃混合気が等容燃焼した結果で生じるとみなして、その圧力上昇分ＤＰを次のように熱力学の計算式により算出する。

いま、未燃燃料量ＭＵＢの未燃ガスが等容燃焼で全て燃えると仮定すると、熱力学より発熱量Ｑは次式で与えられる。

Ｑ＝ＣＦ＃×ＭＵＢ …（補４）
ただし、ＣＦ＃：燃料の低位発熱量、
一方、この発熱量Ｑにより燃焼室５内のガス温度が上昇するので、この温度上昇分をΔＴとすれば次式が成立する。

Ｑ＝Ｃｖ×Ｍ×ΔＴ …（補５）
ただし、Ｍ：燃焼室５内の全てのガスの質量、
Ｃｖ：既燃ガスの定容比熱、
（補４）、（補５）の両式は等しいと置いて温度上昇分ΔＴについて解くと次式が得られる。

ΔＴ＝（ＣＦ＃×ＭＵＢ）／（Ｃｖ×Ｍ） …（補６）
気体の状態方程式ＰＶ＝ｎＲＴの両辺を微分する（ただし、定容変化であるためＶは一定）。

Ｖ×ｄＰ＝ｄｎ×Ｒ×Ｔ＋ｎ×Ｒ×ｄＴ …（補７）
ここで、ノックが発生するような状態ではモル数ｎの変化は小さいため、（補７）式右辺のｄｎ＝０として次式を得る。

ｄＰ＝（ｎ×Ｒ／Ｖ）×ｄＴ …（補８）
（補８）、（補６）の両式より温度上昇分ｄＴ（＝ΔＴ）を消去し圧力上昇分ＤＰについて整理すると最終的に次式を得る。

ｄＰ＝ｎ×Ｒ×ＣＦ＃×ＭＵＢ／（Ｖ×Ｃｖ×Ｍ）…（補９）
すなわち、（補９）式は、未燃燃料量ＭＵＢ、燃焼室５の自着火時期における容積Ｖ、既燃ガスの定容比熱Ｃｖ、燃焼室５内の全てのガスの質量Ｍ、燃焼室５内の全てのガスの総モル数ｎが分かれば、圧力上昇分ＤＰを計算式により求めることができることを示している。

この場合、燃焼室５の自着火時期は公知の手法を用いて求めることができる。公知の手法とは、単位クランク角毎に燃焼室５内の温度及び圧力を算出して図３０よりその温度、圧力に対する１／τを求め、この１／τを積算した値が１となるときのクランク角を自着火時期θｋｎｋとするものである。ここで、図３０のτは燃焼室５内の燃料が自着火にいたるまでの時間である。自着火時期θｋｎｋが求まれば、このθｋｎｋから燃焼室５の自着火時期における容積Ｖｋｎｋが機械的に求まる。

一方、自着火時期θｋｎｋが分かれば図３１より自着火時の燃焼質量割合ＢＲｋｎｋを求めることができ、この自着火時の燃焼質量割合ＢＲｋｎｋと燃料量ＱＩＮＪとから次の計算式により未燃燃料量ＭＵＢを求めることができる。すなわち、燃料量ＱＩＮＪのうちＢＲｋｎｋが既に燃焼している割合であるから、残りの１−ＢＲｋｎｋがまだ燃焼していない割合となる。

ＭＵＢ＝ＱＩＮＪ×（１−ＢＲｋｎｋ） …（補１０）
ただし、図３１では計算を簡単にするため、燃焼開始遅れ期間、初期燃焼期間、主燃焼期間の３つ分け、各期間の特性を直線で近似している。

この場合、図３１の特性はエンジンの負荷や回転速度が相違しても変わらないので、エンジン仕様が定まれば一義的に定まることから、運転条件の相違毎に適合する必要はない。

次に、既燃ガスの定容比熱Ｃｖについても次のように熱力学の計算式により算出することができる。すなわち、定圧比熱Ｃｐの定義はＣｐ＝（∂Ｅ／∂Ｔ）ｐであるから、この式を積分すると次式を得る。

∫ｄＥ＝Ｃｐ×∫ｄＴ …（補１１）
∴Ｅ＝Ｃｐ×Ｔ …（補１２）
（補１２）式より定圧比熱Ｃｐは次式で与えられる。

Ｃｐ＝Ｅ／Ｔ …（補１３）
一方、理想気体で等圧変化のときにはＣｐ−Ｃｖ＝Ｒが成立するので、この式と（補１３）式より定圧比熱Ｃｐを消去し定容比熱Ｃｖについて整理すると最終的に次式が得られる。

Ｃｖ＝Ｅ／Ｔ−Ｒ …（補１４）
ただし、Ｅ：エンタルピ、
Ｔ：燃焼室５内の自着火時時期における平均温度、
上記（補９）式の燃焼室５内の全てのガスの質量Ｍは次式により計算できる。

Ｍ＝ＭＲＥＳ＋ＭＡＣＹＬ＋ＱＩＮＪ …（補１５）
ただし、ＭＲＥＳ：内部不活性ガス量、
ＭＡＣＹＬ：シリンダ新気量、
ＱＩＮＪ：燃料量、
このように、未燃燃料量ＭＵＢ、既燃ガスの定容比熱Ｃｖ、燃焼室５内の全てのガスの質量Ｍも（補１０）、（補１４）、（補１５）の計算式によりそれぞれ求めることができることがわかる。残る未知数は、上記（補９）式の燃焼室５内の全てのガスの総モル数ｎ、（補１４）式のエンタルピＥ及び燃焼室５内の自着火時期における平均温度Ｔ（＝ＴＥ）である。

ここで、上記（補９）式の燃焼室５内の全てのガスの総モル数ｎ及び各成分ガスのモル数は燃焼の基礎式を用いて計算により求めることができ、各成分ガスのモル数と実験式とを用いて（補１４）式のエンタルピＥを計算することができる。燃焼室５内の自着火時期における平均温度ＴＥも熱力学の式を用いて計算により求めることができる。

このように、ノックによる圧力上昇分ＤＰを、テーブルやマップに頼らずにほぼ計算式により求められるように構成しているのであり、これによりテーブルやマップ作成に要する実験の工数と時間を大幅に削減することができる。

後は、このようにして得た圧力上昇分ＤＰとノックとを関連させることであり、これをノック強度推定値に換算する。

次に、エンジンコントローラ３１で実行される上記ノックによる圧力上昇分ＤＰ及びノック強度推定値ＫＩＣの算出を以下のフローチャートを参照しながら詳述する。

図３２は後述するノック強度推定値ＫＩＣの算出に必要となる物理量を予め算出するためのもので、一定時間毎（例えば１０ｍｓ毎）に実行する。図１３のフローに続けて実行する。なお、以下では前述のフローで既に求めている物理量を重複して求めているような部分もあるが、説明する。

ステップ１６１では、図１４のステップ５２、５３で算出されているシリンダ新気量ＭＡＣＹＬ［ｇ］及び内部不活性ガス量ＭＲＥＳ［ｇ］、燃料量ＱＩＮＪ［ｇ］、図１０のステップ２９で算出されている初期燃焼期間ＢＵＲＮ１［ｄｅｇ］、図１２のステップ３４で算出されている主燃焼期間ＢＵＲＮ２［ｄｅｇ］、図１３のステップ４３で算出されている基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］のほか、ＩＧＮＤＥＡＤ［ｄｅｇ］、クランク角Θ［ｄｅｇＡＴＤＣ］を読み込む。なお、燃料量ＱＩＮＪ［ｇ］は燃料噴射パルス幅ＴＩ［ｍｓ］に比例させて求めればよい。

上記のクランク角Θは後述するように圧縮上死点ＴＤＣを基準のゼロとしてこれより遅角側をプラス、これより進角側をマイナスに採った値であり、クランク角センサ（３３、３４）により検出されるクランク角θとは相違する。ただし、２つのクランク角θとΘの間には一定の関係があるので、クランク角θからクランク角Θを求めている。

ステップ１６２ではシリンダ新気量ＭＡＣＹＬ［ｇ］をＷＩＤＲＹ［ｇ］に、内部不活性ガス量ＭＲＥＳ［ｇ］をＭＡＳＳＺ［ｇ］に移す。これらＷＩＤＲＹ、ＭＡＳＳＺは、ノック強度指標ＫＮＫＩの算出にのみ用いるために導入したもので、ＷＩＤＲＹはシリンダ新気量、ＭＡＳＳＺは内部不活性ガス量である。

ステップ１６３では燃焼室５の内部不活性ガス率ＲＴＯＥＧＲを次式により算出する。

ＲＴＯＥＧＲ＝ＭＡＳＳＺ／（ＭＡＳＳＺ＋ＷＩＤＲＹ＋ＱＩＮＪ）…（５９）
ステップ１６４ではクランク角毎の燃焼質量割合ＢＲを算出する。ここでは、図３１に示したように、点火から燃焼開始までの期間、初期燃焼期間、主燃焼期間の３つの期間をそれぞれ直線で近似する。この場合、変数としてのクランク角Θは圧縮上死点ＴＤＣを基準のゼロとしてこれより遅角側をプラス、これより進角側をマイナスに採った値であり、このクランク角Θ［ｄｅｇＡＴＤＣ］を用いると燃焼質量割合ＢＲは次のような一次式となる。

燃焼遅れ期間；
ＢＲ＝０ …（６０）
初期燃焼期間；
ＢＲ＝ＳＳ１×（Θ＋ＭＢＴＣＡＬ−ＩＧＮＤＥＡＤ） …（６１）
主燃焼期間 ；
ＢＲ＝０．０２＋ＳＳ２×（Θ＋ＭＢＴＣＡＬ−ＩＧＮＤＥＡＤ−ＢＵＲＮ１）
…（６２）
ただし、ＳＳ１：０．０２／ＢＵＲＮ１、
ＳＳ２：０．５８／ＢＵＲＮ２、
ステップ１６５では燃焼室５内の燃料が全て燃焼した（つまりＢＲ＝１）ときの各ガス成分のモル数を算出する。ただし、ガス成分としては、燃料のほかは、Ｏ₂、Ｎ₂、ＣＯ₂、ＣＯ、Ｈ₂Ｏに限定する。また、ガソリンの燃料組成をＣ₇Ｈ₁₄で近似している。

まず、燃料量ＱＩＮＪ［ｇ］の燃料が燃焼した後に発生する総排出ガスのモル数ＷＥＤＲＹ［ｍｏｌ］と排出ガス中のＯ₂、Ｎ₂、ＣＯ₂、ＣＯ、Ｈ₂Ｏといったそれぞれのガス成分のモル数ＸＥＯ２［ｍｏｌ］、ＸＥＮ２［ｍｏｌ］、ＸＥＣＯ２［ｍｏｌ］、ＸＥＣＯ［ｍｏｌ］、ＸＥＨ２Ｏ［ｍｏｌ］を次のように算出する。

総排出ガス；ＷＥＤＲＹ＝ＭＩＤＲＹ＃×ＷｌＤＲＹ−ＱｌＮＪ
／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ）×（Ａ＃／４） …（６３．１）
酸素 ；ＸＥＯ２ ＝｛ＭＩＤＲＹ＃×ＷｌＤＲＹ×０．２１−ＱＩＮＪ
／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ）×（Ｂ＃＋Ａ＃／４）｝
／ＷＥＤＲＹ …（６３．２）
二酸化炭素；ＸＥＣＯ２＝｛ＱＩＮＪ／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃）×Ｂ＃｝
／ＷＥＤＲＹ …（６３．３）
一酸化炭素；ＸＥＣＯ ＝０ …（６３．４）
窒素 ；ＸＥＮ２ ＝１−ＸＥＯ２−ＸＥＣＯ２−ＸＥＣＯ …（６３．５）
水 ；ＸＥＨ２Ｏ＝｛ＭＩＤＲＹ＃×ＷＩＤＲＹ×１５／７４５
＋ＱＩＮＪ／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃）
×Ａ＃／２｝／ＷＥＤＲＹ …（６３．６）
ただし、ＭＩＤＲＹ＃：１ｇ当たりの新気ガスのモル数
ＡＨ＃：水素のモル質量、
ＡＣ＃：炭素のモル質量、
Ａ＃、Ｂ＃：定数、
ここで、ガソリンの燃料組成をＣ₇Ｈ₁₄で近似しているため、定数Ａ＃は１４、定数Ｂ＃は７である。

次に、燃焼サイクル初期の各ガス成分のモル数ＷＧＡＳ［ｍｏｌ］、ＷＥＧＲ［ｍｏｌ］、ＷＯ２［ｍｏｌ］、ＷＮ２［ｍｏｌ］、ＷＣＯ２［ｍｏｌ］、ＷＣＯ［ｍｏｌ］、ＷＨ２Ｏ［ｍｏｌ］を次のように算出する。

燃料 ；ＷＧＡＳ＝ＱＩＮＪ／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃）
…（６４．１）
不活性ガス；ＷＥＧＲ＝ＭＩＤＲＹ＃×ＷＩＤＲＹ×ＲＴＯＥＧＲ…（６４．２）
酸素 ；ＷＯ２ ＝ＭＩＤＲＹ＃×ＷＩＤＲＹ×０．２１＋ＷＥＧＲ×ＸＥＯ２
…（６４．３）
窒素 ；ＷＮ２ ＝ＭＩＤＲＹ＃×ＷＩＤＲＹ×０．８９＋ＷＥＧＲ×ＸＥＮ２
…（６４．４）
二酸化炭素；ＷＣＯ２＝ＷＥＧＲ×ＸＥＣＯ２ …（６４．５）
一酸化炭素；ＷＣＯ ＝ＷＥＧＲ×ＸＥＣＯ …（６４．６）
水 ；ＷＨ２Ｏ＝ＭＩＤＲＹ＃×ＷＩＤＲＹ×１５／７４５
＋ＷＥＧＲ×ＸＥＨ２Ｏ …（６４．７）
次に、全て燃焼した（つまりＢＲ＝１）ときの各ガス成分のモル数ＭＬＧＡＳ［ｍｏｌ］、ＭＬＯ２［ｍｏｌ］、ＭＬＮ２［ｍｏｌ］、ＭＬＣＯ２［ｍｏｌ］、ＭＬＣＯ［ｍｏｌ］、ＭＬＨ２Ｏ［ｍｏｌ］を次のように算出する。

燃料 ；ＭＬＧＡＳ＝ＷＧＡＳ−ＱＩＮＪ／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃）
…（６５．１）
酸素 ；ＭＬＯ２ ＝ＷＯ２−（Ｂ＃＋Ａ＃／４）×ＱｌＮＪ
／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃） …（６５．２）
窒素 ；ＭＬＮ２ ＝ＷＮ２ …（６５．３）
二酸化炭素；ＭＬＣＯ２＝ＷＣＯ２＋Ｂ＃×ＱＩＮＪ
／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃） …（６５．４）
一酸化炭素；ＭＬＣＯ ＝ＷＣＯ …（６５．５）
水 ；ＭＬＨ２Ｏ＝ＷＨ２Ｏ＋Ａ＃／２×ＱＩＮＪ
／（Ｂ＃×ＡＣ＃＋Ａ＃×ＡＨ＃） …（６５．６）
これで、燃焼室５内の燃料が全て燃焼した（つまりＢＲ＝１）ときの各ガス成分のモル数の算出を終了するので、ステップ１６６に進み各ガス成分のモル数の総和を燃焼室５内の燃料が全て燃焼したときの総ガスモル数ＭＬＡＬＬとして、つまり次式により総ガスモル数ＭＬＡＬＬを算出する。

ＭＬＡＬＬ＝ＭＬＧＡＳ＋ＭＬＯ２＋ＭＬＮ２＋ＭＬＣＯ２＋ＭＬＣＯ＋ＭＬＨ２Ｏ
…（６５．７）
図３３は一燃焼サイクル毎に自着火時期θｋｎｋ［ｄｅｇＡＴＤＣ］を算出するためのもので、所定のクランク角毎（例えば実機については１０ｄｅｇ程度毎）に実行する。図３２のフローに続けて実行する。

自着火時期θｋｎｋは次のようにして求めることができる。すなわち、ある温度Ｔ及びある圧力のときに燃焼室５内の燃料が自着火にいたるまでの時間をτとしたとき、その逆数である１／τは図３０に示した分布となり、クランク角θについてのその積分値∫１／τｄθ＝１（積分範囲は吸気弁閉時期ＩＶＣから所定のクランク角θ１まで）が成立するときのクランク角θ１が自着火時期θｋｎｋになることが知られているので、これを利用する。

図３３においてまず、ステップ１７１ではクランク角センサ（３３、３４）により検出される現在のクランク角θが吸気弁閉時期ＩＶＣ以降（吸気弁閉時期を含む）か否かをみる。このクランク角θは、燃焼質量割合ＢＲの算出に用いるクランク角Θと相違して、基準の位置は採りやすい上死点や下死点など任意の位置に設ければよい。

現在のクランク角θが吸気弁閉時期ＩＶＣ以降にあればステップ１７２に進み図３２のステップ１６２、１６４で算出されている内部不活性ガス量ＭＡＳＳＺ、シリンダ新気量ＷＩＤＲＹ、燃焼質量割合ＢＲのほか、燃料量ＱＩＮＪ、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、温度センサ４５より検出される排気温度ＴＥＸＨ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］を読み込む。

ステップ１７３では燃焼室５の圧縮開始時温度ＴＣ０［Ｋ］を次式により算出する。

ＴＣ０＝｛（ＷＩＤＲＹ＋ＱＩＮＪ）×ＴＣＯＬ＋ＭＡＳＳＺ×ＴＥＸＨ｝
／（ＷＩＤＲＹ＋ＱＩＮＪ＋ＭＡＳＳＺ） …（６６）
ここでは、不活性ガスと新気の比熱を等しくして式を簡略化している。

ステップ１７４では燃焼室５の圧縮開始時圧力ＰＣ０［Ｐａ］を算出する。これは圧力センサ４４により検出される吸気弁閉時期ＩＶＣのコレクタ内圧力ＰＣＯＬをＰＣ０とすればよい。

ステップ１７５では燃焼室５内の瞬間圧縮比εθを算出する。この瞬間圧縮比εθは燃焼室５の隙間容積Ｖｃを燃焼室５の現在のクランク角θにおける容積で除算した値である。燃焼室５の現在のクランク角θにおける容積はピストン６のストローク位置つまりエンジンのクランク角により定まるので、クランク角をパラメータとするテーブルを予め作成しておき現在のクランク角θからこのテーブルを検索することにより求めればよい。

ステップ１７６、１７７では燃焼室５内の燃料が燃焼したときの平均温度ＴＣ［Ｋ］と平均圧力ＰＣ［Ｐａ］を次式により算出する。

ＴＣ＝ＴＣ０×εθ＾０・３５
＋ＣＦ＃×ＱＩＮＪ×ＢＲ／（ＭＡＳＳＺ＋ＷＩＤＲＹ＋ＱＩＮＪ）
…（６７）
ＰＣ＝ＰＣ０×εθ＾１．３５×ＴＣ／ＴＣ０／εθ＾０．３５…（６８）
ただし、εθ ：瞬間圧縮比、
ＣＦ＃：燃料の低位発熱量、
（６７）、（６８）式は燃焼室５内でガスが断熱圧縮されると共に定容変化で燃焼すると仮定したときの式である。すなわち、（６７）式右辺第１項が断熱圧縮後の温度を、（６８）式右辺のＰＣ０×εθ＾１．３５が断熱圧縮後の圧力を、これに対して（６７）式右辺第２項が定容変化で燃焼により温度上昇した分を、（６８）式右辺のＴＣ／ＴＣ０／εθ＾０．３５が定容変化での燃焼による圧力上昇率を表している。

ステップ１７８では燃焼室５内の未燃混合気の温度Ｔを次式により算出する。

Ｔｕｂ＝ＴＣ０×εθ＾０・３５×（ＰＣ／ＰＣ０／εθ＾１．３５）
＾（０．３５／１．３５） …（６９）
（６９）式は燃焼室５内でガスが断熱圧縮されると共に、（６７）式とは相違して、可逆断熱変化で燃焼すると仮定したときの式である。すなわち、（６９）式右辺のＴＣ０×εθ＾０・３５が断熱圧縮後の温度を、（６９）式右辺の（ＰＣ／ＰＣ０／εθ＾１．３５）＾（０．３５／１．３５）が可逆断熱変化での燃焼による温度上昇率を表している。

なお、未燃混合気の圧力は上記（６８）式の平均圧力ＰＣと等しいと仮定する。

ここで、（６７）式の平均温度ＴＣと、（６９）式の未燃混合気の温度Ｔｕｂとの違いは次のようなものである。すなわち、（６７）式の平均温度ＴＣは燃焼室５内で発生した熱が燃焼室５内の全てのガスを昇温させると仮定したときの温度である。これに対して、（６９）式の未燃混合気温度Ｔｕｂは燃焼室５内でガスが既燃ガスと未燃ガスの２つに分離された状態にあり、燃焼室５内で発生した熱は既燃ガスのみを昇温させると仮定したときの温度である。そして、未燃混合気の自着火により急激な圧力上昇が生じてノックが発生する。

ステップ１７９では未燃混合気温度Ｔｕｂと未燃混合気圧力（＝ＰＣ）から図３０を内容とするマップを検索することにより１／τを算出する。１／τは図３０のように温度、圧力が大きくなるほど大きくなる値である。この１／τの特性は燃料の種類（例えばハイオクタンガソリンとレギュラーガソリン）により異なるので、使用する燃料に合わせて選択すればよい。

ステップ１８０ではこの１／τをＳＵＭに加算する。ＳＵＭは１／τの積算値を表す。この積算値ＳＵＭには吸気弁閉時期ＩＶＣになったときに初期値のゼロを入れている。

ステップ１８１ではこの積算値ＳＵＭと１を比較する。吸気弁閉時期ＩＶＣの直後には積算値ＳＵＭは１に満たないので、そのまま今回の処理を終了する。次回より所定のクランク角毎ここでは１０ｄｅｇ毎にステップ１７２〜１８０の操作を繰り返すと、やがて積算値ＳＵＭ＝１となるので、このときには自着火時期になったと判断し、ステップ１８１よりステップ１８２に進んでそのときのクランク角θを自着火時期θｋｎｋに移す。

ステップ１８３では、自着火時の燃焼質量割合ＢＲｋｎｋを算出する。これは自着火時期θｋｎｋが初期燃焼期間にあるときにはこの自着火時期θｋｎｋを圧縮上死点ＴＤＣを基準とする上記のクランク角Θに変換し、その変換したクランク角Θを上記の（６１）式に、また自着火時期θｋｎｋが主燃焼期間にあるときにはこの自着火時期θｋｎｋを圧縮上死点ＴＤＣを基準とする上記のクランク角Θに変換し、その変換したクランク角Θを上記の（６２）式にそれぞれ代入することによって求めることができる。

ステップ１８４では、燃焼室５の自着火時θｋｎｋにおける平均温度ＴＥを算出する。これは上記（６７）式右辺の燃焼質量割合ＢＲに１．０を代入して得られる燃焼室５の平均温度ＴＣを燃焼室５の自着火時平均温度ＴＥとして求めればよい。

ステップ１８５では、燃焼室５の自着火時θｋｎｋにおける容積Ｖｋｎｋを算出する。燃焼室５の自着火時期θｋｎｋにおける容積Ｖｋｎｋは、燃焼室５の現在のクランク角θにおける容積と同様に、ピストン６のストローク位置つまりエンジンのクランク角により定まるので、クランク角をパラメータとするテーブルを予め作成しておき自着火時θｋｎｋからこのテーブルを検索することにより求めればよい。

ステップ１８６ではノック強度推定値ＫＩＣを算出する。このノック強度推定値ＫＩＣの算出については図３４のフローにより説明する。図３４（図３３ステップ１８６のサブルーチン）においてステップ１９１では図３２のステップ１６２、１６３、１６６で算出されている内部不活性ガス量ＭＡＳＳＺ［ｇ］、シリンダ新気量ＷＩＤＲＹ［ｇ］、総ガスモル数ＭＬＡＬＬ［ｍｏｌ］、図３３のステップ１８３、１８４、１８５で算出されている自着火時の燃焼質量割合ＢＲｋｎｋ、自着火時平均温度ＴＥ、自着火時容積Ｖｋｎｋのほか、燃料量ＱＩＮＪ［ｇ］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］を読み込む。

ステップ１９２では燃料量ＱＩＮＪ［ｇ］と自着火時の燃焼質量割合ＢＲｋｎｋとから自着火時の未燃燃料量ＭＵＢ［ｇ］を次式により算出する。

ＭＵＢ＝ＱＩＮＪ×（１−ＢＲｋｎｋ） …（７０）
（７０）式は上記（補１０）式そのものである。

ステップ１９３ではガスエンタルピ（自着火する燃料ガスのエンタルピ）Ｅ［ｃａｌ／ｍｏｌ］を算出する。このガスエンタルピの算出については図３５のフローにより説明する。図３５（図３４ステップ１９３のサブルーチン）においてステップ２０１では、図３３のステップ１８４で算出されている燃焼室５の自着火時平均温度ＴＥ、図３２のステップ１６５、１６６で算出されている各ガス成分のモル数ＭＬＧＡＳ、ＭＬＯ２、ＭＬＮ２、ＭＬＣＯ２、ＭＬＣＯ、ＭＬＨ２Ｏ、総ガスモル数ＭＬＡＬＬを読み込む。

ステップ２０２では自着火時平均温度ＴＥから各ガス成分のエンタルピＥＯ２、ＥＮ２、ＥＣＯ２、ＥＣＯ、ＥＨ２Ｏを算出する。各ガス成分のエンタルピは次の水谷の実験式（内燃機関ｖｏｌ．１１ Ｎｏ．１２５ｐ７９参照）を用いて算出すればよい。
（１）ＴＥ＜１２００Ｋの場合
Ｅ＝Ａ０＃＋１０００×（Ａ１＃×（ＴＥ／１０００）
＋Ａ２＃／２×（ＴＥ／１０００）＾２
＋Ａ３＃／３×（ＴＥ／１０００）＾３
＋Ａ４＃／４×（ＴＥ／１０００）＾４
＋Ａ５＃／５×（ＴＥ／１０００）＾５）＋ＨＤＬ＃…（７１．１）
（２）ＴＥ＞１２００Ｋの場合
Ｅ＝Ｂ０＃＋１０００×（Ｂ１＃×（ＴＥ／１０００）
＋Ｂ２＃×ＬＮ（ＴＥ／１０００）
−Ｂ３＃／（ＴＥ／１０００）
−Ｂ４＃／２／（ＴＥ／１０００）＾２
−Ｂ５＃／３／（ＴＥ／１０００）＾３）＋ＨＤＬ＃…（７１．２）
ただし、Ａ０＃〜Ａ５＃、Ｂ０〜Ｂ５＃、ＨＤＬ＃は実験により求める適合値、
ステップ２０３では燃料のエンタルピＥＧを次式により算出する。

ＥＧ＝Ｂ＃／ＡＣ＃×ＥＣＯ２＋Ａ＃／ＡＨ＃×ＥＨ２０／２
＋（Ｂ＃／ＡＣ＃＋Ａ＃／ＡＨ＃／４）×ＥＯ２…（７１．３）
ステップ２０４では各ガス成分の平均エンタルピＥを次式により計算して図３５の処理を終了し、図３４のステップ１９４に戻る。

Ｅ＝（ＭＬＧＡＳ×ＥＧ＋ＭＬＯ２×ＥＯ２＋ＭＬＮ２×ＥＮ２
＋ＭＬＣＯ２×ＥＣＯ２＋ＭＬＣＯ×ＥＣＯ＋ＭＬＨ２Ｏ×ＥＨ２Ｏ）
／ＭＬＡＬＬ …（７１．４）
図３４のステップ１９４では既燃ガスの定容比熱Ｃｖ［Ｊ／Ｋ・ｇ］をガスエンタルピＥ、燃焼室５の自着火時平均温度ＴＥを用いて次式により計算する。

Ｃｖ＝Ｅ／ＴＥ−Ｒ＃…（７２）
ただし、Ｒ＃：普遍気体定数、
（７２）式は上記（補１４）式においてＴ→ＴＥ、Ｒ→Ｒ＃の置き換えにより得られる式である。

ステップ１９５では自着火による圧力上昇分つまりノックによる圧力上昇分ＤＰ［Ｐａ］を次式により算出する。

ＤＰ＝（ＷＡＬＬ×ＭＵＢ×Ｒ＃×ＣＦ＃）
／｛Ｃｖ×Ｖｋｎｋ×（ＭＡＳＳＺ＋ＱＩＮＪ＋ＷＩＤＲＹ）｝
…（７３）
ただし、ＣＦ＃：燃料の低位発熱量、
（７３）式の圧力上昇分ＤＰは図２９に示したようにノック発生により燃焼室５内の圧力がステップ的に大きくなるので、この圧力上昇分を計算式により求めるようにしたものである。

（７３）式は上記（補９）式においてｄＰ→ＤＰ、ｎ→ＭＬＡＬＬ、Ｒ→Ｒ＃、Ｖ→Ｖｋｎｋ、Ｍ→ＭＡＳＳＺ＋ＷＩＤＲＹ＋ＱＩＮＪの置き換えにより得られる式である。

ステップ１９６ではノック強度推定値基本値ＫＩＣ０を次式により算出する。

ＫＩＣ０＝相関係数１×ＤＰ …（７４）
ここで、（７４）式右辺の相関係数１はノック強度との相関を表すための係数である。この場合、ノックによる圧力上昇分ＤＰが大きいほどノック強度推定値基本値ＫＩＣ０が大きくなるようにしている。

ステップ１９７ではエンジン回転速度ＮＲＰＭから図３６を内容とするテーブルを検索することにより回転速度補正係数ＫＮを次式により算出し、ステップ１９８でこの回転速度補正係数ＫＮをノック強度推定値基本値ＫＩＣ０に乗算した値をノック強度推定値ＫＩＣとして、つまり次式によりノック強度推定値指標ＫＩＣを算出する。

ＫＩＣ＝ＫＩＣ０×ＫＮ …（７５）
ここで、回転速度補正係数ＫＮは、エンジン回転速度ＮＲＰＭが低いときのほうが回転速度が高いときよりノックによる圧力振動をドライバーが強く感じるので、この違いをノック強度推定値に反映させるためのものである。すなわち、図３６のようにＫＮの値は、基準回転速度ＮＲＰＭ０のときを１．０としてこれより低い回転速度域では１．０を超える値、この逆に基準回転速度ＮＲＰＭ０より高い回転速度域で１．０未満の値である。実際の値は実験により適合する。

このようにしてノック強度推定値ＫＩＣの算出を終了したら図３３に戻り図３３の処理をすべて終了する。その後は、現在のクランク角θが次の燃焼サイクルの吸気弁閉時期ＩＶＣの近くにくるまでそのまま待機し、次の燃焼サイクルの吸気弁閉時期ＩＶＣの近くにくると再び図３３の処理を実行する。こうして一燃焼サイクル毎に自着火時期θｋｎｋ、ノック強度推定値ＫＩＣが求められる。

図３７は燃焼室５内のデポジット量ＶＤＥＰを算出するためのもので、一定時間毎（例えば１０ｍｓ毎）に実行する。

ステップ２２１ではシリンダブロックに取り付けられたノックセンサ４７（実ノック強度検出手段）からの電圧信号に基づいてノックが発生しているか否かをみる。ノックは特定の周波数（５ｋ〜２０ｋＨｚ）で振動する現象として現れるので、ノックセンサ４７はこの周波数帯の振動を電圧信号に変えて出力する。この電圧信号がある基準値よりも大きくなったときにノックが発生したと判断しステップ２２２に進んでノック強度ＫＩＤを検出する。ノック強度が大きいほどノックセンサ４７からの出力が大きくなるのでその値を用いて次式により実ノック強度ＫＩＤを検出する。

ＫＩＤ＝相関係数２×ノックセンサ出力 …（７６）
ここで、（７６）式右辺の相関係数２は、（７４）式の相関係数１と同じにノック強度との相関を表すための係数で、ノック強度が大きいほど実ノック強度ＫＩＤも大きくなる。

ステップ２２３では図３３のステップ１８６で算出されているノック強度推定値ＫＩＣを読み込み、このノック強度推定値ＫＩＣと実ノック強度ＫＩＤをステップ２２４において比較する。ノック強度推定値ＫＩＣは推定によりノック強度を算出した値で、本来ノックセンサ４７で検出された値と等しくなるはずである。しかしながら、燃焼室５内にデポジットが堆積し燃焼室５の容積が減少するとノック発生時期が早くなり、実ノック強度も強くなる。ノック強度推定値ＫＩＣはあくまで燃焼室５内にデポジットが堆積していない状態での燃焼室容積（Ｖｋｎｋ）に基づいて計算されているので、燃焼室５内にデポジットが堆積している場合には実ノック強度ＫＩＤよりも強いノック強度推定値ＫＩＣを算出することになりＫＩＤ≠ＫＩＣとなる。このときには、ステップ２２４よりステップ２２５に進んで実ノック強度ＫＩＤとノック強度推定値ＫＩＣの差ＫＩＳを計算し、この差ＫＩＳから図３８を内容とするテーブルを検索することにより燃焼室５内に堆積しているデポジット量ＶＤＥＰを算出する。

ここで、デポジット量ＶＤＥＰは燃焼室５内でのデポジット堆積による燃焼室５の容積増加分を表し、図３８に示したように差ＫＩＳが大きくなるほど大きくなる値である。

一方、ＫＩＤ＝ＫＩＣであるときには燃焼室５内にデポジットが堆積していないと判断しそのまま今回の処理を終了する。

このようにして算出したデポジット量ＶＤＥＰは燃焼開始時容積Ｖ０の補正に用いる。すなわち、図１０においてステップ１６３、１７１では燃焼室５の燃焼開始時における容積Ｖ０からこのデポジット量ＶＤＥＰだけ減算する。

次に、図３９は点火時期指令値ＱＡＤＶを算出するためのもので、一定時間毎（例えば１０ｍｓ毎）に実行する。図１３のフローに続けて実行する。

ステップ２１１では、図３３のステップ１８６により算出されているノック強度推定値ＫＩＣ、図１３のステップ４３により算出されている基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を読み込む。

ステップ２１２ではノック強度推定値ＫＩＣと所定のスライスレベルＳＬを比較する。
ここで用いられるノック強度推定値ＫＩＣは前回の燃焼サイクルでの値である。ノック強度推定値ＫＩＣがスライスレベルＳＬより大きいときには前回の燃焼サイクルにおいてノックが発生したことを意味するので、今回の燃焼サイクルではノックを回避するため、ステップ２１３に進み基本点火時期ＭＢＴＣＡＬより所定値（遅角量）［ｄｅｇ］遅らせた値を点火時期指令値ＱＡＤＶ［ｄｅｇＢＴＤＣ］とする。これに対して、ノック強度推定値ＫＩＣがスライスレベルＳＬ以下のときには前回の燃焼サイクルにおいてノックが発生していないことを意味し、従って今回の燃焼サイクルにおいてノックを回避する必要がないのでステップ２１２よりステップ２１４に進み基本点火時期ＭＢＴＣＡＬをそのまま点火時期指令値ＱＡＤＶとする。

ここで、本実施形態の作用効果を説明する。

自着火時期（ノック発生時期）に基づいて燃焼室５の自着火時期における容積Ｖｋｎｋ（ノック発生時期の燃焼室容積）を算出することができ、自着火時の燃焼質量割合ＢＲｋｎｋをこの燃焼室５の自着火時期における容積Ｖｋｎｋで割った値（ＢＲｋｎｋ／Ｖｋｎｋ）は、ノック強度（つまりＤＰ）と相関がある。すなわち、上記（７０）、（７３）式よりＤＰは（１−ＢＲｋｎｋ）／Ｖｋｎｋに比例するので、ＤＰとＢＲｋｎｋ／Ｖｋｎｋとの間には相関がある。そこで、（１−ＢＲｋｎｋ）／Ｖｋｎｋを用いて自着火による圧力上昇分ＤＰ、最終的にはノック強度推定値を算出し、このノック強度推定値と実ノック強度を比べることで燃焼室５内のデポジット量を算出できる。

本実施形態（請求項４に記載の発明）はこうした考え方に基づき、図３１に示したように燃焼質量割合ＢＲを算出し（図３２のステップ１６４）、自着火時期θｋｎｋ（ノック発生時期）を算出（推定）し（図３３のステップ１８２）、この自着火時期θｋｎｋと燃焼質量割合ＢＲから自着火時の燃焼質量割合ＢＲｋｎｋを算出し（図３３のステップ１８３）、自着火時期θｋｎｋに基づいて燃焼室５の自着火時期における容積Ｖｋｎｋ（ノック発生時期の燃焼室容積）を算出し（図３３のステップ１８５）、この燃焼室５の自着火時期における容積Ｖｋｎｋと自着火時の燃焼質量割合ＢＲｋｎｋ（ノック発生時期の燃焼質量割合）に基づいてノック強度推定値ＫＩＣを算出し（図３４のステップ１９２、１９５〜１９８）、このノック強度推定値ＫＩＣと実ノック強度ＫＩＤに基づいて燃焼室５内のデポジット量ＶＤＥＰを算出する（図３７のステップ２２５、２２６）ので、吸気系にはデポジットが堆積していても燃焼室５内にはデポジットが堆積していない場合や吸気系にはデポジットが堆積していないのに燃焼室５内にはデポジットが堆積している場合においても、燃焼室５内のデポジット量を精度よく算出できる。

また、デポジット量の算出結果を情報として知らせるようにすることで、サービス工場における清掃等の対応が可能となり、エンジンの効率低下や排気悪化等を事前に防ぐことができる。

一方、燃焼室５内にデポジットが堆積すると圧縮比が上がりノッキングが発生しやすくなる。この場合に本実施形態（請求項９に記載の発明）によれば算出したデポジット量に基づいてノックが生じないように点火時期や吸気弁開時期を遅角側に制御するので、燃焼室内のデポジットの堆積に伴うノッキングの発生を防止できる。

ノック強度推定値ＫＩＣはデポジット量ＶＤＥＰが多いほど実ノック強度ＫＩＤに比べ、大きくなる傾向にある。したがってＫＩＤとＫＩＣの差ＫＩＳが大きいほど燃焼室５内のデポジット量が多いといえる。本実施形態（請求項５に記載の発明）によれば、その差ＫＩＳを用いることで容易に燃焼室内のデポジット量を算出できる（図３８）。

本実施形態（請求項７に記載の発明）によれば、自着火時期θｋｎｋを燃焼質量割合ＢＲに基づいて推定するので（図３３のステップ１７２〜１８２）、筒内圧センサなどノック発生時期を検出するセンサを使わなくても自着火時期θｋｎｋを求めることができる。

ただし、筒内圧センサなどノック発生時期を検出するセンサを設けてもかまわない（請求項２に記載の発明）。

実施形態ではデポジット量に基づいてノックが生じないように点火時期を遅角側に制御する場合で説明したが、デポジット量に基づいてノックが生じないように吸気弁開時期を遅角側に制御するようにしてもかまわない（請求項１１に記載の発明）。

実施形態では、自着火時期θｋｎｋ（ノック発生時期）に基づいて燃焼室５の自着火時期における容積Ｖｋｎｋ（ノック発生時期の燃焼室容積）を算出し、この燃焼室５の自着火時期における容積Ｖｋｎｋと自着火時の燃焼質量割合ＢＲｋｎｋ（ノック発生時期の燃焼質量割合）に基づいてノック強度推定値ＫＩＣを算出し、このノック強度推定値ＫＩＣと実ノック強度ＫＩＤに基づいて燃焼室５内のデポジット量ＶＤＥＰを算出する場合で説明したが（請求項４に記載の発明）、自着火時期θｋｎｋ（ノック発生時期）と自着火時の燃焼質量割合ＢＲｋｎｋ（ノック発生時期の燃焼質量割合）に基づいてノック強度推定値ＫＩＣを算出し、このノック強度推定値ＫＩＣと実ノック強度ＫＩＤに基づいて燃焼室５内のデポジット量ＶＤＥＰを算出するようにしてもかまわない（請求項３に記載の発明）。さらに、自着火時期θｋｎｋ、自着火時の燃焼質量割合ＢＲｋｎｋ及び実ノック強度ＫＩＤに基づいて燃焼室５内のデポジット量ＶＤＥＰを算出するようにしてもかまわない（請求項２に記載の発明）。

また、燃焼質量割合に代えて熱発生率（燃焼質量割合の積分値である）を用いてもかまわない（請求項１、６、８、１２に記載の発明）。

実施形態では、上記（１２）式、（１９）式で説明したように、燃焼室５の燃焼開始時期における容積Ｖ０からデポジット量ＶＤＥＰを減算する場合で説明したが、上記（２）式の容積ＶＩＶＣから、あるいは上記（４３）式、（２２）式において燃焼室５の圧縮上死点時における容積ＶＴＤＣからこのデポジット量ＶＤＥＰを減算するようにしてもかまわない。

請求項２、９に記載の発明において、燃焼質量割合算出手段の機能は図３２のステップ１６４により、ノック発生時期検出・推定手段の機能は図３３のステップ１７２〜１８２により、ノック発生時燃焼質量割合算出手段の機能は図３３のステップ１８３により、デポジット量算出手段の機能は図３７のステップ２２６により、点火時期制御手段の機能は図３９のステップ２１２、２１３によりそれぞれ果たされている。

一実施形態のエンジンの制御システム図。 エンジンコントローラで実行される点火時期制御のブロック図。 燃焼室の圧力変化図。 燃焼質量割合の変化を説明する特性図。 物理量の算出を説明するためのフローチャート。 エンジンのクランクシャフトとコネクティングロッドの位置関係を説明するダイアグラム。 水温補正係数の特性図。 当量比補正係数の特性図。 基準クランク角の特性図。 初期燃焼期間の算出を説明するためのフローチャート。 温度上昇率の特性図。 主燃焼期間の算出を説明するためのフローチャート。 基本点火時期の算出を説明するためのフローチャート。 内部不活性ガス率の算出を説明するためのフローチャート。 内部不活性ガス量の算出を説明するためのフローチャート。 ＥＶＣ時不活性ガス量の算出を説明するためのフローチャート。 オーバーラップ中吹き返し不活性ガス量の算出を説明するためのフローチャート。 過給判定フラグ、チョーク判定フラグの設定を説明するためのフローチャート。 過給無しかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 過給無しかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 過給有りかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 過給有りかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 排気弁閉時期における燃焼室容積の特性図。 不活性ガスのガス定数の特性図。 オーバーラップ中の積算有効面積の特性図。 オーバーラップ中の積算有効面積の説明図。 不活性ガスの比熱比の特性図。 混合気の比熱比の特性図。 ノック発生時の燃焼室内の圧力履歴を示す特性図。 １／τの特性図。 直線で近似した場合の燃焼質量割合の変化を示す特性図。 ノック強度推定値の算出に必要となる物理量の算出を説明するためのフローチャート。 一燃焼サイクル毎の自着火時期の算出を説明するためのフローチャート。 一燃焼サイクル毎のノック強度推定値の算出を説明するためのフローチャート。 ガスエンタルピの算出を説明するためのフローチャート。 回転速度補正係数の特性図。 燃焼室内のデポジット量の算出を説明するためのフローチャート。 燃焼室内のデポジット量の特性図。 点火時期指令値の算出を説明するためのフローチャート。 点火時期指令値算出部のブロック図。

符号の説明

１ エンジン
５ 燃焼室
１１ 点火装置（火花点火手段）
１５ 吸気弁
２１ 燃料インジェクタ
２７ 吸気ＶＴＣ機構（可変バルブタイミング機構）
３１ エンジンコントローラ
４７ ノックセンサ（実ノック強度検出手段）

Claims (11)

1. 熱発生率を算出する熱発生率算出手段と、
   ノック発生時期を検出または推定するノック発生時期検出・推定手段と、
   実ノック強度を検出する実ノック強度検出手段と、
   これら熱発生率、ノック発生時期及び実ノック強度に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段と
   を備える
   を備えることを特徴とするエンジンのデポジット量検出装置。
2. 燃焼質量割合を算出する燃焼質量割合算出手段と、
   ノック発生時期を検出または推定するノック発生時期検出・推定手段と、
   この燃焼質量割合と前記ノック発生時期に基づいてノック発生時期の燃焼質量割合を算出するノック発生時燃焼質量割合算出手段と、
   実際のノック強度を検出する実ノック強度検出手段と、
   これらノック発生時期、ノック発生時の燃焼質量割合及び実ノック強度に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段と
   を備えることを特徴とするエンジンのデポジット量検出装置。
3. 前記ノック発生時期と前記ノック発生時期の燃焼質量割合に基づいてノック強度推定値を算出し、このノック強度推定値と前記実ノック強度に基づいて燃焼室内のデポジット量を算出することを特徴とする請求項２に記載のエンジンのデポジット量検出装置。
4. 前記ノック発生時期に基づいてノック発生時期の燃焼室容積を算出し、このノック発生時期の燃焼室容積と前記ノック発生時期の燃焼質量割合に基づいて前記ノック強度推定値を算出することを特徴とする請求項３に記載のエンジンのデポジット量検出装置。
5. 前記燃焼室内のデポジット量は、前記ノック強度推定値と前記実ノック強度の差が大きいほど大きいことを特徴とする請求項３に記載のエンジンのデポジット量検出装置。
6. 前記ノック発生時期を前記熱発生率に基づいて推定することを特徴とする請求項１に記載のエンジンのデポジット量検出装置。
7. 前記ノック発生時期を前記燃焼質量割合に基づいて推定することを特徴とする請求項２に記載のエンジンのデポジット量検出装置。
8. 熱発生率を算出する熱発生率算出手段と、
   ノック発生時期を検出または推定するノック発生時期検出・推定手段と、
   実ノック強度を検出する実ノック強度検出手段と、
   これら熱発生率、ノック発生時期及び実ノック強度に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段と、
   このデポジット量に基づいてノックが生じないように点火時期を遅角側に制御する点火時期制御手段と
   を備えることを特徴とするエンジンのノック制御装置。
9. 燃焼質量割合を算出する手段と、
   ノック発生時期を検出または推定するノック発生時期検出・推定手段と、
   この燃焼質量割合と前記ノック発生時期に基づいてノック発生時期の燃焼質量割合を算出するノック発生時燃焼質量割合算出手段と、
   実際のノック強度を検出する実ノック強度検出手段と、
   これらノック発生時期、ノック発生時の燃焼質量割合及び実ノック強度に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段と、
   このデポジット量に基づいてノックが生じないように点火時期を遅角側に制御する点火時期制御手段と
   を備えることを特徴とするエンジンのノック制御装置。
10. 吸気弁開閉時期を可変に調整し得る可変バルブタイミング機構と、
    熱発生率を算出する熱発生率算出手段と、
    ノック発生時期を検出または推定するノック発生時期検出・推定手段と、
    実ノック強度を検出する実ノック強度検出手段と、
    これら熱発生率、ノック発生時期及び実ノック強度に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段と、
    このデポジット量に基づいてノックが生じないように吸気弁開時期を遅角側に制御する吸気弁開時期制御手段と
    を備えることを特徴とするエンジンのノック制御装置。
11. 吸気弁開閉時期を可変に調整し得る可変バルブタイミング機構と、
    燃焼質量割合を算出する手段と、
    ノック発生時期を検出または推定するノック発生時期検出・推定手段と、
    この燃焼質量割合と前記ノック発生時期に基づいてノック発生時期の燃焼質量割合を算出するノック発生時燃焼質量割合算出手段と、
    実際のノック強度を検出する実ノック強度検出手段と、
    これらノック発生時期、ノック発生時の燃焼質量割合及び実ノック強度）に基づいて燃焼室内のデポジット量を算出するデポジット量算出手段と、
    このデポジット量に基づいてノックが生じないように吸気弁開時期を遅角側に制御する吸気弁開時期制御手段と
    を備えることを特徴とするエンジンのノック制御装置。

Images