(57)【要約】【課題】未燃ガス密度を用いることなく燃焼期間を算出すると共に、特に初期燃焼期間における燃焼ガス体積相当容積を制御上支障なく算出することを可能とする。【解決手段】層流燃焼速度、燃焼ガス体積相当容積、燃焼ガス質量及び反応確率に基づいて燃焼開始から所定クランク角までの燃焼期間を算出する手段（３１）と、この燃焼期間に基づいてＭＢＴの得られる基本点火時期を算出する手段（３１）と、この基本点火時期で火花点火を行う手段（１１）と、前記燃焼ガスの燃焼開始時期推定値を算出する手段（３１）とを備え、前記燃焼ガス体積相当容積算出手段が燃焼ガス体積相当容積を算出する際にこの燃焼開始時期推定値を用いる。【選択図】        図１
2004332653.tif 000002
【特許請求の範囲】
【請求項１】
燃焼ガスの層流状態での燃焼速度である層流燃焼速度を算出する層流燃焼速度算出手段と、
燃焼室内の燃焼ガス体積に相当する容積を算出する燃焼ガス体積相当容積算出手段と、
所定クランク角までに前記燃焼室内で燃焼する燃焼ガス質量を算出する燃焼ガス質量算出手段と、
所定運転条件での燃焼ガスの燃焼のしやすさを示す反応確率を算出する反応確率算出手段と、
これら層流燃焼速度、燃焼ガス体積相当容積、燃焼ガス質量及び反応確率に基づいて燃焼開始から所定クランク角までの燃焼期間を算出する燃焼期間算出手段と、
この燃焼期間に基づいてＭＢＴの得られる基本点火時期を算出する基本点火時期算出手段と、
この基本点火時期で火花点火を行う火花点火手段と、
前記燃焼ガスの燃焼開始時期推定値を算出する燃焼開始時期推定値算出手段と
を備え、
前記燃焼ガス体積相当容積算出手段が燃焼ガス体積相当容積を算出する際にこの燃焼開始時期推定値を用いることを特徴とするエンジンの制御装置。
【請求項２】
前記燃焼開始時期推定値を、前記基本点火時期の算出値の１サイクル前の値に基づいて算出することを特徴とする請求項１に記載のエンジンの制御装置。
【請求項３】
前記燃焼開始時期推定値を、前記燃焼室内の不活性ガス状態を表す指標とエンジン回転速度に基づいて算出することを特徴とする請求項１に記載のエンジンの制御装置。
【請求項４】
前記燃焼室内の不活性ガス状態を表す指標は総ガス質量と燃料質量の比であることを特徴とする請求項３に記載のエンジンの制御装置。
【請求項５】
前記燃焼室内の不活性ガス状態を表す指標は不活性ガス率であることを特徴とする請求項３に記載のエンジンの制御装置。
【請求項６】
前記所定クランク角は、ＭＢＴで点火を行った際に混合気の燃焼圧力が最大となる基準クランク角であることを特徴とする請求項１に記載のエンジンの制御装置。
【請求項７】
前記燃焼開始から基準クランク角までの燃焼期間を、燃焼開始から火炎核が形成されるまでの初期燃焼期間と、火炎核の形成後から基準クランク角までの主燃焼期間とに分割する場合に、前記燃焼開始時期推定値算出手段が、前記主燃焼期間に比例させて初期燃焼期間推定値を算出する手段と、この初期燃焼期間推定値と主燃焼期間とに基づいて前記燃焼開始時期推定値を算出する手段とからなることを特徴とする請求項６に記載のエンジンの制御装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、エンジン（内燃機関）の点火時期制御装置、特に最大の軸トルクを発生するのに必要な最小点火進角値（いわゆるＭＢＴ）となるように点火時期を制御するものに関する。
【０００２】
【従来の技術】
シリンダ内総ガス質量ＭＡＳＳＣを未燃ガス密度基本値ＤＥＮＳ及び層流火炎速度ＳＬＶで割った値に所定の着火遅れ時間Ｂ１を加算し、この加算値をクランク角に単位換算した値をＭＢＴの得られる基本点火時期として演算するものがある（特許文献１参照）。
【０００３】
【特許文献１】
特開平１０−３０５３５号公報
【０００４】
【発明が解決しようとする課題】
ところで、燃焼解析により得た知見を図４を参照しながら説明すると、図４は吸気下死点から膨張行程前半にかけての燃焼質量割合ＢＲの変化を示している。ここで、燃焼質量割合ＢＲは燃焼室に供給された燃料に対する燃焼ガス質量の比率を表し、これが燃焼解析により検出することが可能となった物理量である。この燃焼ガス質量ＢＲは燃焼開始時に０％であり、完全燃焼によって１００％に達する。この場合に、基準クランク角θＰＭＡＸにおける燃焼質量割合ＢＲは一定で約６０％であることを新たに見出した。
【０００５】
なお、図３は吸気下死点から膨張行程前半にかけての燃焼室圧力の変化を示し、ＭＢＴで混合気に点火した場合に混合気の燃焼圧力が最大値Ｐｍａｘとなるクランク角が基準クランク角θＰＭＡＸである。この基準クランク角θＰＭＡＸは燃焼方式によらずほぼ一定であり、一般に圧縮上死点後１２〜１５度、最大で圧縮上死点後１０〜２０度の範囲にあることが知られている。
【０００６】
ここで、燃焼を開始してから基準クランク角θＰＭＡＸに達するまでの期間を燃焼期間として定義すれば、この燃焼期間は燃焼室の総ガス質量を層流燃焼速度で除算することによって求めることが基本的に可能であるため、上記従来装置では、総ガス質量ＭＡＳＳＣを未燃ガス密度基本値ＤＥＮＳと乱流火炎速度基本値ＦＬＭＴとで除算して燃焼期間を求めている。従来装置における未燃ガス密度基本値ＤＥＮＳは未燃ガス質量を未燃ガス体積で割って得られる値であるため、理論的にはこれら未燃ガス質量及び未燃ガス体積を検出すれば未燃ガス密度基本値ＤＥＮＳを正確に求めることができるのであるが、実際には燃焼室内にける未燃ガス体積を推定することは困難であるので、従来装置においては充填効率ＩＴＡＣに基づいて未燃ガス密度基本値ＤＥＮＳを求めている。
【０００７】
しかしながら、従来装置のように質量に相当する充填効率ＩＴＡＣのみの関数で未燃ガス密度基本値ＤＥＮＳを算出するのでは、運転条件により変化する未燃ガス体積分だけ算出精度が悪くなり、算出精度が落ちる分だけ燃焼期間の算出精度が低下し、基本点火時期の演算が不正確となり、ＭＢＴが得られない事態が生じ得る。
【０００８】
そこで本発明は、燃焼ガス体積はそのときの燃焼室容積で近似し、代わりに燃焼ガス質量を導入し、これら燃焼室容積、燃焼ガス質量に基づいて燃焼期間を算出すると共に、燃焼室内の燃焼ガス体積に相当する容積を算出する際に、燃焼開始時期推定値を用いることにより、最近の燃焼解析結果を盛り込んだ新たな点火時期制御方法を提供することを目的とする。
【０００９】
【課題を解決するための手段】
本発明は、燃焼ガスの層流状態での燃焼速度である層流燃焼速度、燃焼室内の燃焼ガス体積に相当する容積、所定クランク角までに前記燃焼室内で燃焼する燃焼ガス質量、所定運転条件での燃焼ガスの燃焼のしやすさを示す反応確率をそれぞれ算出し、これら層流燃焼速度、燃焼ガス体積相当容積、燃焼ガス質量及び反応確率に基づいて燃焼開始から所定クランク角までの燃焼期間を算出し、この燃焼期間に基づいてＭＢＴの得られる基本点火時期を算出するように構成し、このようにして得た基本点火時期で火花点火を行うと共に、前記燃焼ガスの燃焼開始時期推定値を算出し、前記燃焼ガス体積相当容積を算出する際にこの燃焼開始時期推定値を用いるようにしたものである。
【００１０】
【発明の効果】
本発明によれば、層流燃焼速度、燃焼ガス体積に近似させた燃焼ガス体積相当容積、新たに導入した燃焼ガス質量及び反応確率に基づいて燃焼期間を算出するようにしたので、従来装置のように未燃ガス密度を用いることなく燃焼期間を算出することが可能となり、これにより正確かつ容易にＭＢＴの得られる基本点火時期を算出することができる。
【００１１】
また、特に初期燃焼期間における燃焼ガス体積相当容積を算出するには燃焼開始時期が必要となる。この場合に、点火時期と燃焼開始時期との間には密接な関係があるものの（図４参照）、これから点火時期を算出しようとしている段階で、点火時期の後にくる燃焼開始時期を知ることはできないのであるが、本発明によれば、燃焼ガスの燃焼開始時期推定値を算出し、燃焼ガス体積相当容積を算出する際にこの燃焼開始時期推定値を用いるようにしているので、特に初期燃焼期間における燃焼ガス体積相当容積を、点火時期制御上支障なく算出することができる。
【００１２】
【発明の実施の形態】
以下、図面に基づき本発明の実施形態について説明する。
【００１３】
図１は、本発明のシステムを説明するための概略図である。
【００１４】
空気は吸気コレクタ２に蓄えられた後、吸気マニホールド３を介して各気筒の燃焼室５に導入される。燃料は各気筒の吸気ポート４に配置された燃料インジェクタ２１より噴射供給される。空気中に噴射された燃料は気化しつつ空気と混合してガス（混合気）を作り、燃焼室５に流入する。この混合気は吸気弁１５が閉じることで燃焼室５内に閉じこめられ、ピストン６の上昇によって圧縮される。
【００１５】
この圧縮混合気に対して高圧火花により点火を行うため、パワートランジスタ内蔵の点火コイルを各気筒に配した電子配電システムの点火装置１１を備える。すなわち、点火装置１１は、バッテリからの電気エネルギーを蓄える点火コイル１３と、点火コイル１３の一次側への通電、遮断を行うパワートランジスタと、燃焼室５の天井に設けられ点火コイル１３の一次電流の遮断によって点火コイル１３の二次側に発生する高電圧を受けて、火花放電を行う点火プラグ１４とからなっている。
【００１６】
圧縮上死点より少し手前で点火プラグ１４により火花が飛ばされ圧縮混合気に着火されると、火炎が広がりやがて爆発的に燃焼し、この燃焼によるガス圧がピストン６を押し下げる仕事を行う。この仕事はクランクシャフト７の回転力として取り出される。燃焼後のガス（排気）は排気弁１６が開いたとき排気通路８へと排出される。
【００１７】
排気通路８には三元触媒９を備える。三元触媒９は排気の空燃比が理論空燃比を中心とした狭い範囲（ウインドウ）にあるとき、排気に含まれるＨＣ、ＣＯ、ＮＯｘといった有害三成分を同時に効率よく除去できる。空燃比は吸入空気量と燃料量の比であるので、エンジンの１サイクル（４サイクルエンジンではクランク角で７２０°区間）当たりに燃焼室５に導入される吸入空気量と、燃料インジェクタ２１からの燃料噴射量との比が理論空燃比となるように、エンジンコントローラ３１ではエアフローメータ３２からの吸入空気流量の信号とクランク角センサ（３３、３４）からの信号に基づいて燃料インジェクタ２１からの燃料噴射量を定めると共に、三元触媒９の上流に設けたＯ２センサ３５からの信号に基づいて空燃比をフィードバック制御している。
【００１８】
吸気コレクタ２の上流には絞り弁２３がスロットルモータ２４により駆動される、いわゆる電子制御スロットル２２を備える。運転者が要求するトルクはアクセルペダル４１の踏み込み量（アクセル開度）に現れるので、エンジンコントローラ３１ではアクセルセンサ４２からの信号に基づいて目標トルクを定め、この目標トルクを実現するための目標空気量を定め、この目標空気量が得られるようにスロットルモータ２４を介して絞り弁２３の開度を制御する。
【００１９】
吸気弁用カムシャフト２５、排気弁用カムシャフト２６及びクランクシャフト７の各前部にはそれぞれカムスプロケット、クランクスプロケットが取り付けられ、これらスプロケットにタイミングチェーン（図示しない）を掛け回すことで、カムシャフト２５、２６がエンジンのクランクシャフト７により駆動されるのであるが、このカムスプロケットと吸気弁用カムシャフト２５との間に介在して、作動角一定のまま吸気弁用カムの位相を連続的に制御し得る吸気バルブタイミングコントロール機構（以下、「吸気ＶＴＣ機構」という。）２７と、カムスプロケットと排気弁用カムシャフト２６との間に介在して、作動角一定のまま排気弁用カムの位相を連続的に制御し得る排気バルブタイミングコントロール機構（以下、「排気ＶＴＣ機構」という。）２８とを備える。吸気弁１５の開閉時期や排気弁１６の開閉時期を変えると燃焼室５に残留する不活性ガスの量が変化する。燃焼室５内の不活性ガスの量が増えるほどポンピングロスが減って燃費がよくなるので、運転条件によりどのくらいの不活性ガスが燃焼室５内に残留したらよいかを目標吸気弁閉時期や目標排気弁閉時期にして予め定めており、エンジンコントローラ３１ではそのときの運転条件（エンジンの負荷と回転速度）より目標吸気弁閉時期と目標排気弁閉時期を定め、それら目標値が得られるように吸気ＶＴＣ機構２７、排気ＶＴＣ機構２８の各アクチュエータを介して吸気弁閉時期と排気弁閉時期を制御する。
【００２０】
吸気温度センサ４３からの吸気温度の信号、吸気圧力センサ４４からの吸気圧力の信号、排気温度センサ４５からの排気温度の信号、排気圧力センサ４６からの排気圧力の信号が、水温センサ３７からの冷却水温の信号と共に入力されるエンジンコントローラ３１では、パワートランジスタ１３を介して点火プラグ１４の一次側電流の遮断時期である点火時期を制御する。
【００２１】
図２はエンジンコントローラ３１内で行われる点火時期制御のブロック図で、大きくは点火時期演算部５１と点火時期制御部６１とからなる。点火時期演算部５１はさらに初期燃焼期間算出部５２、主燃焼期間算出部５３、燃焼期間算出部５４、基本点火時期算出部５５及び燃焼開始時期推定値算出部５６からなる。
【００２２】
初期燃焼期間算出部５２では、混合気が着火してから火炎核が形成されるまでの期間を初期燃焼期間ＢＵＲＮ１として算出する。主燃焼期間算出部５３では、火炎核が形成されてから燃焼圧力が最大値Ｐｍａｘに達するまでの期間を主燃焼期間ＢＵＲＮ２として算出する。燃焼期間算出部５４では、これら初期燃焼期間ＢＵＲＮ１と主燃焼期間ＢＵＲＮ２との合計を、燃焼開始より最大燃焼圧力Ｐｍａｘに至るまでの燃焼期間ＢＵＲＮとして算出する。基本点火時期算出部５５では、この燃焼期間ＢＵＲＮに基づいてＭＢＴの得られる点火時期（この点火時期を「基本点火時期」という。）ＭＢＴＣＡＬを算出する。
【００２３】
点火時期制御部６１ではこのようにして算出された基本点火時期を点火時期指令値とし、この指令値で点火プラグ１４が燃焼室５内の混合気に対して着火するように、イグニッションコイル１３への通電角と非通電角を制御する。
【００２４】
上記のように燃焼期間ＢＵＲＮを初期燃焼期間ＢＵＲＮ１と主燃焼期間ＢＵＲＮ２に分けて算出し、燃焼期間ＢＵＲＮに応じて基本点火時期ＭＢＴＣＡＬを求めるようにしたのは、燃焼解析より得られた結果に基づくものである。以下、燃焼解析に基づくこの点火時期制御をさらに説明する。
【００２５】
図３に示すようにＭＢＴ（最大トルクの得られる最小進角値）で混合気に点火した場合に混合気の燃焼圧力が最大値Ｐｍａｘとなるクランク角を基準クランク角θＰＭＡＸ［ｄｅｇＡＴＤＣ］とする。基準クランク角θＰＭＡＸは燃焼方式によらずほぼ一定であり、一般に圧縮上死点後１２〜１５度、最大で圧縮上死点後１０〜２０度の範囲にある。
【００２６】
図４に火花点火エンジンにおける燃焼室内の燃焼解析により得られた燃焼質量割合ＢＲ（燃焼ガス質量割合）の変化を示す。燃焼室に供給された燃料に対する燃焼質量の比率を表す燃焼質量割合ＢＲは、燃焼開始時に０％であり、完全燃焼によって１００％に達する。基準クランク角θＰＭＡＸにおける燃焼質量割合は一定で約６０％であることが実験により確かめられている。
【００２７】
燃焼質量割合ＢＲが０％から基準クランク角θＰＭＡＸ相当の約６０％に達するまでの変化代に相当する燃焼期間は、燃焼開始直後で燃焼質量割合にも燃焼圧力にもほとんど変化のない期間である初期燃焼期間と、燃焼質量割合と燃焼圧力が急激に増加する主燃焼期間とに分けられる。初期燃焼期間は、燃焼開始から火炎核が形成されるまでの段階であり、火炎核が形成されるのは燃焼質量割合が０％から２％〜１０％まで変化したときである。この初期燃焼期間中は、燃焼圧力や燃焼温度の上昇速度が小さく、燃焼質量割合の変化に対して初期燃焼期間は長い。初期燃焼期間の長さは燃焼室内の温度や圧力の変化の影響を受けやすい。
【００２８】
一方、主燃焼期間においては、火炎核から外側へと火炎が伝播するのであり、その火炎速度（つまり燃焼速度）が急上昇する。そのため、主燃焼期間の燃焼質量割合の変化は初期燃焼期間の燃焼質量割合の変化に比べて大きい。
【００２９】
エンジンコントローラ３１では、燃焼質量割合が２％に達する（変化する）までを初期燃焼期間ＢＵＲＮ１［ｄｅｇ］とし、初期燃焼期間ＢＵＲＮ１の終了後、基準クランク角θＰＭＡＸに至るまでの区間（燃焼室量割合でいえば２％より約６０％に達するまでの間）を主燃焼期間ＢＵＲＮ２［ｄｅｇ］として区別する。そして、初期燃焼期間ＢＵＲＮ１に主燃焼期間ＢＵＲＮ２を加えた合計である燃焼期間ＢＵＲＮ［ｄｅｇ］を算出し、この燃焼期間ＢＵＲＮから基準クランク角θＰＭＡＸ［ｄｅｇＡＴＤＣ］を差し引き、さらに後述する点火無駄時間相当クランク角ＩＧＮＤＥＡＤ［ｄｅｇ］を加えたクランク角位置を、ＭＢＴの得られる点火時期である基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］として設定する。
【００３０】
火炎核の形成される初期燃焼期間での燃焼室５内の圧力、温度は、点火時の圧力、温度とほぼ等価になるが、これから点火時期を算出しようとしているのに、最初から正確な点火時期を設定することはできない。そこで、図２に示したように燃焼開始時期推定値算出部５６で基本点火時期の前回値を燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ［ｄｅｇＢＴＤＣ］として算出し、この値を初期燃焼期間算出部５２に対して与えるようにし、初期燃焼期間算出部５２において初期燃焼期間の算出をサイクリックに繰り返すことで、精度の高い結果を時間遅れなしに出すようにしている。
【００３１】
次に、エンジンコントローラ３１で実行される上記の基本点火時期ＭＢＴＣＡＬの算出を以下のフローチャートを参照しながら詳述する。
【００３２】
図５は点火時期の算出に必要な各種の物理量を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。
【００３３】
まずステップ１１では、吸気弁閉時期ＩＶＣ［ｄｅｇＢＴＤＣ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、内部不活性ガス率ＭＲＥＳＦＲ［％］、温度センサ３７により検出される冷却水温ＴＷＫ［Ｋ］、目標当量比ＴＦＢＹＡ、クランク角センサにより検出されるエンジン回転速度ＮＲＰＭ［ｒｐｍ］、点火無駄時間ＤＥＡＤＴＩＭＥ［μｓｅｃ］を読み込む。
【００３４】
ここで、クランク角センサはクランクシャフト７のポジションを検出するポジションセンサ３３と、吸気用カムシャフト２５ポジションを検出するフェーズセンサ３４とからなり、これら２つのセンサ３３、３４からの信号に基づいてエンジン回転速度ＮＲＰＭ［ｒｐｍ］が算出されている。
【００３５】
吸気弁閉時期ＩＶＣは吸気ＶＴＣ機構２７に与える指令値から既知である。あるいはフェーズセンサ３４により実際の吸気弁閉時期を検出してもかまわない。
【００３６】
内部不活性ガス率ＭＲＥＳＦＲは燃焼室内に残留する不活性ガス量を燃焼室内の総ガス量で除した値で、その算出については後述する。点火無駄時間ＤＥＡＤＴＩＭＥは一定値である。
【００３７】
目標当量比ＴＦＢＹＡは図示しない燃料噴射量の算出フローにおいて算出されている。目標当量比ＴＦＢＹＡは無名数であり、理論空燃比を１４．７とすると、次式により表される値である。
【００３８】
ＴＦＢＹＡ＝１４．７／目標空燃比…（１）
例えば（１）式より目標空燃比が理論空燃比のときＴＦＢＹＡ＝１．０となり、目標空燃比が例えば２２．０といったリーン側の値であるとき、ＴＦＢＹＡは１．０未満の正の値である。
【００３９】
ステップ１２では、燃焼室５の吸気弁閉時期ＩＶＣにおける容積（つまり圧縮開始時期での容積）ＶＩＶＣ［ｍ３］を算出する。燃焼室５の吸気弁閉時期における容積ＶＩＶＣは、ピストン６のストローク位置によって決まる。ピストン６のストローク位置はエンジンのクランク角位置によって決まる。
【００４０】
図６を参照して、エンジンのクランクシャフト７１の回転中心７２がシリンダの中心軸７３からオフセットしている場合を考える。コネクティングロッド７４、コネクティングロッド７４とクランクシャフト７１との結節点７５、コネクティングロッド７４とピストンをつなぐピストンピン７６が図に示す関係にあるとする。このときの、燃焼室５の吸気弁閉時期における容積ＶＩＶＣは次式（２）〜（６）で表すことができる。
【００４１】
2004332653.tif 000003
ただし、Ｖｃ：隙間容積［ｍ３］、
ε                  ：圧縮比、
Ｄ                  ：シリンダボア径［ｍ］、
ＳＴ              ：ピストンの全ストローク［ｍ］、
Ｈｉｖｃ            ：吸気弁閉時期におけるピストンピン７６のＴＤＣからの距離［ｍ］、
Ｈｘ：      ：ピストンピン７６のＴＤＣからの距離の最大値と最小値の差［ｍ］、
ＣＮＤ          ：コネクティングロッド７４の長さ［ｍ］、
ＣＲｏｆｆ        ：結節点７５のシリンダ中心軸７３からのオフセット距離［ｍ］、
ＰＩＳｏｆｆ    ：クランクシャフト回転中心７２のシリンダ中心軸７３からのオフセット距離［ｍ］、
θｉｖｃ            ：吸気弁閉時期のクランク角［ｄｅｇＡＴＤＣ］、
θｏｆｆ            ：ピストンピン７６とクランクシャフト回転中心７２とを結ぶ線がＴＤＣにおいて垂直線となす角度［ｄｅｇ］、
Ｘ                  ：結節点７５とピストンピン７６との水平距離［ｍ］、
吸気弁閉時期のクランク角θｉｖｃは前述のように、エンジンコントローラ３１から吸気ＶＴＣ機構２７への指令信号によって決まるので、既知である。式（２）〜（６）にこのときのクランク角θｉｖｃ（＝ＩＶＣ）を代入すれば、燃焼室５の吸気弁閉時期における容積ＶＩＶＣを算出することができる。したがって、実用上は燃焼室５の吸気弁閉時期における容積ＶＩＶＣは吸気弁閉時期ＩＶＣをパラメータとするテーブルで設定したものを用いる。吸気ＶＴＣ機構２７を備えないときには定数で与えることができる。
【００４２】
ステップ１３では、燃焼室５の吸気弁閉時期ＩＶＣにおける温度（つまり圧縮開始時期温度）ＴＩＮＩ［Ｋ］を算出する。燃焼室５に流入するガスの温度は、燃焼室５に流入する新気と燃焼室５に残留する不活性ガスとが混じったガスの温度であり、燃焼室５に流入する新気の温度は吸気コレクタ２内の新気温度ＴＣＯＬに等しく、また燃焼室５内に残留する不活性ガスの温度は排気ポート部近傍の排気温度ＴＥＸＨで近似できるので、燃焼室５の吸気弁閉時期ＩＶＣにおける温度ＴＩＮＩは吸気弁閉時期ＩＶＣになったタイミングでの、吸気コレクタ２内の新気温度ＴＣＯＬ、排気温度ＴＥＸＨ、燃焼室５内に残留する不活性ガスの割合である内部不活性ガス率ＭＲＥＳＦＲから次式により求めることができる。
【００４３】
ＴＩＮＩ＝ＴＥＸＨ×ＭＲＥＳＦＲ＋ＴＣＯＬ×（１−ＭＲＥＳＦＲ）…（７）
ステップ１４では燃焼室５の吸気弁閉時期ＩＶＣにおける圧力（つまり圧縮開始時期圧力）ＰＩＮＩ［Ｐａ］を算出する。すなわち、吸気弁閉時期ＩＶＣになったタイミングでのコレクタ内圧力ＰＣＯＬを吸気弁閉時期ＩＶＣにおける圧力ＰＩＮＩとして取り込む。
【００４４】
ステップ１５では、燃焼室５内の混合気の燃えやすさを表す反応確率ＲＰＲＯＢＡ［％］を算出する。反応確率ＲＰＲＯＢＡは無次元の値であり、残留不活性ガス率ＭＲＥＳＦＲ、冷却水温ＴＷＫ［Ｋ］、目標当量比ＴＦＢＹＡの３つのパラメータに依存するので、次式により表すことができる。
【００４５】
ＲＰＲＯＢＡ＝ｆ３（ＭＲＥＳＦＲ、ＴＷＫ、ＴＦＢＹＡ）…（８）
具体的に説明すると、ＭＲＥＳＦＲ、ＴＷＫ、ＴＦＢＹＡの３つのパラメータの組み合わせによって得られる反応確率の最大値を１００％とし、これらのパラメータと反応確率ＲＰＲＯＢＡの関係を実験的に求め、求めた反応確率ＲＰＲＯＢＡをパラメータに応じたテーブルとしてエンジンコントローラ３１のメモリに予め格納しておく。ステップ１４ではパラメータに応じてこのテーブルを検索することにより反応確率ＲＰＲＯＢＡを求める。
【００４６】
具体的には、冷却水温ＴＷＫに応じて図７に示すような特性を有する水温補正係数のテーブルと、同様に設定された内部不活性ガス率補正係数のテーブル（図示しない）と、目標当量比ＴＦＢＹＡに応じて図８に示すような特性を有する当量比補正係数のテーブルを予めメモリに格納しておく。各補正係数の最大値はそれぞれ１．０であり、３種類の補正係数の積に反応確率の最大値１００％を掛け合わせることで、反応確率ＲＰＲＯＢＡを算出する。
【００４７】
各テーブルを説明すると、図７に示す水温補正係数は冷却水温ＴＷＫが高いほど大きく、冷却水温ＴＷＫが８０℃以上では１．０になる。図８に示す当量比補正係数は目標当量比ＴＦＢＹＡが１．０のとき、つまり理論空燃比のときに最大値の１．０となり、目標当量比が１．０より大きくても小さくても当量比補正係数は減少する。内部不活性ガス率補正係数は図示しないが、内部不活性ガス率ＭＲＥＳＦＲがゼロの場合に１．０となる。
【００４８】
ステップ１６では、基準クランク角θＰＭＡＸ［ｄｅｇＡＴＤＣ］を算出する。前述のように基準クランク角θＰＭＡＸはあまり変動しないが、それでもエンジン回転速度ＮＲＰＭの上昇に応じて進角する傾向があるため、基準クランク角θＰＭＡＸはエンジン回転速度ＮＲＰＭの関数として次式で表すことができる。
【００４９】
θＰＭＡＸ＝ｆ４（ＮＲＰＭ）…（９）
具体的にはエンジン回転速度ＮＲＰＭから、エンジンコントローラ３１のメモリに予め格納された図９に示す特性のテーブルを検索することにより基準クランク角θＰＭＡＸを求める。算出を容易にするために、基準クランク角θＰＭＡＸを一定とみなすことも可能である。
【００５０】
最後にステップ１７では、点火無駄時間相当クランク角ＩＧＮＤＥＡＤ［ｄｅｇ］を算出する。点火無駄時間相当クランク角ＩＧＮＤＥＡＤは、エンジンコントローラ３１から点火コイル１３の一次電流を遮断する信号を出力したタイミングから点火プラグ１４が実際に点火するまでのクランク角区間で、次式により表すことができる。
【００５１】
ＩＧＮＤＥＡＤ＝ｆ５（ＤＥＡＤＴＩＭＥ、ＮＲＰＭ）…（１０）
ここでは、点火無駄時間ＤＥＡＤＴＩＭＥを２００μｓｅｃとする。（１０）式は、エンジン回転速度ＮＲＰＭから点火無駄時間ＤＥＡＤＴＩＭＥに相当するクランク角である点火無駄時間相当クランク角ＩＧＮＤＥＡＤを算出するためのものである。
【００５２】
図１０は初期燃焼期間ＢＵＲＮ１［ｄｅｇ］を算出するためのもの、また図１２は主燃焼期間ＢＵＲＮ２［ｄｅｇ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１０、図１２は図５に続けて実行する。
図１０、図１２はどちらを先に実行してもかまわない。
【００５３】
まず図１０から説明すると、ステップ１６１では、燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ［ｄｅｇＢＴＤＣ］、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ３］、図５のステップ１３で算出されている燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出されている燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、図５のステップ１５で算出されている反応確率ＲＰＲＯＢＡ［％］を読み込む。
【００５４】
ここで、燃焼開始時期推定値ＭＢＴＣＡＬ Ｅは、基本点火時期ＭＢＴＣＡＬの［ｄｅｇＢＴＤＣ］の１サイクル前の値に基づいて算出される値であり、その算出については図１３により後述する。
【００５５】
ステップ１６２では燃焼室５の燃焼開始時期における容積Ｖ０［ｍ３］を算出する。前述したように、ここでの燃焼開始時期は基本点火時期の１サイクル前の値に基づいて算出している。すなわち、基本点火時期の１サイクル前の値に基づいて算出される燃焼開始時期推定値ＭＢＴＣＡＬ Ｅから次式により燃焼室５の燃焼開始時期における容積Ｖ０を算出する。
【００５６】
Ｖ０＝ｆ６（ＭＢＴＣＡＬ Ｅ）…（１１）
具体的には燃焼開始時期推定値ＭＢＴＣＡＬ Ｅにおけるピストン６のストローク位置と、燃焼室５のボア径から、燃焼室５のＭＢＴＣＡＬ Ｅにおける容積Ｖ０を算出する。図５のステップ１２では、燃焼室５の吸気弁閉時期ＩＶＣにおける容積ＶＩＶＣを、吸気弁閉時期をパラメータとする吸気弁閉時期容積のテーブルを検索することにより求めたが、ここではＭＢＴＣＡＬ Ｅをパラメータとする前回燃焼開始時期容積のテーブルを検索することにより、燃焼室５の燃焼開始時期推定値ＭＢＴＣＡＬ Ｅにおける容積Ｖ０を求めればよい。
【００５７】
ステップ１６３では燃焼開始時期における有効圧縮比Ｅｃを算出する。有効圧縮比Ｅｃは無次元の値であり、次式に示すように燃焼室５の燃焼開始時期における容積Ｖ０を燃焼室５の吸気弁閉時期における容積ＶＩＶＣで除した値である。
【００５８】
Ｅｃ＝ｆ７（Ｖ０、ＶＩＶＣ）＝Ｖ０／ＶＩＶＣ…（１２）
ステップ１６４では吸気弁閉時期ＩＶＣから燃焼開始時期に至る間の燃焼室５内の温度上昇率ＴＣＯＭＰを次式に示すように有効圧縮比Ｅｃに基づいて算出する。
【００５９】
ＴＣＯＭＰ＝ｆ８（Ｅｃ）＝Ｅｃ＾（κ−１）…（１３）
ただし、κ：比熱比、
（１３）式は断熱圧縮されるガスの温度上昇率の式である。なお、（１３）式右辺の「＾」は累乗計算を表している。この記号は後述する式でも使用する。
【００６０】
κは断熱圧縮されるガスの定圧比熱を定容比熱で除した値で、断熱圧縮されるガスが空気であればκ＝１．４であり、簡単にはこの値を用いればよい。ただし、混合気に対してκの値を実験的に求めることで、一層の算出精度の向上が可能である。
【００６１】
図１１は（１３）式を図示したものである。従って、このような特性のテーブルを予めエンジンコントローラ３１のメモリに格納しておき、有効圧縮比Ｅｃに基づき当該テーブルを検索することにより温度上昇率ＴＣＯＭＰを求めることも可能である。
【００６２】
ステップ１６５では、燃焼室５の燃焼開始時期における温度Ｔ０［Ｋ］を、燃焼室５の吸気弁閉時期における温度ＴＩＮＩに温度上昇率ＴＣＯＭＰを乗じることで、つまり
Ｔ０＝ＴＩＮＩ×ＴＣＯＭＰ…（１４）
の式により算出する。
【００６３】
ステップ１６６、１６７はステップ１６４、１６５と同様である。すなわち、ステップ１６６では吸気弁閉時期ＩＶＣから燃焼開始時期に至る間の燃焼室５内の圧力上昇率ＰＣＯＭＰを次式に示すように有効圧縮比Ｅｃに基づいて算出する。
【００６４】
ＰＣＯＭＰ＝ｆ９（Ｅｃ）＝Ｅｃ＾κ…（４１）
ただし、κ：比熱比、
（４１）式も（１３）式と同じに断熱圧縮されるガスの圧力上昇率の式である。（４１）式右辺の「＾」も（１３）式と同じに累乗計算を表している。
【００６５】
κは上記（１３）式で用いている値と同じで、断熱圧縮されるガスが空気であればκ＝１．４であり、簡単にはこの値を用いればよい。ただし、混合気に対してその組成、温度からκの値を求めることで、一層の算出精度の向上が可能である。
【００６６】
図１１と同様の特性のテーブルを予めエンジンコントローラ３１のメモリに格納しておき、有効圧縮比Ｅｃに基づき当該テーブルを検索することにより圧力上昇率ＰＣＯＭＰを求めることも可能である。
【００６７】
ステップ１６７では、燃焼室５の燃焼開始時期における圧力Ｐ０［Ｐａ］を、燃焼室５の吸気弁閉時期における圧力ＰＩＮＩに圧力上昇率ＰＣＯＭＰを乗じることで、つまり
Ｐ０＝ＰＩＮＩ×ＰＣＯＭＰ…（４２）
の式により算出する。
【００６８】
ステップ１６８では、初期燃焼期間における層流燃焼速度ＳＬ１［ｍ／ｓｅｃ］を次式（公知）により算出する。
【００６９】
2004332653.tif 000004
ただし、Ｔｓｔｄ    ：基準温度［Ｋ］、
Ｐｓｔｄ    ：基準圧力［Ｐａ］、
ＳＬｓｔｄ：基準温度Ｔｓｔｄと基準圧力Ｐｓｔｄにおける基準層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ０      ：燃焼室５の燃焼開始時期における温度［Ｋ］、
Ｐ０      ：燃焼室５の燃焼開始時期における圧力［Ｐａ］、
層流燃焼速度（層流火炎速度）は気体の流れがない状態での火炎の伝播速度のことであり、燃焼室５内の圧縮速度、燃焼室５内の吸気流速に因らず、燃焼室５の温度及び圧力の関数となることが知られていることから、初期燃焼期間における層流燃焼速度を燃焼開始時温度Ｔ０と燃焼開始時圧力Ｐ０の関数として、また後述するように主燃焼期における層流燃焼速度を圧縮上死点時温度ＴＴＤＣと圧縮上死点圧力ＰＴＤＣの関数としている。これは、層流燃焼速度は一般的に、エンジン負荷、燃焼室５内の不活性ガス率、吸気弁閉時期、比熱比、吸気温度により変化するのであるが、これらは燃焼室５内の温度Ｔと圧力Ｐに影響する因子であるので、層流燃焼速度は最終的に燃焼室５内の温度Ｔと圧力Ｐにより規定できるとするものである。
【００７０】
上記の（１５）式において基準温度Ｔｓｔｄと基準圧力Ｐｓｔｄと基準層流燃焼速度ＳＬｓｔｄは実験により予め定められる値である。
【００７１】
燃焼室５の通常の圧力である２ｂａｒ以上の圧力下では、（１５）式の圧力項（Ｐ０／Ｐｓｔｄ）−０．１６は小さな値となる。従って、圧力項（Ｐ０／Ｐｓｔｄ）−０．１６を一定値として、基準層流燃焼速度ＳＬｓｔｄを基準温度Ｔｓｔｄのみで規定することも可能である。
【００７２】
従って、基準温度Ｔｓｔｄが５５０［Ｋ］で、基準層流燃焼速度ＳＬｓｔｄが１．０［ｍ／ｓｅｃ］で、圧力項が０．７である場合の燃焼開始時期における温度Ｔ０と層流燃焼速度ＳＬ１との関係は近似的に次式で定義することができる。
【００７３】
2004332653.tif 000005
ステップ１６９では、初期燃焼期間におけるガス流動の乱れ強さＳＴ１を算出する。このガス流動の乱れ強さＳＴ１は無次元の値であり、燃焼室５に流入する新気の流速と燃料インジェクタ２１の噴射燃料のペネトレーションとに依存する。
【００７４】
燃焼室５に流入する新気の流速は、吸気通路の形状と、吸気弁１５の作動状態と、吸気弁１５を設ける吸気ポート４の形状に依存する。噴射燃料のペネトレーションは燃料インジェクタ２１の噴射圧力と、燃料噴射期間と、燃焼噴射タイミングに依存する。
【００７５】
最終的に、初期燃焼期間におけるガス流動の乱れ強さＳＴ１は、エンジン回転速度ＮＲＰＭの関数として次式で表すことができる。
【００７６】
ＳＴ１＝ｆ１２（ＮＲＰＭ）＝Ｃ１×ＮＲＰＭ…（１７）
ただし、Ｃ１：定数、
乱れ強さＳＴ１を回転速度ＮＲＰＭをパラメータとするテーブルから求めることも可能である。
【００７７】
ステップ１７０では層流燃焼速度Ｓ１と乱れ強さＳＴ１から、初期燃焼期間におけるガスの燃焼速度ＦＬＡＭＥ１［ｍ／ｓｅｃ］を次式により算出する。
【００７８】
ＦＬＡＭＥ１＝ＳＬ１×ＳＴ１…（１８）
燃焼室５内にガス乱れがあるとガスの燃焼速度が変化する。（１８）式はこのガス乱れに伴う燃焼速度への寄与（影響）を考慮したものである。
【００７９】
ステップ１７１では、次式により初期燃焼期間ＢＵＲＮ１［ｄｅｇ］を算出する。
【００８０】
ＢＵＲＮ１＝｛（ＮＲＰＭ×６）×（ＢＲ１×Ｖ０）｝／（ＲＰＲＯＢＡ×ＡＦ１×ＦＬＡＭＥ１）…（１９）
ただし、ＡＦ１：火炎核の反応面積（固定値）［ｍ２］、
この（１９）式および後述する（２２）式は、燃焼ガス質量を燃焼速度で割ると燃焼期間が得られるとする次の基本式より導いたものであるが、（１９）、（２２）式右辺の分子、分母ががただちに燃焼ガス質量、燃焼速度を表すものではない。
【００８１】
燃焼期間［ｓｅｃ］＝シリンダ内総質量［ｇ］／（未燃ガス密度［ｇ／ｍ３］×火炎表面積［ｍ２］×火炎速度［ｍ／ｓｅｃ］）…（補１）
（補１）式右辺分母の未燃ガス密度は、未燃ガス質量［ｇ］を未燃ガス体積［ｍ３］で割った値であるので、従来装置のように質量に相当する充填効率ＩＴＡＣのみの関数では未燃ガス密度を正確に計算できているとはいえない。そこで、（補１）式に対して実験結果とを照らし合わせつつ所定の近似を導入して初めて得られたのが上記（１９）式及び後述する（２２）式に示す実験式である。
【００８２】
ここで、（１９）式右辺のＢＲ１は燃焼開始時期より初期燃焼期間ＢＵＲＮ１の終了時期までの燃焼質量割合の変化代であり、ここではＢＲ１＝２％に設定している。（１９）式右辺の（ＮＲＰＭ×６）は単位をｒｐｍからクランク角（ｄｅｇ）に変換するための処理である。火炎核の反応面積ＡＦ１は実験的に設定される。
【００８３】
また、初期燃焼期間中はほぼ燃焼室容積は変わらないとみなすことができる。従って、初期燃焼期間ＢＵＲＮ１を算出するに際して最初の燃焼室容積である燃焼開始時の燃焼室容積Ｖ０を採用している。
【００８４】
次に図１２のフローに移ると、ステップ１８１では図１０のステップ１１と同様に、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ３］、図５のステップ１３で算出されている燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出されている燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、図５のステップ１５で算出されている反応確率ＲＰＲＯＢＡ［％］を読み込み、さらにシリンダ新気量ＭＡＣＹＬ［ｇ］、目標当量比ＴＦＢＹＡ、内部不活性ガス量ＭＲＥＳ［ｇ］、外部不活性ガス量ＭＥＧＲ［ｇ］を読み込む。
【００８５】
ここで、図１には外部ＥＧＲ装置は示していないが、図１２に関する限り外部ＥＧＲ装置を備えているエンジンを前提として説明する。この場合に、外部不活性ガス量ＭＥＧＲは例えば公知の手法（特開平１０−１４１１５０号公報参照）を用いて算出すればよい。なお、図１に示す本実施形態のように外部ＥＧＲ装置を備えていないエンジンを対象とするときには外部不活性ガス量ＭＥＧＲ＝０で扱えば足りる。シリンダ新気量ＭＡＣＹＬ、内部不活性ガス量ＭＲＥＳの算出については図１４以降で後述する。
【００８６】
ステップ１８２、１８３は図１０のステップ１６３、１６４と同様である。すなわち、ステップ１８２で圧縮上死点時期における有効圧縮比Ｅｃ ２を算出する。有効圧縮比Ｅｃ ２も上記（１２）式の有効圧縮比Ｅｃと同様に無次元の値であり、次式に示すように燃焼室５の圧縮上死点時における容積ＶＴＤＣを燃焼室５の吸気弁閉時期における容積ＶＩＶＣで除した値である。
【００８７】
Ｅｃ ２＝ｆ１３（ＶＴＤＣ、ＶＩＶＣ）＝ＶＴＤＣ／ＶＩＶＣ…（４３）
（４３）式において燃焼室５の圧縮上死点時における容積ＶＴＤＣは運転条件によらず一定であり、予めエンジンコントローラ３１のメモリに格納しておけばよい。
【００８８】
ステップ１８３では吸気弁閉時期ＩＶＣから圧縮上死点に至る間の燃焼室５内の断熱圧縮による温度上昇率ＴＣＯＭＰ ２を次式に示すように有効圧縮比Ｅｃ ２に基づいて算出する。
【００８９】
2004332653.tif 000006
ただし、κ：比熱比、
図１１と同様の特性のテーブルを予めエンジンコントローラ３１のメモリに格納しておき、有効圧縮比Ｅｃ ２から当該テーブルを検索することにより温度上昇率ＴＣＯＭＰ ２を求めることも可能である。
【００９０】
ステップ１８４ではシリンダ新気量ＭＡＣＹＬ、目標当量比ＴＦＢＹＡ、内部不活性ガス量ＭＲＥＳ、外部不活性ガス量ＭＥＧＲから次式により燃焼室５の総ガス質量ＭＧＡＳ［ｇ］を算出する。
【００９１】
ＭＧＡＳ＝ＭＡＣＹＬ×（１＋ＴＦＢＹＡ／１４．７）＋ＭＲＥＳ＋ＭＥＧＲ…（４５）
（４５）式右辺の括弧内の「１」は新気分、「ＴＦＢＹＡ／１４．７」は燃料分である。
【００９２】
ステップ１８５ではこの燃焼室５の総ガス質量ＭＧＡＳと、シリンダ新気量ＭＡＣＹＬ、目標当量比ＴＦＢＹＡを用い、次式により混合気の燃焼による温度上昇量（燃焼上昇温度）ＴＢＵＲＮ［Ｋ］を算出する。
【００９３】
ＴＢＵＲＮ＝｛ＭＡＣＹＬ×ＴＦＢＹＡ／１４．７×ＢＲｋ×Ｑ｝／（Ｃｖ×ＭＧＡＳ）…（４６）
ただし、Ｑ：燃料の定発熱量、
ＢＲｋ：シリンダ内燃料の燃焼質量割合、
Ｃｖ：定積比熱、
（４６）式右辺の分子はシリンダ内燃料による発生総熱量［Ｊ］、分母は単位発生熱量当たりの温度上昇率［Ｊ／Ｋ］を意味している。すなわち、（４６）式は熱力学の公式に当てはめた近似式である。
【００９４】
ここで、シリンダ内燃料の燃焼質量割合ＢＲｋとしては予め実験等で適合しておく。簡易的には例えば６０％／２＝３０％を設定する。これは、本実施形態では燃焼質量割合ＢＲが約６０％に達するまでを燃焼期間として扱うので、そのちょうど中間の３０％をＢＲｋとして設定するものである。
【００９５】
燃料の定発熱量Ｑは燃料の種類により異なる値であるので、燃料の種類に応じ予め実験等で求めておく。定積比熱Ｃｖは２〜３の値であり予め実験等で代表値を適合しておく。ただし、混合気に対してその組成、温度から定積比熱Ｃｖの値を求めることで、一層の算出精度の向上が可能である。
【００９６】
ステップ１８６では、燃焼室５の圧縮上死点における温度ＴＴＤＣ［Ｋ］を、燃焼室５の吸気弁閉時期における温度ＴＩＮＩに圧縮上死点までの温度上昇率ＴＣＯＭＰ ２を乗じその乗算値に上記の燃焼上昇温度ＴＢＵＲＮを加算することで、つまり次式により算出する。
【００９７】
ＴＴＤＣ＝ＴＩＮＩ×ＴＣＯＭＰ ２＋ＴＢＵＲＮ…（４７）
ステップ１８７では、この燃焼室５の圧縮上死点における温度ＴＴＤＣと容積ＶＴＤＣ及び燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ、容積ＶＩＶＣ及び温度ＴＩＮＩから次式により燃焼室５の圧縮上死点における圧力ＰＴＤＣ［Ｋ］を算出する。
【００９８】
ＰＴＤＣ＝ＰＩＮＩ×ＶＩＶＣ×ＴＴＤＣ／（ＶＴＤＣ×ＴＩＮＩ）…（４８）
（４８）式は状態方程式を用いて得たものである。すなわち、吸気弁閉時期における圧力、容積及び温度（ＰＩＮＩ、ＶＩＶＣ、ＴＩＮＩ）を用いて次の状態方程式が成立する。
【００９９】
ＰＩＮＩ×ＶＩＶＣ＝ｎ×Ｒ×ＴＩＮＩ…（補２）
ただし、ｎ：モル数、
Ｒ：ガス定数、
圧縮上死点近傍では容積はほぼ等しいので、圧縮上死点での圧力、容積及び温度（ＰＴＤＣ、ＶＴＤＣ、ＴＴＤＣ）を用いて次の状態方程式が成立する。
【０１００】
ＰＴＤＣ×ＶＴＤＣ＝ｎ×Ｒ×ＴＴＤＣ…（補３）
この（補３）式と上記（補２）との両式からｎ×Ｒを消去しＰＴＤＣについて解くと、上記（４８）式が得られる。
【０１０１】
ステップ１８８では図１０のステップ１６８と同様にして、次式（公知）により、主燃焼期間における層流燃焼速度ＳＬ２［ｍ／ｓｅｃ］を算出する。
【０１０２】
2004332653.tif 000007
ただし、Ｔｓｔｄ    ：基準温度［Ｋ］、
Ｐｓｔｄ    ：基準圧力［Ｐａ］、
ＳＬｓｔｄ：基準温度Ｔｓｔｄと基準圧力Ｐｓｔｄにおける基準層流燃焼速度［ｍ／ｓｅｃ］、
ＴＴＤＣ：燃焼室５の圧縮上死点における温度［Ｋ］、
ＰＴＤＣ：燃焼室５の圧縮上死点における圧力［Ｐａ］、
（４９）式の解説は上記（１６）式と同様ある。すなわち、（４９）式の基準温度Ｔｓｔｄと基準圧力Ｐｓｔｄと基準層流燃焼速度ＳＬｓｔｄは実験により予め定められる値である。燃焼室５の通常の圧力である２ｂａｒ以上の圧力下では、（４９）式の圧力項（ＰＴＤＣ／Ｐｓｔｄ）−０．１６は小さな値となる。従って、圧力項（ＰＴＤＣ／Ｐｓｔｄ）−０．１６を一定値として、基準層流燃焼速度ＳＬｓｔｄを基準温度Ｔｓｔｄのみで規定することも可能である。よって、基準温度Ｔｓｔｄが５５０［Ｋ］で、基準層流燃焼速度ＳＬｓｔｄが１．０［ｍ／ｓｅｃ］で、圧力項が０．７である場合の圧縮上死点における温度ＴＴＤＣと層流燃焼速度ＳＬ２との関係は近似的に次式で定義することができる。
【０１０３】
2004332653.tif 000008
ステップ１８９期間では主燃焼期間におけるガス流動の乱れ強さＳＴ２を算出する。このガス流動の乱れ強さＳＴ２も初期燃焼期間におけるガス流動の乱れ強さＳＴ１と同様に、エンジン回転速度ＮＲＰＭの関数として次式で表すことができる。
【０１０４】
ＳＴ２＝ｆ１７（ＮＲＰＭ）＝Ｃ２×ＮＲＰＭ…（２０）
ただし、Ｃ２：定数、
乱れ強さＳＴ２を回転速度をパラメータとするテーブルから求めることも可能である。
【０１０５】
ステップ１９０では、層流燃焼速度ＳＬ２［ｍ／ｓｅｃ］と主燃焼期間におけるガス流動の乱れ強さＳＴ２とから、主燃焼期間における燃焼速度ＦＬＡＭＥ２［ｍ／ｓｅｃ］を次式により算出する。
【０１０６】
ＦＬＡＭＥ２＝ＳＬ２×ＳＴ２…（２１）
ただし、ＳＬ２：層流燃焼速度［ｍ／ｓｅｃ］、
（２１）式は（１８）式と同様、ガス乱れに伴う燃焼速度への寄与を考慮したものである。
【０１０７】
ステップ１９１では、主燃焼期間ＢＵＲＮ２［ｄｅｇ］を（１９）式に類似した次式で算出する。
【０１０８】
ＢＵＲＮ２＝｛（ＮＲＰＭ×６）×（ＢＲ２×ＶＴＤＣ）｝／（ＲＰＲＯＢＡ×ＡＦ２×ＦＬＡＭＥ２）…（２２）
ただし、ＡＦ２：火炎核の反応面積［ｍ２］
ここで、（２２）式右辺のＢＲ２は主燃焼期間の開始時期より終了時期までの燃焼質量割合の変化代である。初期燃焼期間の終了時期に燃焼質量割合ＢＲが２％になり、その後、主燃焼期間が開始し、燃焼質量割合ＢＲが６０％に達して主燃焼期間が終了すると考えているので、ＢＲ２＝６０％−２％＝５８％を設定している。ＡＦ２は火炎核の成長行程における平均の反応面積であり、（１９）式のＡＦ１と同様に、予め実験的に定めた固定値とする。
【０１０９】
主燃焼期間では圧縮上死点を挟んで燃焼室容積が変化する。つまり、主燃焼期間の開始時期と、主燃焼期間の終了時期のほぼ中央に圧縮上死点位置が存在するとみなすことができる。また、圧縮上死点付近ではクランク角が変化しても燃焼室容積があまり変化しない。そこで主燃焼期間での燃焼室容積としてはこの圧縮上死点での燃焼室容積ＶＴＤＣで代表させることとしている。
【０１１０】
図１３は基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１０、図１２のうち遅く実行されるフローに続けて実行する。
【０１１１】
ステップ４１では、図１０のステップ１７１で算出されている初期燃焼期間ＢＵＲＮ１、図１２のステップ１９１で算出されている主燃焼期間ＢＵＲＮ２、図５のステップ１６で算出されている点火時期無駄時間相当クランク角ＩＧＮＤＥＡＤ、図５のステップ１７で算出されている基準クランク角θＰＭＡＸを読み込む。
【０１１２】
ステップ４２では、初期燃焼期間ＢＵＲＮ１と主燃焼期間ＢＵＲＮ２の合計を燃焼期間ＢＵＲＮ［ｄｅｇ］として算出する。
【０１１３】
ステップ４３では次式により基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出する。
【０１１４】
ＭＢＴＣＡＬ＝ＢＵＲＮ−θＰＭＡＸ＋ＩＧＮＤＥＡＤ…（２３）
ステップ４４では、この基本点火時期ＭＢＴＣＡＬから点火無駄時間相当クランク角ＩＧＮＤＥＡＤを差し引いた値を燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ［ｄｅｇＢＴＤＣ］として算出する。
【０１１５】
このようにして算出した基本点火時期ＭＢＴＣＡＬは、点火時期指令値として点火レジスタに移され、実際のクランク角がこの点火時期指令値と一致したタイミングでエンジンコントローラ３１より一次電流を遮断する点火信号が点火コイル１３に出力される。
【０１１６】
また、今サイクルの点火時期指令値としてステップ４３で算出された基本点火時期ＭＢＴＣＡＬが用いられたとすると、次サイクルの点火時期になるまでの間、ステップ４４で算出された燃焼開始時期推定値ＭＢＴＣＡＬ Ｅが図１０のステップ１６２において用いられる。
【０１１７】
次に、図１４は燃焼室５内の内部不活性ガス率ＭＲＥＳＦＲを算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。このフローは上記図５のフローに先立って実行する。
【０１１８】
ステップ５１ではエアフローメータ３２の出力と目標当量比ＴＦＢＹＡを読み込む。ステップ５２ではエアフロメータ３２の出力に基づいて、燃焼室５に流入する新気量（シリンダ新気量）ＭＡＣＹＬを算出する。このシリンダ新気量ＭＡＣＹＬの算出方法については公知の方法を用いればよい（特開２００１−５００９１号公報参照）。
【０１１９】
ステップ５３では、燃焼室５内の内部不活性ガス量ＭＲＥＳを算出する。この内部不活性ガス量ＭＲＥＳの算出については、図１５のフローにより説明する。
【０１２０】
図１５（図１４ステップ５３のサブルーチン）においてステップ６１では、燃焼室５内の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬを算出する。この不活性ガス量ＭＲＥＳＣＹＬの算出についてはさらに図１６のフローにより説明する。
【０１２１】
図１６（図１５ステップ６１のサブルーチン）においてステップ７１では、排気弁閉時期ＥＶＣ［ｄｅｇＢＴＤＣ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、圧力センサ４６により検出される排気圧力ＰＥＸＨ［ｋＰａ］を読み込む。
【０１２２】
ここで、吸気弁閉時期ＩＶＣが吸気ＶＴＣ機構２７に与える指令値から既知であったように、排気弁閉時期ＥＶＣも排気ＶＴＣ機構２８に与える指令値から既知である。
【０１２３】
ステップ７２では燃焼室５の排気弁閉時期ＥＶＣにおける容積ＶＥＶＣを算出する。これは吸気弁閉時期ＩＶＣにおける容積ＶＩＶＣと同様に、排気弁閉時期をパラメータとするテーブルを検索することにより求めればよい。すなわち、排気弁ＶＴＣ機構２８を備える場合には、排気弁閉時期ＥＶＣから図２３に示すテーブルを検索することにより、燃焼室５の排気弁閉時期ＥＶＣにおける容積ＶＥＶＣを求めればよい。排気ＶＴＣ機構２８を備えないときには定数で与えることができる。
【０１２４】
また、図示しないが圧縮比を変化させる機構を有する場合には、圧縮比の変化量に応じた排気弁閉時期における燃焼室容積ＶＥＶＣをテーブルから求める。排気ＶＴＣ機構２８に加えて圧縮比を変化させる機構をも有する場合には、排気弁閉時期と圧縮比変化量とに応じたマップを検索することにより排気弁閉時期における燃焼室容積を求める。
【０１２５】
ステップ７３では、目標当量比ＴＦＢＹＡから図２４に示すテーブルを検索することにより、燃焼室５内の不活性ガスのガス定数ＲＥＸを求める。図２４に示すように、不活性ガスのガス定数ＲＥＸは目標当量比ＴＦＢＹＡが１．０のとき、つまり理論空燃比のとき最も小さく、これより大きくても小さくても大きくなる。
【０１２６】
ステップ７４では、排気温度ＴＥＸＨに基づいて燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣを推定する。簡単には排気温度ＴＥＸＨをそのままＴＥＶＣとおけばよい。なお、燃焼室５の排気弁閉時期における温度ＴＥＶＣは、インジェクタ２１の燃料噴射量に応じた熱量により変化するため、このような特性をも加味すれば、ＴＥＶＣの算出精度が向上する。
【０１２７】
ステップ７５では、排気圧力ＰＥＸＨに基づいて燃焼室５の排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣを算出する。簡単には排気圧力ＰＥＸＨをＰＥＶＣと置けばよい。
【０１２８】
ステップ７６では、燃焼室５の排気弁閉時期ＥＶＣにおける容積ＶＥＶＣ、排気弁閉時期ＥＶＣにおける温度ＴＥＶＣ、排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣ及び不活性ガスのガス定数ＲＥＸから、燃焼室５の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬを次式により算出する。
【０１２９】
ＭＲＥＳＣＹＬ＝（ＰＥＶＣ×ＶＥＶＣ）／（ＲＥＸ×ＴＥＶＣ）…（２４）
このようにして燃焼室５の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬの算出を終了したら図１５に戻り、ステップ６２で吸排気弁１５、１６のオーバーラップ（図では「Ｏ／Ｌ」と略記する）中に排気側から吸気側へ吹き返す不活性ガス量であるオーバーラップ中吹き返し不活性ガス量ＭＲＥＳＯＬを算出する。
【０１３０】
この不活性ガス量ＭＲＥＳＯＬの算出については図１７のフローにより説明する。
【０１３１】
図１７（図１５ステップ６２のサブルーチン）においてステップ８１では、吸気弁開時期ＩＶＯ［ｄｅｇＢＴＤＣ］と、排気弁閉時期ＥＶＣ［ｄｅｇＢＴＤＣ］、図１６のステップ７４で算出されている燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣを読み込む。
【０１３２】
ここで、吸気弁開時期ＩＶＯは、吸気弁閉時期ＩＶＣより吸気弁１５の開き角だけ前の時期となるので、吸気弁閉時期ＩＶＣより吸気弁１５の開き角（予め分かっている）とから求めることができる。
【０１３３】
ステップ８２では吸気弁開時期ＩＶＯと排気弁閉時期ＥＶＣとから、吸排気弁のオーバーラップ量ＶＴＣＯＬ［ｄｅｇ］を次式により算出する。
【０１３４】
ＶＴＣＯＬ＝ＩＶＯ＋ＥＶＣ…（２５）
例えば、吸気ＶＴＣ機構２７用アクチュエータへの非通電時に吸気弁開時期ＩＶＯが吸気上死点位置にあり、吸気ＶＴＣ機構２７用アクチュエータへの通電時に吸気弁開時期が吸気上死点より進角する特性であり、かつ排気ＶＴＣ機構２８用アクチュエータへの非通電時に排気弁閉時期ＥＶＣが排気上死点にあり、排気弁ＶＴＣ機構２８用アクチュエータへの通電時に排気弁閉時期ＥＶＣが排気上死点より進角する特性である場合には、ＩＶＯとＥＶＣの合計が吸排気弁のオーバーラップ量ＶＴＣＯＬとなる。
【０１３５】
ステップ８３では、吸排気弁のオーバーラップ量ＶＴＣＯＬから、図２５に示すテーブルを検索することによりオーバーラップ中の積算有効面積ＡＳＵＭＯＬを算出する。図２５に示すようにオーバーラップ中の積算有効面積ＡＳＵＭＯＬは吸排気弁のオーバーラップ量ＶＴＣＯＬが大きくなるほど大きくなる値である。
【０１３６】
ここで、図２６は、吸排気弁のオーバーラップ中の積算有効面積ＡＳＵＭＯＬの説明図であり、横軸はクランク角、縦軸は吸気弁１２と排気弁１５とのそれぞれの開口面積を示している。オーバーラップ中の任意の時点における有効開口面積は、排気弁開口面積と吸気弁開口面積とのうち小さい方とする。オーバーラップ中の全期間における積算有効面積ＡＳＵＭＯＬは、吸気弁１５及び排気弁１６が開いている期間の積分値（図中の斜線部）である。
【０１３７】
このようにオーバーラップ中積算有効面積ＡＳＵＭＯＬを算出することで、吸気弁１５と排気弁１６とのオーバーラップ量を１つのオリフィス（流出孔）であると近似することができ、排気系の状態と吸気系の状態とからこの仮想オリフィスを通過するガス流量を簡略的に算出し得る。
【０１３８】
ステップ８４では、目標当量比ＴＦＢＹＡと、燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣとから、図２７に示すマップを検索することにより、燃焼室５に残留する不活性ガスの比熱比ＳＨＥＡＴＲを算出する。図２７に示したように、燃焼室に残留する不活性ガスの比熱比ＳＨＥＡＴＲは目標当量比ＴＦＢＹＡが１．０の近傍にあるときが最も小さくなり、それより大きくても小さくても大きくなる。また、目標当量比ＴＦＢＹＡが一定の条件では、燃焼室５の排気弁閉時期ＥＶＣにおける温度ＴＥＶＣが高くなるほど小さくなる。
【０１３９】
ステップ８５では過給判定フラグＴＢＣＲＧ及びチョーク判定フラグＣＨＯＫＥを設定する。この過給判定フラグＴＢＣＲＧ及びチョーク判定フラグＣＨＯＫＥの設定については図１８のフローにより説明する。
【０１４０】
図１８（図１７ステップ８５のサブルーチン）においてステップ１０１では、吸気圧力センサ４４により検出される吸気圧力ＰＩＮと、図１６のステップ７５で算出されている燃焼室５の排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣを読み込む。
【０１４１】
ステップ１０２では、吸気圧力ＰＩＮと、燃焼室５の排気弁閉時期ＥＶＣにおける圧力ＰＥＶＣとから、次式により吸気排気圧力比ＰＩＮＢＹＥＸを算出する。
【０１４２】
ＰＩＮＢＹＥＸ＝ＰＩＮ／ＰＥＶＣ…（２６）
この吸気排気圧力比ＰＩＮＢＹＥＸは無名数であり、これと１をステップ１０３で比較する。吸気排気圧力比ＰＩＮＢＹＥＸが１以下の場合には過給無しと判断し、ステップ１０４に進んで過給判定フラグＴＢＣＲＧ（ゼロに初期設定）＝０とする。
【０１４３】
吸気排気圧力比ＰＩＮＢＹＥＸが１より大きい場合には過給有りと判断し、ステップ１０５へ進んで過給判定フラグＴＢＣＲＧ＝１とする。
【０１４４】
ステップ１０６では、図１４のステップ５１で読み込まれている目標当量比ＴＦＢＹＡから図２８に示すテーブルを検索することにより、混合気の比熱比ＭＩＸＡＩＲＳＨＲを求め、これをステップ１０７で不活性ガスの比熱比ＳＨＥＡＴＲと入れ換える。図２８に示したように、混合気の比熱比ＭＩＸＡＩＲＳＨＲは、目標当量比ＴＦＢＹＡが小さくなるほど大きくなる値である。
【０１４５】
ステップ１０６、１０７において、不活性ガスの比熱比ＳＨＥＡＴＲを混合気の比熱比ＭＩＸＡＩＲＳＨＲに置き換えるのは、ターボ過給や慣性過給等の過給時を考慮したものである。すなわち、過給時には吸排気弁のオーバーラップ中のガス流れが吸気系から排気系へ向かう（吹き抜ける）ので、この場合においては、上記の仮想オリフィスを通過するガスの比熱比を不活性ガスの比熱比から混合気の比熱比に変更することで、吹き抜けるガス量を精度良く推定し、内部不活性ガス量を精度良く算出するためである。
【０１４６】
ステップ１０８では、図１７のステップ８４または図１８のステップ１０６、１０７で算出している不活性ガスの比熱比ＳＨＥＡＴＲに基づき、最小と最大とのチョーク判定しきい値ＳＬＣＨＯＫＥＬ、ＳＬＣＨＯＫＥＨを次式により算出する。
【０１４７】
2004332653.tif 000009
これらのチョーク判定しきい値ＳＬＣＨＯＫＥＬ、ＳＬＣＨＯＫＥＨは、チョークする限界値を算出している。
【０１４８】
ステップ１０８において、（２７ａ）右辺、（２７ｂ）右辺の各累乗計算が困難な場合には、（２７ａ）、（２７ｂ）式の算出結果を、最小チョーク判定しきい値ＳＬＣＨＯＫＥＬのテーブルと最大チョーク判定しきい値ＳＬＣＨＯＫＥＨのテーブルとしてそれぞれエンジンコントローラ３１のメモリに予め記憶しておき、不活性ガスの比熱比ＳＨＥＡＴＲから当該テーブルを検索することにより求めてもよい。
【０１４９】
テップ１０９では、吸気排気圧力比ＰＩＮＢＹＥＸが、最小チョーク判定しきい値ＳＬＣＨＯＫＥＬ以上でかつ最大チョーク判定しきい値ＳＬＣＨＯＫＥＨ以下の範囲内にあるか否か、すなわちチョーク状態にないか否かを判定する。吸気排気圧力比ＰＩＮＢＹＥＸが範囲内にある場合にはチョーク無しと判断し、ステップ１１０に進んでチョーク判定フラグＣＨＯＫＥ（ゼロに初期設定）＝０とする。
【０１５０】
吸気排気圧力比Ｐ１ＮＢＹＥＸが範囲内にない場合にはチョーク有りと判断し、ステップ１１１に進んでチョーク判定フラグＣＨＯＫＥ＝１とする。
【０１５１】
このようにして過給判定フラグとチョーク判定フラグの設定を終了したら図１７に戻り、ステップ８６〜８８で次の４つの場合分けを行う。
【０１５２】
〈１〉過給判定フラグＴＢＣＲＧ＝０かつチョーク判定フラグＣＨＯＫＥ＝０のとき
〈２〉過給判定フラグＴＢＣＲＧ＝０かつチョーク判定フラグＣＨＯＫＥ＝０のとき
〈３〉過給判定フラグＴＢＣＲＧ＝０かつチョーク判定フラグＣＨＯＫＥ＝１のとき
〈４〉過給判定フラグＴＢＣＲＧ＝１かつチョーク判定フラグＣＨＯＫＥ＝０のとき
そして、上記〈１〉のときにはステップ８９に進んで、過給無しかつチョーク無し時のオーバーラップ中の平均吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ１を、上記〈２〉のときにはステップ９０に進んで過給無しかつチョーク有り時のオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ２を、上記〈３〉のときにはステップ９１に進んで過給有りかつチョーク無し時のオーバーラップ中の平均吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ３を、上記〈４〉のときにはステップ９２に進んで過給有りかつチョーク有り時の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ４をそれぞれ算出し、算出結果をオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。
【０１５３】
ここで、過給無しかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ１の算出について図１９のフローにより説明する
図１９（図１７ステップ８９のサブルーチン）においてステップ１２１では、図１６のステップ７３、７５で算出されている不活性ガスのガス定数ＲＥＸ、燃焼室５の排気弁閉時期における圧力ＰＥＶＣを読み込む。
【０１５４】
ステップ１２２では、不活性ガスのガス定数ＲＥＸと、図１７のステップ８１で読み込まれている燃焼室５の排気弁閉時期における温度ＴＥＶＣとに基づき、後述するガス流量の算出式に用いる密度項ＭＲＳＯＬＤを次式により算出する。
【０１５５】
ＭＲＳＯＬＤ＝ＳＱＲＴ｛１／（ＲＥＸ×ＴＥＶＣ）｝…（２８）
ここで、（２８）式右辺の「ＳＱＲＴ」はすぐ右のカッコ内の値の平方根を計算させる関数である。
【０１５６】
なお、密度項ＭＲＳＯＬＤの平方根計算が困難な場合は、（２８）式の算出結果をマップとしてエンジンコントローラ３１のメモリに予め記憶しておき、ガス定数ＲＥＸと燃焼室５の排気弁閉時期における温度ＴＥＶＣとからそのマップを検索することにより求めてもよい。
【０１５７】
ステップ１２３では、図１７のステップ８４で算出されている不活性ガスの比熱比ＳＨＥＡＴＲと、図１８のステップ１０２で算出されている吸気排気圧力比ＰＩＮＢＹＥＸとに基づき、後述するガス流量の算出式に用いる圧力差項ＭＲＳＯＬＰを次式により算出する。
【０１５８】
ＭＲＳＯＬＰ＝ＳＱＲＴ［ＳＨＥＡＴＲ／（ＳＨＥＡＴＲ−１）×｛ＰＴＮＢＹＥＸ＾（２／ＳＨＥＡＴＲ）−ＰＴＮＢＹＥＸ＾（（ＳＨＥＡＴＲ＋１）／ＳＨＥＡＴＲ）｝］…（２９）
ステップ１２４では、これら密度項ＭＲＳＯＬＤ、圧力差項ＭＲＳＯＬＰと、燃焼室５の排気弁閉時期における圧力ＰＥＶＣとから、過給無しかつチョーク無し時のオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ１を次式（ガス流量の算出式）により算出し、その算出値をステップ１２５でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。
【０１５９】
ＭＲＥＳＯＬｔｍｐ１＝１．４×ＰＥＶＣ×ＭＲＳＯＬＤ×ＭＲＳＯＬＰ…（３０）
次に、過給無しかつチョーク有り時の吹き返し不活性ガス流量の算出について図２０のフローにより説明する
図２０（図１７ステップ９０のサブルーチン）においてステップ１３１、１３２では、図１９のステップ１２１、１２２と同様にして、不活性ガスのガス定数ＲＥＸ、燃焼室５の排気弁閉時期における圧力ＰＥＶＣを読み込み、これらから前述の（２８）式により密度項ＭＲＳＯＬＤを算出する。
【０１６０】
ステップ１３３では、図１７のステップ８４で算出されている不活性ガスの比熱比ＳＨＥＡＴＲに基づき、チョーク時圧力差項ＭＲＳＯＬＰＣを次式により算出する。
【０１６１】
ＭＲＳＯＬＰＣ＝ＳＱＲＴ［ＳＨＥＡＴＲ×｛２／（ＳＨＥＡＴＲ＋１）｝＾｛（ＳＨＥＡＴＲ＋１）／〔ＳＨＥＡＴＲ−１）｝］…（３１）
なお、（３１）式の累乗計算と平方根計算とが困難な場合には、（３１）式の算出結果を、チョーク時圧力差項ＭＲＳＯＬＰＣのテーブルとしてエンジンコントローラ３１のメモリに予めに記憶しておき、不活性ガスの比熱比ＳＨＥＡＴＲからそのテーブルを検索することにより求めてもよい。
【０１６２】
ステップ１３４では、これら密度項ＭＲＳＯＬＤ、チョーク時圧力差項ＭＲＳＯＬＰＣと、燃焼室５の排気弁閉時期における圧力ＰＥＶＣとから、過給無しかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ２を次式により算出し、その算出値をステップ１３５でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。
【０１６３】
ＭＲＥＳＯＬｔｍｐ２＝ＰＥＶＣ×ＭＲＳＯＬＤ×ＭＲＳＯＬＰＣ…（３２）
次に、過給有りかつチョーク無し時の吹き返しガス流量の算出について図２１のフローにより説明する
図２１（図１７ステップ９１のサブルーチン）においてステップ１４１では、吸気圧力センサ４４により検出される吸気圧力ＰＩＮを読み込む。
【０１６４】
ステップ１４２では、図１８のステップ１０６、１０７で算出されている不活性ガスの比熱比ＳＨＥＡＴＲと、図１８のステップ１０２で算出されている吸気排気圧力比ＰＩＮＢＹＥＸとから、過給時圧力差項ＭＲＳＯＬＰＴを次式により算出する。
【０１６５】
ＭＲＳＯＬＰＴ＝ＳＱＲＴ［ＳＨＥＡＴＲ／（ＳＨＥＡＴＲ−１）×｛ＰＩＮＢＹＥＸ＾（−２／ＳＨＥＡＴＲ）−ＰＩＮＢＹＥＸ＾（−（ＳＨＥＡＴＲ＋１）／ＳＨＥＡＴＲ）｝］…（３３）
なお、（３３）式の累乗計算と平方根計算とが困難な場合は、（３３）式の算出結果を、過給時圧力差項ＭＲＳＯＬＰＴのマップとしてエンジンコントローラ３１のメモリに予め記憶しておき、不活性ガスの比熱比ＳＨＥＡＴＲと吸気排気圧力比ＰＩＮＢＹＥＸとからそのマップを検索することにより求めてもよい。
【０１６６】
ステップ１４３では、この過給時圧力差項ＭＲＳＯＬＰＴと吸気圧力ＰＩＮとに基づいて、過給有りかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ３を次式により算出し、その算出値をステップ１４４でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。
【０１６７】
ＭＲＥＳＯＬｔｍｐ３＝−０．１５２×ＰＩＮ×ＭＲＳＯＬＰＴ…（３４）
ここで、（３４）式の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ３は負の値とすることで、オーバーラップ中に吸気系から排気系へ吹き抜ける混合気のガス流量を表すことができる。
【０１６８】
次に、過給有りかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出について図２２のフローにより説明する
図２２（図１７ステップ９２のサブルーチン）においてステップ１５１、１５２では、図２１のステップ１４１と同じく吸気圧力センサ４４により検出される吸気圧力ＰＩＮを読み込むと共に、図２０のステップ１３２と同じくチョーク時圧力差項ＭＲＳＯＬＰＣを前述の（３１）式により算出する。
【０１６９】
ステップ１５３では、このチョーク時圧力差項ＭＲＳＯＬＰＣと吸気圧力ＰＩＮとに基づいて、過給有りかつチョーク有り時のオーバーラップ中吹き返しガス流量ＭＲＥＳＯＬｔｍｐ４を次式により算出し、その算出値をステップ１５４でオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐに移す。
【０１７０】
ＭＲＥＳＯＬｔｍｐ４＝−０．１０８×ＰＩＮ×ＭＲＳＯＬＰＣ…（３５）
ここで、（３５）式の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐ４も、ＭＲＥＳＯＬｔｍｐ３と同様、負の値とすることで、オーバーラップ中に吸気側から排気側へ吹き抜ける混合気のガス流量を表すことができる。
【０１７１】
このようにして、過給の有無とチョークの有無との組み合わせにより場合分けした、オーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐの算出を終了したら図１７に戻り、ステップ９３においてこのオーバーラップ中の吹き返し不活性ガス流量ＭＲＥＳＯＬｔｍｐとオーバーラップ期間中の積算有効面積ＡＳＵＭＯＬとから、次式によりオーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬを算出する。
【０１７２】
ＭＲＥＳＯＬ＝（ＭＲＥＳＯＬｔｍＰ×ＡＳＵＭＯＬ×６０）／（ＮＲＰＭ×３６０）…（３６）
このようにしてオーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬの算出を終了したら図１５に戻り、ステップ６３において燃焼室５内の排気弁閉時期ＥＶＣにおける不活性ガス量ＭＲＥＳＣＹＬと、このオーバーラップ中吹き返しガス量ＭＲＥＳＯＬとを加算して、つまり次式により内部不活性ガス量ＭＲＥＳを算出する。
【０１７３】
ＭＲＥＳ＝ＭＲＥＳＣＹＬ＋ＭＲＥＳＯＬ…（３７）
前述のように、過給有り時にはオーバーラップ中吹き返し不活性ガス流量（ＭＲＥＳＯＬｔｍｐ３、ＭＲＥＳＯＬｔｍｐ４）が負となるため、上記（３６）式のオーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬも負となり、このとき（３７）式によれば、オーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬの分だけ内部不活性ガス量が減じられる。
【０１７４】
このようにして内部不活性ガス量ＭＲＥＳの算出を終了したら図１４に戻り、ステップ５４においてこの内部不活性ガス量ＭＲＥＳと、目標当量比ＴＦＢＹＡとを用いて、次式により内部不活性ガス率ＭＲＥＳＦＲ（燃焼室５内の総ガス量に対する内部不活性ガス量の割合）を算出する。
【０１７５】
ＭＲＥＳＦＲ＝ＭＲＥＳ／｛ＭＲＥＳ＋ＭＡＣＹＬ×（１＋ＴＦＢＹＡ／１４．７）｝…（３８）
これで内部不活性ガス率ＭＲＥＳＦＲの算出を総て終了する。
【０１７６】
このように本実施形態によれば、内部不活性ガス量ＭＲＥＳを、燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬと、吸排気弁のオーバーラップ中の吹き返しガス量ＭＲＥＳＯＬとで構成し（図１５のステップ６３参照）、この場合に、燃焼室５の排気弁閉時期における温度ＴＥＶ及び圧力ＰＥＶＣを算出し（図１６のステップ７４、７５）、これら温度ＴＥＶＣ、圧力ＰＥＶＣと不活性ガスのガス定数ＲＥＸとに基づいて状態方程式（上記（２４）式）により燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬを算出する（図１６のステップ７６参照）ようにしたので、特に、燃焼室５内部の状態量（ＰＥＶＣ、ＶＥＶＣ、ＴＥＶＣ）が刻々と変化する過渡運転時においても、運転条件に関わらず精度良く燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬを算出（推定）できる。
【０１７７】
また、燃焼室５の排気弁閉時期における温度ＴＥＶＣ及び圧力ＰＥＶＣ、不活性ガスのガス定数ＲＥＸ及び比熱比ＳＨＥＡＴＲ、吸気圧力ＰＩＮに基づいてオーバーラップ中の吹き返し不活性ガス流量（ＭＲＥＳＯＬｔｍｐ１、ＭＲＥＳＯＬｔｍｐ２）を算出し（図１９、図２０参照）、このガス流量にオーバーラップ中の積算有効面積ＡＳＵＭＯＬを乗算して、オーバーラップ中の吹き返しガス量ＭＲＥＳＯＬを算出する（図１７のステップ９３参照）ようにしたので、精度良くオーバーラップ中吹き返しガス量ＭＲＥＳＯＬを算出（推定）できる。
【０１７８】
このように、燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬ、オーバーラップ中吹き返しガス量ＭＲＥＳＯＬとも精度良く算出（推定）できると、これらの和である内部不活性ガス量ＭＲＥＳも精度良く算出（推定）できることになり、この精度良く推定することが可能となった内部不活性ガス量ＭＲＥＳに基づいて算出される内部不活性ガス率ＭＲＥＳＦＲを、点火時期の算出に用いる燃焼室５内の吸気弁閉時期ＩＶＣにおける温度ＴＩＮＩに活かすことで（図５のステップ１３参照）、燃焼室５内の吸気弁閉時期ＩＶＣにおける温度ＴＩＮＩを精度良く算出できる。また、精度良く推定することが可能となった内部不活性ガス量ＭＲＥＳを、燃料噴射量、バルブ開閉タイミング（オーバーラップ量）などにも活かすことで、エンジンを適切に制御することが可能である。
【０１７９】
また、不活性ガスのガス定数ＲＥＸや不活性ガスの比熱比ＳＨＥＡＴＲは、目標当量比ＴＦＢＹＡに応じた値としているので（図２４、図２７参照）、理論空燃比を外れた空燃比での運転時（例えば理論空燃比よりもリーンな空燃比で運転を行うリーン運転時、冷間始動時のようにエンジンが元々不安定な状態を安定させるために理論空燃比の空燃比よりもリッチ側の空燃比で運転するエンジン始動直後、同じく大きな出力が要求されるために理論空燃比の空燃比よりもリッチ側の空燃比で運転する全負荷運転時）にも、燃焼室５の排気弁閉時期における不活性ガス量ＭＲＥＳＣＹＬ、オーバーラップ中吹き返しガス量ＭＲＥＳＯＬ、これらの合計である内部不活性ガス量ＭＲＥＳ、これに基づく内部不活性ガス率ＭＲＥＳＦＲを精度良く算出できる。
【０１８０】
また、オーバーラップ期間の積算有効面積ＡＳＵＭＯＬを仮想オリフィスの面積とし、この仮想オリフィスを排気が燃焼室５から吸気系へと吹き抜けると仮定しているので、オーバーラップ中の吹き返し不活性ガス量ＭＲＥＳＯＬの算出が簡略化されている。
【０１８１】
ここで、本実施形態の作用効果を説明する。
【０１８２】
本実施形態（請求項１に記載の発明）によれば、燃焼ガスの層流状態での燃焼速度である層流燃焼速度（ＳＬ１、ＳＬ２）、燃焼ガス体積に近似させた燃焼ガス体積相当容積（Ｖ０、ＶＴＤＣ）、基準クランク角θＰＭＡＸ（所定クランク角）までの燃焼質量割合の変化代（ＢＲ１、ＢＲ２）及び所定運転条件での燃焼ガスの燃焼のしやすさを示す反応確率ＲＰＲＯＢＡに基づいて燃焼開始から基準クランク角θＰＭＡＸまでの燃焼期間ＢＵＲＮ（＝ＢＵＲＮ１＋ＢＵＲＮ２）を算出するようにしたので、従来装置のように未燃ガス密度を用いることなく燃焼期間ＢＵＲＮを算出することが可能となり、これにより正確かつ容易にＭＢＴの得られる基本点火時期ＭＢＴＣＡＬを算出することができる。
【０１８３】
特に、初期燃焼期間ＢＵＲＮ１の算出に必要となる燃焼開始時容積Ｖ０を算出するには燃焼開始時期が必要となる。この場合に、点火時期と燃焼開始時期との間には密接な関係があるものの（図４参照）、これから点火時期を算出しようとしている段階で、点火時期の後にくる燃焼開始時期を知ることはできない。本実施形態（請求項１に記載の発明）によれば、この燃焼開始時容積Ｖ０（初期燃焼期間における燃焼ガス体積相当容積）を算出する際に、前回のサイクルで燃焼時に用いた点火時期から点火無駄時間相当クランク角ＩＧＮＤＥＡＤを差し引いた値を燃焼開始時期推定値ＭＢＴＣＡＬ Ｅとして用いるようにしているので、今回のサイクルでの点火時期制御上必要となる燃焼開始時容積Ｖ０を支障なく算出することができる。
【０１８４】
また、定常運転時であれば今回のサイクルと次回のサイクルとで点火時期、燃焼開始時期は変わらないとして点火時期より燃焼開始時期を推定できるのであり、燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを、基本点火時期ＭＢＴＣＡＬの１サイクル前の算出値に基づいて算出する本実施形態（請求項２に記載の発明）によれば、燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを簡易に得ることができる。
【０１８５】
基準クランク角θＰＭＡＸまでの燃焼質量割合は約６０％であることを新たに見出していることから、これに対応して燃焼期間ＢＵＲＮの終期を基準クランク角θＰＭＡＸとする本実施形態（請求項７に記載の発明）によれば、燃焼質量割合をわざわざ算出する必要がなく、これにより燃焼期間ＢＵＲＮの算出が容易になる。また、燃焼ガスによる燃焼圧力が最大となるときのクランク角である基準クランク角θＰＭＡＸはほぼ一定であるので、基準クランク角θＰＭＡＸを決定するために図示トルクもしくは図示平均有効圧力を算出する必要が無く演算負荷を軽減できる。
【０１８６】
図２９、図３０は基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出するための第２、第３実施形態の各フローチャートで、第１実施形態の図１３と置き換わるものである。図１３と同一部分には同一のステップ番号をつけている。
【０１８７】
第２実施形態の図２９から説明すると、第１実施形態の図１３と相違する部分はステップ２０１、２０２である。すなわち、ステップ２０１では、初期燃焼期間推定値ＢＵＲＮ１ Ｅを主燃焼期間ＢＵＲＮ２に比例させて、つまり次式により算出する。
【０１８８】
ＢＵＲＮ１＿Ｅ＝ａ１×ＢＵＲＮ２＋ｂ１…（５１）
ただし、ａ１、ｂ１：定数、
（５１）式は実験式である。運転条件を相違させて実験した結果、得られた初期燃焼期間と主燃焼期間のデータを概括すると、運転条件が相違しても初期燃焼期間と主燃焼期間との間に比例関係があるとみなせることから、前回のサイクルでの燃焼時の主燃焼期間ＢＵＲＮ２から今回のサイクルでの燃焼時の初期燃焼期間を推定するようにしたものである。（５１）式の定数ａ１、ｂ１の値はエンジン機種により相違するので適合等により決定する。
【０１８９】
ステップ２０２ではこの初期燃焼期間推定値ＢＵＲＮ１ Ｅを用いて次式により今回のサイクルでの燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ［ｄｅｇＢＴＤＣ］を算出する。
【０１９０】
ＭＢＴＣＡＬ Ｅ＝（ＢＵＲＮ１＿Ｅ＋ＢＵＲＮ２）−ＰＭＡＸ…（５２）
このように前回のサイクルでの主燃焼期間に比例させて今回のサイクルでの初期燃焼期間推定値を算出し、この今回のサイクルでの初期燃焼期間推定値と、主燃焼期間とに基づいて今回のサイクルでの燃焼開始時期推定値を算出する第２実施形態（請求項６に記載の発明）でも、第１実施形態と同様の作用効果が得られる。
【０１９１】
次に第３実施形態の図３０ではステップ２１１〜２１５が第１実施形態の図１３と相違する。すなわち、ステップ２１１では図１２のステップ１８１と同様に、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、シリンダ新気量ＭＡＣＹＬ［ｇ］、目標当量比ＴＦＢＹＡ、内部不活性ガス量ＭＲＥＳ［ｇ］、外部不活性ガス量ＭＥＧＲ［ｇ］を読み込み、これらを用いステップ２１２で上記の（４５）式により総ガス質量ＭＧＡＳ［ｇ］を、またステップ２１３で次式により燃料質量ＭＦＵＥＬ［ｇ］をそれぞれ算出する。
【０１９２】
ＭＦＵＥＬ＝ＭＡＣＹＬ×ＴＦＢＹＡ／１４．７…（５３）
ステップ２１４ではこれら総ガス質量ＭＧＡＳ、燃料質量ＭＦＵＥＬから次式により総ガス質量と燃料質量の比である総ガス燃料比ＧＢＹＦ［無名数］を算出する。
【０１９３】
ＧＢＹＦ＝（ＭＧＡＳ／ＭＦＵＥＬ）…（５４）
この比ＧＢＹＦは燃焼室５内の不活性ガス状態を表す指標である。すなわち、同じ燃料量のときにこの比ＧＢＹＦが大きいほど不活性ガス量が多いことを表す。
【０１９４】
ステップ２１５ではこの総ガス燃料比ＧＢＹＦとエンジン回転速度ＮＲＰＭから図３１を内容とするマップを検索することにより燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを求める。
【０１９５】
燃焼開始時期推定値ＭＢＴＣＡＬ Ｅは、図３１に示したように総ガス燃料比ＧＢＹＦが一定であればエンジン回転速度ＮＲＰＭが低いほど、またエンジン回転速度ＮＲＰＭが一定であれば総ガス燃料比ＧＢＹＦが大きいほど（つまり不活性ガス量が少ないほど、空燃比でいえばリーンになるほど）、進む特性である。これは、総ガス燃料比ＧＢＹＦが一定であれば低回転速度になるほど、またエンジン回転速度ＮＲＰＭが一定であれば総ガス燃料比ＧＢＹＦが大きいほど燃焼期間が長引くので、ＭＢＴが得られるようにするにはその分早めに燃やす、つまり燃焼開始時期を早める必要があるためである。図３１の実際の値は予め適合等で求めておく。
【０１９６】
第１実施形態のように、前回のサイクルでの点火時期に基づいて今回のサイクルでの燃焼開始時期を推定する場合、加速時のように過渡時になると前回のサイクルでの点火時期より推定される燃焼開始時期と、今回のサイクルでの実際の燃焼開始時期との間にずれが生じる。例えば、アクセルペダルをステップ的に踏み込む場合を考えると、踏み込む前と後では踏み込んだ後のほうが燃料が多くなり燃焼速度が速くなるので、ＭＢＴを得るための点火時期としてはアクセルペダルを踏み込む前のほうが踏み込んだ後より進角側にある。従って、アクセルペダルを踏み込んだ後には、実際の燃焼開始時期は踏み込む前より遅角側にくるはずである。このため、アクセルペダルを踏み込む前に算出した点火時期（つまり前回のサイクルでの点火時期）よりアクセルペダルを踏み込んだ後の燃焼開始時期（つまり今回のサイクルでの燃焼開始時期）を推定したのでは、その燃焼開始時期推定値が実際の燃焼開始時期より進角側となってしまい、実際の燃焼開始時期からずれ分だけ燃焼開始時容積Ｖ０（燃焼ガス体積相当容積）の算出値に誤差が生じ、ＭＢＴが得られなくなるのである。
【０１９７】
これに対して第３実施形態（請求項３に記載の発明）によれば、燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを、燃焼室５内の不活性ガス状態を表す指標である総ガス燃料比ＧＢＹＦとエンジン回転速度ＮＲＰＭとに基づいて算出するようにしているので、アクセルペダルを急激に踏み込むような加速時においても、燃焼開始時期を精度よく推定できる。
【０１９８】
図３２は基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出するための第４実施形態のフローチャートで、第３実施形態の図３０と置き換わるものである。図３０と同一部分には同一のステップ番号をつけている。
【０１９９】
第３実施形態では燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを、総ガス燃料比ＧＢＹＦとエンジン回転速度ＮＲＰＭとから図３１を内容とする１つのマップを用いて求めたが、第４実施形態はこのマップに代えて２つのテーブルを用いて求めるようにしたものである。すなわち、ステップ２２１ではエンジン回転速度ＮＲＰＭから図３３を内容とするテーブルを検索することにより、基本燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ１［ｄｅｇＢＴＤＣ］を、またステップ２２２では総ガス燃料比ＧＢＹＦから図３４を内容とするテーブルを検索することにより、不活性ガス燃料比補正係数ＭＢＴＣＡＬ Ｅ２［無名数］を求め、ステップ２２３においてこれらの積を燃焼開始時期推定値ＭＢＴＣＡＬ Ｅとして、つまり次式により燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを算出する。
【０２００】
ＭＢＴＣＡＬ Ｅ＝ＭＢＴＣＡＬ Ｅ１×ＭＢＴＣＡＬ Ｅ２…（５５）
所定の運転条件になると、前述の吸気ＶＴＣ機構２７や排気ＶＴＣ機構２８を働かせ、吸気弁１５や排気弁１６の開閉時期を変えることにより内部不活性ガスの量を増やしポンピングロスを減らして燃費をよくするのであるが、基本燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ１はこれら吸気ＶＴＣ機構２７や排気ＶＴＣ機構２８を働かせていないときの燃焼開始時期推定値である。なお、外部ＥＧＲ装置を備えるときには、この外部ＥＧＲ装置についても働かせていないときの燃焼開始時期推定値である。図３３に示したように基本燃焼開始時期推定値ＭＢＴＣＡＬ Ｅ１も図３１と同様にエンジン回転速度ＮＲＰＭが低いほど遅角する値である。
【０２０１】
一方、前述の吸気ＶＴＣ機構２７、排気ＶＴＣ機構２８、図示しない外部ＥＧＲ装置のいずれかを働かせるときには、３つとも働かせないときより燃焼室５内の不活性ガス量が増し燃焼期間が長引くことになるので、その分だけ燃焼開始時期を早める必要がある。補正係数ＭＢＴＣＡＬ Ｅ２はこうした不活性ガス量が燃焼期間つまり燃焼開始時期に影響する分を考慮するもので、図３４のように総ガス燃料比ＧＢＹＦが大きくなるほど（空燃比がリーンになるほど）大きくなる値である。図３３、図３４の実際の値は予め適合等で求めておく。
【０２０２】
このように燃焼開始時期推定値ＭＢＴＣＡＬ Ｅを２つのテーブルを用いて求めるようにした第４実施形態でも、第３実施形態と同様の作用効果が得られる。
【０２０３】
第３、第４実施形態で示したように、エンジン回転速度ＮＲＰＭと総ガス燃料比ＧＢＹＦに依存して変化する燃焼期間に基づけば燃焼開始時期を推定できることからすると、初期燃焼期間もこれら（ＮＲＰＭ、ＧＢＹＦ）に依存するので、図３５のようにエンジン回転速度ＮＲＰＭと総ガス燃料比ＧＢＹＦから初期燃焼期間推定値ＢＵＲＮ１ Ｅを求めることも可能である（第５実施形態）。従って、このものでは、図３５により得た初期燃焼期間推定値ＢＵＲＮ１ Ｅを第２実施形態の図２９ステップ２０２で用いればよい。なお、図３５の実際の値は予め適合等で求めておく。
【０２０４】
図３１、図３４、図３５においては総ガス燃料比ＧＢＹＦを用いたが、総ガス燃料比ＧＢＹＦに代えて、図１４ステップ５４で得られている内部不活性ガス率ＭＲＥＳＦＲを用いることもできる。図示しない外部ＥＧＲ装置をさらに備えるエンジンでは外部不活性ガス率をこの内部不活性ガス率ＭＲＥＳＦＲに加えた総不活性ガス率を用いればよい。外部不活性ガス率の算出方法は特開平１０−３１８１１０号公報に記載されている。
【０２０５】
実施形態では、初期燃焼期間を燃焼質量割合の変化代としてゼロから２％まで（つまりＢＲ１＝２％）、主燃焼期間を燃焼質量割合の変化代として２％から６０％まで（つまりＢＲ２＝５８％）と規定したが、本発明は必ずしもこの数値に限定されるものでない。
【０２０６】
実施形態では、燃焼ガス質量割合で説明したが、燃焼ガス質量そのものを用いてもかまわない。
【０２０７】
実施形態では、燃焼開始より所定クランク角までに燃焼室内で燃焼する燃焼ガス質量の変化代を、初期燃焼期間に相当する分と主燃焼期間に相当する分の２つに分割する場合で説明したが、これに限られるものでなく、所定クランク角までに燃焼室内で燃焼する燃焼ガス質量の変化代を３以上の複数に分割し、その分割されたそれぞれの燃焼ガス質量の変化代に対応する分割燃焼期間を算出し、それら全ての分割燃焼期間を合計した値を燃焼開始から所定クランク角までの燃焼期間として算出してもかまわない。このものによれば、点火からの燃焼ガス質量の変化が一様でない場合においても正確に燃焼期間を算出できる。
【０２０８】
請求項１に記載の発明において、層流燃焼速度算出手段の機能は図１０のステップ１６８、図１２のステップ１８８により、燃焼ガス体積相当容積算出手段の機能は図１０のステップ１６２により、燃焼ガス質量算出手段の機能は図１０のステップ１７１、図１２のステップ１９１により、反応確率算出手段の機能は図５のステップ１５により、燃焼期間算出手段の機能は図１０のステップ１７１、図１２のステップ１９１、図１３のステップ４２により、基本点火時期算出手段の機能は図１３のステップ４３により、燃焼開始時期推定値算出手段の機能は図１３のステップ４４により、燃焼ガス体積相当容積算出手段が燃焼ガス体積相当容積を算出する際にこの燃焼開始時期推定値を用いる点は図１０のステップ１６２によりそれぞれ果たされている。
【図面の簡単な説明】
【図１】一実施形態のエンジンの制御システム図。
【図２】エンジンコントローラで実行される点火時期制御のブロック図。
【図３】燃焼室の圧力変化図。
【図４】燃焼質量割合の変化を説明する特性図。
【図５】物理量の算出を説明するためのフローチャート。
【図６】エンジンのクランクシャフトとコネクティングロッドの位置関係を説明するダイアグラム。
【図７】水温補正係数の特性図。
【図８】当量比補正係数の特性図。
【図９】基準クランク角の特性図。
【図１０】初期燃焼期間の算出を説明するためのフローチャート。
【図１１】温度上昇率の特性図。
【図１２】主燃焼期間の算出を説明するためのフローチャート。
【図１３】基本点火時期の算出を説明するためのフローチャート。
【図１４】内部不活性ガス率の算出を説明するためのフローチャート。
【図１５】内部不活性ガス量の算出を説明するためのフローチャート。
【図１６】ＥＶＣ時不活性ガス量の算出を説明するためのフローチャート。
【図１７】オーバーラップ中吹き返し不活性ガス量の算出を説明するためのフローチャート。
【図１８】過給判定フラグ、チョーク判定フラグの設定を説明するためのフローチャート。
【図１９】過給無しかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。
【図２０】過給無しかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。
【図２１】過給有りかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。
【図２２】過給有りかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。
【図２３】排気弁閉時期における燃焼室容積の特性図。
【図２４】不活性ガスのガス定数の特性図。
【図２５】オーバーラップ中の積算有効面積の特性図。
【図２６】オーバーラップ中の積算有効面積の説明図。
【図２７】不活性ガスの比熱比の特性図。
【図２８】混合気の比熱比の特性図。
【図２９】第２実施形態の基本点火時期の算出を説明するためのフローチャート。
【図３０】第３実施形態の基本点火時期の算出を説明するためのフローチャート。
【図３１】第３実施形態の燃焼開始時期推定値の特性図。
【図３２】第４実施形態の基本点火時期の算出を説明するためのフローチャート。
【図３３】第４実施形態の基本燃焼開始時期推定値の特性図。
【図３４】第４実施形態の総ガス燃料比補正係数の特性図。
【図３５】第５実施形態の初期燃焼期間推定値の特性図。
【符号の説明】
１  エンジン
５  燃焼室
１１  点火装置（火花点火手段）
１５  吸気弁
２１  燃料インジェクタ
２７  吸気ＶＴＣ機構
３１  エンジンコントローラ
３３、３４  クランク角センサ
４３  吸気温度センサ
４４  吸気圧力センサ
４５  排気温度センサ
４６  排気圧力センサ
【図１】
2004332653.tif 000010
【図２】
2004332653.tif 000011
【図３】
2004332653.tif 000012
【図４】
2004332653.tif 000013
【図５】
2004332653.tif 000014
【図６】
2004332653.tif 000015
【図７】
2004332653.tif 000016
【図８】
2004332653.tif 000017
【図９】
2004332653.tif 000018
【図１０】
2004332653.tif 000019
【図１１】
2004332653.tif 000020
【図１２】
2004332653.tif 000021
【図１３】
2004332653.tif 000022
【図１４】
2004332653.tif 000023
【図１５】
2004332653.tif 000024
【図１６】
2004332653.tif 000025
【図１７】
2004332653.tif 000026
【図１８】
2004332653.tif 000027
【図１９】
2004332653.tif 000028
【図２０】
2004332653.tif 000029
【図２１】
2004332653.tif 000030
【図２２】
2004332653.tif 000031
【図２３】
2004332653.tif 000032
【図２４】
2004332653.tif 000033
【図２５】
2004332653.tif 000034
【図２６】
2004332653.tif 000035
【図２７】
2004332653.tif 000036
【図２８】
2004332653.tif 000037
【図２９】
2004332653.tif 000038
【図３０】
2004332653.tif 000039
【図３１】
2004332653.tif 000040
【図３２】
2004332653.tif 000041
【図３３】
2004332653.tif 000042
【図３４】
2004332653.tif 000043
【図３５】
2004332653.tif 000044