JP4120620B2 - エンジンの制御装置 - Google Patents

Description

本発明は、エンジン（内燃機関）の軸トルク演算装置及びエンジンの制御装置に関する。

燃焼室内の圧力を検出する圧力センサの出力に基づいて図示トルクを算出するものがある（特許文献１参照）。
特開平８−３１２４２６号公報

ところで、特許文献１に記載の技術では、燃焼室内の圧力を検出する圧力センサが必要であり、コストがかかる。

そこで本発明は、圧力センサにより燃焼室内の圧力を直接検出することなくエンジンの発生する軸トルクを推定可能な軸トルク演算装置を提供すると共に、その推定した軸トルクを用いたエンジンの制御装置を提供することを目的とする。

本発明は、エンジンの運転状態を検出し、この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合を推定し、この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定すると共に、燃焼中のクランク角における燃焼質量割合の近似式を算出し、この燃焼質量割合の近似式に基づいて図示正仕事推定値を算出し、この図示正仕事推定値にポンピングロスと摩擦損失を加算した値を軸トルク推定値として算出するように構成する。

本発明は、これに加えて前記推定した軸トルクに基づいてトルク制御を行うようにも構成する。

本発明によれば、エンジンの運転状態を検出し、この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合（Ｘ１’、Ｘ２’）を推定し、この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定すると共に、燃焼中のクランク角における燃焼質量割合の近似式を算出し、この燃焼質量割合の近似式に基づいて図示正仕事推定値を算出し、この図示正仕事推定値にポンピングロスと摩擦損失を加算した値を軸トルク推定値として算出するように構成したので、圧力センサにより燃焼室内の圧力を直接検出することなくエンジンの発生する軸トルクを推定できるため、安価な軸トルク検出装置を実現できる。

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

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

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

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

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

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

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

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

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

この場合に、本発明では、圧力センサにより燃焼室内の圧力を直接検出することなくエンジンの発生する軸トルクを推定できるようにするため、燃焼室内のガスの燃焼中のクランク角における燃焼質量割合を推定し、この推定した燃焼質量割合に基づいてエンジンの発生する軸トルク推定値を算出する。

さらに本発明では、この軸トルク推定値に基づいてトルク制御（吸入空気量の増加と点火時期の遅角補正）を行う。

ここでのトルク制御を簡単に説明しておくと、点火は基本的にＭＢＴの得られる点火時期で行うのであるが、点火時期がＭＢＴの得られる点火時期よりリタードされることがあり、このときに発生する軸トルクは、ＭＢＴの得られる点火時期で点火したときの軸トルクより低下する。そこで、ＭＢＴの得られる点火時期よりリタードされた点火時期での軸トルクを大きくしてＭＢＴの得られる点火時期での軸トルクへと近づけるため吸入空気量を増加させる。

一方、アクセル開度とエンジン回転速度からドライバの要求する要求軸トルクが定まっており、ＭＢＴの得られる点火時期での軸トルクがこの要求軸トルクを超えるようだと軸トルクの出過ぎであり、運転上の違和感を生じさせかねない。そこで、ＭＢＴの得られる点火時期での軸トルクが要求軸トルクを超えるときには点火時期をリタードさせる。

図２は点火時期制御の全体の流れを示すフローである。このフローは一定時間毎に実行するフローではなく操作の流れを示している。

ステップ１では基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を、ステップ２〜７では点火時期リタード量ＲＥＴ(new)［ｄｅｇ］をそれぞれ算出する。

ここでは、基本点火時期ＭＢＴＣＡＬの算出を先に説明する。まず、燃焼解析に基づく点火時期制御を概説する（基本的な考え方は特開２００３−１４８２３６公報に記載されている）。

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

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

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

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

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

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

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

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

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

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

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

内部不活性ガス率ＭＲＥＳＦＲは燃焼室内に残留する不活性ガス量を燃焼室内の総ガス量で除した値で、その算出については後述する。

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

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

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

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

ＶＩＶＣ＝ｆ１（θ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
…（５０）
ステップ１８９では主燃焼期間におけるガス流動の乱れ強さＳＴ２を算出する。このガス流動の乱れ強さＳＴ２も初期燃焼期間におけるガス流動の乱れ強さＳＴ１と同様に、エンジン回転速度ＮＲＰＭの関数として次式で表すことができる。

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

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

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

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

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

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

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

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

上記の点火無駄時間相当クランク角ＩＧＮＤＥＡＤの算出については図２９のフローにより説明する（詳細は特願２００３−１０９０４０を参照）。図２９のフローは、図１０、図１２のフローで示した燃焼期間ＢＵＲＮ１、ＢＵＲＮＮ２を算出した後、図１３のフローの実行前に実行する。

ここでは、点火無駄時間ＤＥＡＤＴＩＭＥ［μｓｅｃ］をまず算出し、その点火無駄時間ＤＥＡＤＴＩＭＥをクランク角に換算する。点火無駄時間ＤＥＡＤＴＩＭＥは点火指令を受けてから燃焼室５内で実際に燃焼が開始するまでの遅れ時間である。点火無駄時間に影響する因子は燃焼室５の圧力及び温度、混合気の空燃比の３つでありしかもこれら各因子の点火無駄時間への影響は互いに独立であるとして算出する。

図２９においてステップ２０１では図１３のステップ４３で算出されている基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、図５のステップ１３で算出されている燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出されている燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、目標当量比ＴＦＢＹＡ、エンジン回転速度ＮＲＰＭ［ｒｐｍ］を読み込む。

ステップ２０２〜２０５では図１０のステップ１６２〜１６７と同様にして燃焼開始時期における有効圧縮比Ｅｃ（＝Ｖ０／ＶＩＶＣ）を算出し、その圧縮比Ｅｃを用いて燃焼室５の燃焼開始時期における温度Ｔ０と、燃焼室５の燃焼開始時期における圧力Ｐ０とを算出する。

Ｔ０＝ＴＩＮＩ×Ｅｃ＾（κ−１） …（１０−１）
Ｐ０＝ＰＩＮＩ×Ｅｃ＾κ …（１０−２）
ただし、κ：比熱比、
ステップ２０６では燃焼室５の燃焼開始時期における圧力Ｐ０の関数として点火無駄時間基本値ＤＥＡＤＴＩＭＥ０［μｓｅｃ］を次式により算出する。

ＤＥＡＤＴＩＭＥ０＝ｆ２１（Ｐ０） …（１０−３）
例えば燃焼室５の燃焼開始時期における圧力Ｐ０から図３０を内容とするテーブルを検索することにより点火無駄時間基本値ＤＥＡＤＴＩＭＥ０を求めればよい。点火無駄時間基本値ＤＥＡＤＴＩＭＥ０は理論空燃比のときかつ燃焼室５の燃焼開始時期における温度が基準温度である場合の値で、図３０のように圧力Ｐ０が高いほど大きくなる。高圧で点火無駄時間が長くなるのは、高圧下では混合気が気化しにくくなるためである。

ステップ２０７〜２０９は点火無駄時間に対する２つの補正係数を算出する部分である。すなわち、燃焼室５の燃焼開始時期における温度が基準温度より低下したり、混合気の空燃比が理論空燃比より外れるときには、点火無駄時間基本値ＤＥＡＤＴＩＭＥ０が実際と合わなくなるので、温度補正係数Ｋｄｔｍｔと空気過剰率補正係数Ｋｄｔｍｌｍｂとを導入している。

まずステップ２０７では燃焼室５の燃焼開始時期における温度Ｔ０の関数として点火無駄時間の温度補正係数Ｋｄｔｍｔを次式により算出する。

Ｋｄｔｍｔ＝ｆ２２（Ｔ０） …（１０−４）
例えば燃焼室５の燃焼開始時期における温度Ｔ０から図３１を内容とするテーブルを検索することにより温度補正係数Ｋｄｔｍｔを求めればよい。温度補正係数Ｋｄｔｍｔは図３１のように基準温度以上のとき１．０であり、基準温度より低くなるほど大きくなる値である。すなわち、温度Ｔ０が低いほど点火無駄時間は長くなる。

ステップ２０８では目標当量比ＴＦＢＹＡの逆数を目標空気過剰率ＴＬＡＭＢＤＡとしてつまり次式により目標空気過剰率ＴＬＡＭＢＤＡ求める。

ＴＬＡＭＢＤＡ＝１／ＴＦＢＹＡ …（１０−５）
ステップ２０９では、この目標空気過剰率ＴＬＡＭＢＤＡの関数として点火無駄時間の空気過剰率補正係数Ｋｄｔｍｌｍｂを次式により算出する。

Ｋｄｔｍｌｍｂ＝ｆ２３（ＴＬＡＭＢＤＡ） …（１０−６）
例えば目標空気過剰率ＴＬＡＭＢＤＡから図３２を内容とするテーブルを検索することにより空気過剰率補正係数Ｋｄｔｍｌｍｂを求めればよい。空気過剰率補正係数Ｋｄｔｍｌｍｂは図３２のように目標空気過剰率ＴＬＡＭＢＤＡが１のとき（つまり理論空燃比のとき）最小の１．０となり、これより大きくても小さくても大きくなる値である。すなわち、空燃比が理論空燃比よりリッチ側にずれてもリーン側にずれても点火無駄時間は長くなる。

ステップ２１０では点火無駄時間基本値ＤＥＡＤＴＩＭＥ０を２つの補正係数Ｋｄｔｍｔ、Ｋｄｔｍｌｍｂで補正することにより、つまり次式により点火無駄時間ＤＥＡＤＴＩＭＥ［μｓｅｃ］を算出する。

ＤＥＡＤＴＩＭＥ＝ＤＥＡＤＴＩＭＥ０×Ｋｄｔｍｔ×Ｋｄｔｍｌｍｂ
…（１０−７）
ここでは、圧力Ｐ０によるものを基本値（ＤＥＡＤＴＩＭＥ０）として温度Ｔ０と空燃比（ＴＬＡＭＢＤＡ）により補正する形式としたが、これに限られるものでなく、温度Ｔ０によるものを基本値として圧力Ｐ０と空燃比（ＴＬＡＭＢＤＡ）により補正したり、空燃比（ＴＬＡＭＢＤＡ）によるものを基本値として温度Ｔ０と圧力Ｐ０により補正するようにしてもかまわない。また、点火無駄時間に影響する因子はこれに限られるものでなく、充填効率や不活性ガス率（内部不活性ガス率ＭＲＥＳや外部不活性ガス率）があるので、これらの因子をも考慮して点火無駄時間を算出することができる。ただし、５つの因子全てを考慮する必要は必ずしもなく、少なくとも１つの因子をパラメータとして点火無駄時間を算出すればよい。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

これで燃焼室５内の内部不活性ガス率ＭＲＥＳＦＲの算出についての説明を終了し、次にはエンジンの発生する軸トルク推定値の算出について説明する。

図２のステップ２では軸トルク推定値Ｔｒｑｅｓｔ１、Ｔｒｑｅｓｔ２を算出する。ここでは、基本点火時期ＭＢＴＣＡＬで点火したときの軸トルク推定値と、基本点火時期ＭＢＴＣＡＬよりリタード（遅角）された点火時期で点火したときの軸トルク推定値との２種類の軸トルク推定値を算出するので、基本点火時期ＭＢＴＣＡＬで点火したときの軸トルク推定値を第１軸トルク推定値Ｔｒｑｅｓｔ１、基本点火時期ＭＢＴＣＡＬよりリタードされた点火時期で点火したときの軸トルク推定値を第２軸トルク推定値Ｔｒｑｅｓｔ２として区別する。

また、各軸トルク推定値を求める過程で初期燃焼期間ＢＵＲＮ１、点火無駄時間相当クランク角ＩＧＮＤＥＡＤ、主燃焼期間ＢＵＲＮ２を算出するが、これらも第１軸トルク推定値Ｔｒｑｅｓｔ１を求める過程で算出するＢＵＲＮ１、ＩＧＮＤＥＡＤ、ＢＵＲＮ２に「ｒｅｆ１」を付して「ＢＵＲＮ１ｒｅｆ１」、「ＩＧＮＤＥＡＤｒｅｆ１」、「ＢＵＲＮ２ｒｅｆ１」とし、これらをそれぞれ第１初期燃焼期間、第１点火無駄時間相当クランク角、第１主燃焼期間という。また、第２軸トルク推定値Ｔｒｑｅｓｔ２を求める過程で算出するＢＵＲＮ１、ＩＧＮＤＥＡＤ、ＢＵＲＮ２に「ｒｅｆ２」を付して「ＢＵＲＮ１ｒｅｆ２」、「ＩＧＮＤＥＡＤｒｅｆ２」、「ＢＵＲＮ２ｒｅｆ２」とし、これらをそれぞれ第２初期燃焼期間、第２点火無駄時間相当クランク角、第２主燃焼期間という。

さらに、各軸トルク推定値を求める過程で算出する燃焼質量割合の近似式、図示正仕事推定値をそれぞれ算出するが、これらも第１軸トルク推定値Ｔｒｑｅｓｔ１を求める過程で算出する燃焼質量割合の近似式、図示正仕事推定値にそれぞれ「１」を付して「Ｘ１’」、「Ｗｐｏｓ１」とし、これらをそれぞれ燃焼質量割合の第１近似式、第１図示正仕事推定値という。また、第２軸トルク推定値Ｔｒｑｅｓｔ２を求める過程で算出する燃焼質量割合の近似式、図示正仕事推定値にそれぞれ「２」を付して「Ｘ２’」、「Ｗｐｏｓ２」とし、これらをそれぞれ燃焼質量割合の第２近似式、第２図示正仕事推定値という。

上記の第１軸トルク推定値Ｔｒｑｅｓｔ１の算出については図３３（Ａ）のフローにより、上記の第２軸トルクＴｒｑｅｓｔ２の算出については図３３（Ｂ）のフローにより説明する。この場合に、点火時期が相違するだけで軸トルク推定値の算出方法そのものは変わらないので、第１軸トルク推定値Ｔｒｑｅｓｔ１の算出について詳述し、その後に第２軸トルク推定値Ｔｒｑｅｓｔ２の算出に言及する。

図３３（Ａ）（図２のステップ２のサブルーチン）においてステップ２２１では図１３のステップ４３で算出されている基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を読み込む。

ステップ２２２では点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］の値を点火時期前回値θｉｇｎｚ［ｄｅｇＢＴＤＣ］に移すと共に、基本点火時期ＭＢＴＣＡＬを点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］に入れる。

ステップ２２３では第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１、第１主燃焼期間ＢＵＲＮ２ｒｅｆ１、第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１の各値を第１初期燃焼期間前回値ＢＵＲＮ１ｒｅｆ１ｚ、第１主燃焼期間前回値ＢＵＲＮ２ｒｅｆ１ｚ、第１点火無駄時間相当クランク角前回値ＩＧＮＤＥＡＤｒｅｆ１ｚに移す。

ここで、ＢＵＲＮ１ｒｅｆ１、ＢＵＲＮ２ｒｅｆ１、ＩＧＮＤＥＡＤｒｅｆ１の値は後述するステップ２２４〜２２６で算出される値である。従って、初回はこれらＢＵＲＮ１ｒｅｆ１、ＢＵＲＮ２ｒｅｆ１、ＩＧＮＤＥＡＤｒｅｆ１は算出されていないので、ＢＵＲＮ１ｒｅｆ１ｚには図１０のステップ１７１により算出されているＢＵＲＮ１の値を、ＢＵＲＮ２ｒｅｆ１ｚには図１２のステップ１９１により算出されているＢＵＲＮ２の値を、ＩＧＮＤＥＡＤｒｅｆ１ｚには図２９のステップ２１１により算出されているＩＧＮＤＥＡＤの値を入れておく。

ステップ２２４〜２２６では第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１、第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１、第１主燃焼期間ＢＵＲＮ２ｒｅｆ１を算出する。これらの算出については図３４、図３５、図３６のフローにより説明する。

まず、図３４から説明すると、第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１の算出方法の基本的な考え方は図１０に示した初期燃焼期間の算出方法と同じである。

図３４（図３３（Ａ）のステップ２２４のサブルーチン）においてステップ２４１では図１４のステップ５２、５３で算出されているシリンダ新気量ＭＡＣＹＬ［ｋｇ］及び内部不活性ガス量ＭＲＥＳ［ｋｇ］、燃料量ＱＩＮＪ［ｋｇ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、温度センサ４５より検出される排気温度ＴＥＸＨ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、図５のステップ１５で算出されている反応確率ＲＰＲＯＢＡ［％］、図３３（Ａ）のステップ２２２、２２３で得られている点火時期前回値θｉｇｎｚ［ｄｅｇＢＴＤＣ］、第１初期燃焼期間前回値ＢＵＲＮ１ｒｅｆ１ｚ［ｄｅｇ］、第１主燃焼期間前回値ＢＵＲＮ２ｒｅｆ１ｚ［ｄｅｇ］、第１点火無駄時間相当クランク角前回値ＩＧＮＤＥＡＤＲｆ１ｚ［ｄｅｇ］を読み込む。

ステップ２４２では点火時期前回値θｉｇｎｚを用いて燃焼室５の燃焼質量割合１％時における容積Ｖｍｂ１ｐ［ｍ³］を次式により算出する。

Ｖｍｂ１ｐ＝ｆ２５（θｉｇｎｚ−ＩＧＮＤＥＡＤｒｅｆ１ｚ
−ＢＵＲＮ１ｒｅｆ１ｚ／２）
…（６１）
ステップ２４３では燃焼質量割合１％時の有効圧縮比Ｅｍｂ１ｐを次式により算出する。

Ｅｍｂ１ｐ＝Ｖｍｂ１ｐ／ＶＩＶＣ …（６２）
ステップ２４４ではシリンダ新気量ＭＡＣＹＬ［ｋｇ］をＷＩＤＲＹ［ｋｇ］に、内部不活性ガス量ＭＲＥＳ［ｋｇ］をＭＡＳＳＺ［ｋｇ］に移す。これらＷＩＤＲＹ、ＭＡＳＳＺは、軸トルク推定値の算出にのみ用いるために導入したもので、ＷＩＤＲＹはシリンダ新気量、ＭＡＳＳＺは内部不活性ガス量である。

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

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

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

ステップ２４７、２４８では燃焼質量割合１％時の瞬間圧縮比であるＥｍｂ１ｐと、燃焼質量割合１％時の燃焼質量割合であるＢＲ１ｐとを用いて、燃焼室５内の燃料が燃焼したときの平均温度ＴＣ［Ｋ］と平均圧力ＰＣ［Ｐａ］を次式により算出する。

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

ステップ２４９〜２５２により第１初期燃焼期間を算出するための操作は図１０のステップ１６８〜１７１により初期燃焼期間を算出するための操作と同様である。異なるのは図１０のステップ１６８において燃焼室５の燃焼開始時期における温度Ｔ０と圧力Ｐ０とを用いて層流燃焼速度ＳＬ１を算出したのに対して、図３４のステップ２４９では上記（６４）式、（６５）式で求めた温度ＴＣと圧力ＰＣを用いて層流燃焼速度ＳＬ１を算出する点、また図１０のステップ１７１において初期燃焼期間ＢＵＲＮ１の終了時期における燃焼質量割合であるＢＲ１（＝２％）と、燃焼室５の燃焼開始時期における容積であるＶ０とを用いて初期燃焼期間ＢＵＲＮ１を算出したのに対して、図３４のステップ２５２では初期燃焼期間の中間時期における燃焼質量割合である１％（＝ＢＲ１ｐ）と、燃焼室５の燃焼質量割合１％時における容積であるＶｍｂ１ｐとを用いて第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１を算出する点である。すなわち、ステップ２４９では第１初期燃焼期間における層流燃焼速度ＳＬ１［ｍ／ｓｅｃ］を、またステップ２５２では第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１［ｄｅｇ］をそれぞれ次式により算出する。

ＳＬ１＝ＳＬstd×（ＴＣ／Ｔstd）^2.18×（ＰＣ／Ｐstd）^-0.16 …（６６）
ただし、Ｔstd ：基準温度［Ｋ］、
Ｐstd ：基準圧力［Ｐａ］、
ＳＬstd：基準温度Ｔstdと基準圧力Ｐstdにおける基準層流燃焼
速度［ｍ／ｓｅｃ］、
ＢＵＲＮ１ｒｅｆ１＝｛（ＮＲＰＭ×６）×ＢＲ１ｐ×Ｖｍｂ１ｐ｝
／（ＲＰＲＯＢＡ×ＡＦ１×ＦＬＡＭＥ１） …（６７）
ただし、ＡＦ１ ：火炎核の反応面積（固定値）［ｍ²］、
図３５（図３３（Ａ）のステップ２２５のサブルーチン）のフローは第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤを算出するためのもので、その算出方法は前述した図２９に示す点火無駄時間相当クランク角ＩＧＮＤＥＡＤの算出方法と基本的に変わらない。

図２９に示す点火無駄時間相当クランク角ＩＧＮＤＥＡＤの算出方法と異なるのは、図２９においては燃焼室５の燃焼開始時期における容積Ｖ０を算出するのに前回燃焼開始時期ＭＢＴＣＹＣＬを用いたのに対して（図２９のステップ２０１、２０２参照）、図３５では点火時期θｉｇｎを用いる（図３５のステップ２６１、２６２参照）点だけである。従って異なる点を主に説明すると、図３５のステップ２６１において図３３（Ａ）のステップ２２２で得られている点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、図５のステップ１３で算出されている燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出されている燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、目標当量比ＴＦＢＹＡ、エンジン回転速度ＮＲＰＭ［ｒｐｍ］を読み込み、このうち点火時期θｉｇｎを用いステップ２６２において燃焼室５の燃焼開始時期における容積Ｖ０［ｍ³］を算出する。

ステップ２６３〜２６５では図２９のステップ２０３〜２０５と同様にして燃焼室５の燃焼開始時期における温度Ｔ０［Ｋ］と、燃焼室５の燃焼開始時期における圧力Ｐ０［Ｐａ］とを算出する。

図３６（図３３（Ａ）のステップ２２６のサブルーチン）のフローは第１主燃焼期間を算出するためのもので、ここでの第１主燃焼期間の算出方法の基本的な考え方は図１２に示した主初期燃焼期間の算出方法と同じである。

ステップ２８１では図１４のステップ５２、５３で算出されているシリンダ新気量ＭＡＣＹＬ［ｋｇ］及び内部不活性ガス量ＭＲＥＳ［ｋｇ］、燃料量ＱＩＮＪ［ｋｇ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、温度センサ４５より検出される排気温度ＴＥＸＨ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、図５のステップ１２で算出されている燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、エンジン回転速度ＮＲＰＭ［ｒｐｍ］、図５のステップ１５で算出されている反応確率ＲＰＲＯＢＡ［％］、図３３（Ａ）のステップ２２２で得られている点火時期θｉｇｎ、点火時期前回値θｉｇｎｚ、図３３（Ａ）のステップ２２３で得られている第１初期燃焼期間前回値ＢＵＲＮ１ｒｅｆ１ｚ［ｄｅｇ］、第１主燃焼期間前回値ＢＵＲＮ２ｒｅｆ１ｚ［ｄｅｇ］、第１点火無駄時間相当クランク角前回値ＩＧＮＤＥＡＤｒｅｆ１ｚ［ｄｅｇ］、図３４のステップ２５２で算出されている第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１［ｄｅｇ］、図３５のステップ２７１で算出されている第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１［ｄｅｇ］を読み込む。

ステップ２８２〜２８５は燃焼室５の燃焼質量割合３１％時における容積Ｖｍｂ３１ｐを算出する部分である。まずステップ２８２では燃焼質量割合が３１％のときのクランク角の前回値θｍｂ３１ｐｚ［ｄｅｇＢＴＤＣ］を次式により算出する。

θｍｂ３１ｐｚ＝θｉｇｎｚ−ＩＧＮＤＥＡＤｒｅｆ１ｚ−ＢＵＲＮ１ｒｅｆ１ｚ
−ＢＵＲＮＮ２ｒｅｆ１ｚ／２
…（６８）
ステップ２８３では燃焼質量割合が２％のときのクランク角の前回値と、燃焼質量割合が２％のときのクランク角の今回値との差θ２ｐｄｉｆ［ｄｅｇ］を次式により算出する。

θ２ｐｄｉｆ＝（θｉｇｎｚ−ＩＧＮＤＥＡＤｒｅｆ１ｚ−ＢＵＲＮ１ｒｅｆ１ｚ）
−（θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１−ＢＵＲＮ１ｒｅｆ１） …（６９）
ステップ２８４では燃焼質量割合が３１％のときのクランク角の今回値θｍｂ３１ｐ［ｄｅｇＢＴＤＣ］を次式により算出する。

θｍｂ３１ｐ＝θｍｂ３１ｐｚ−θ２ｐｄｉｆ …（７０）
ステップ２８５ではこのクランク角θｍｂ３１ｐから燃焼室５の燃焼質量割合３１％時における容積Ｖｍｂ３１ｐ［ｍ³］を次式により算出する。この操作は図３４のステップ２４２の操作と同様である。

Ｖｍｂ３１ｐ＝ｆ２５（θｍｂ３１ｐ） …（７１）
ステップ２８６ではこの燃焼室５の燃焼質量割合３１％時における容積Ｖｍｂ３１ｐから燃焼質量割合３１％時の有効圧縮比Ｅｍｂ３１ｐを次式により算出する。

Ｅｍｂ３１ｐ＝Ｖｍｂ３１ｐ／ＶＩＶＣ …（７２）
ステップ２８７〜２９１によりＴＣとＰＣを算出するための操作は図３４のステップ２４４〜２４８によりＴＣとＰＣを算出するための操作と同様である。異なるのは図３４のステップ２４７、２４８においては燃焼質量割合１％時の圧縮比であるＥｍｂ１ｐと、燃焼質量割合１％時の燃焼質量割合であるＢＲ１ｐとを用いたのに対して、図３６のステップ２９０、２９１では燃焼質量割合３１％時の圧縮比であるＥｍｂ３１ｐと、燃焼質量割合３１％時の燃焼質量割合であるＢＲ３１ｐとを用いる点である。すなわち、図３６のステップ２９０、２９１では燃焼室５内の燃料が燃焼したときの平均温度ＴＣ［Ｋ］と平均圧力ＰＣ［Ｐａ］を次式により算出する。

ＴＣ＝ＴＣ０×Ｅｍｂ３１ｐ＾０．３５
＋ＣＦ＃×ＱＩＮＪ×ＢＲ３１ｐ／（ＭＡＳＳＺ＋ＷＩＤＲＹ＋ＱＩＮＪ）
…（７３）
ＰＣ＝ＰＣ０×Ｅｍｂ３１ｐ＾１．３５×ＴＣ／ＴＣ０／Ｅｍｂ３１ｐ＾０．３５
…（７４）
ただし、ＣＦ＃：燃料の低位発熱量、
ステップ２９２〜２９５により第１主燃焼期間を算出するための操作は図１２のステップ１８８〜１９１により主燃焼期間を算出するための操作と同様である。異なるのは図１２のステップ１８８において燃焼室５の圧縮上死点時期における温度ＴＴＤＣと圧力ＰＴＤＣとを用いたのに対して、図３６のステップ２９２では上記（７３）式、（７４）式で求めた温度ＴＣと圧力ＰＣとを用いる点、また図１２のステップ１９１において主燃焼期間ＢＵＲＮ２の終了時期における燃焼質量割合であるＢＲ２（＝５８％）と、燃焼室５の圧縮上死点時期における容積であるＶＴＤＣとを用いたのに対して、図３６のステップ２９５では主燃焼期間の中間時期の近傍における燃焼質量割合である３１％（＝ＢＲ３１ｐ）と、燃焼室５の燃焼質量割合３１％時における容積であるＶｍｂ３１ｐとを用いる点である。すなわち、ステップ２９２では第１主燃焼期間における層流燃焼速度ＳＬ２［ｍ／ｓｅｃ］を、またステップ２９５では第１主燃焼期間ＢＵＲＮ２ｒｅｆ１［ｄｅｇ］をそれぞれ次式により算出する。

ＳＬ２＝ＳＬstd×（ＴＣ／Ｔstd）^2.18×（ＰＣ／Ｐstd）^-0.16 …（７５）
ただし、Ｔstd ：基準温度［Ｋ］、
Ｐstd ：基準圧力［Ｐａ］、
ＳＬstd：基準温度Ｔstdと基準圧力Ｐstdにおける基準層流燃焼
速度［ｍ／ｓｅｃ］、
ＢＵＲＮ２ｒｅｆ１＝｛（ＮＲＰＭ×６）×ＢＲ３１ｐ×Ｖｍｂ３１ｐ｝
／（ＲＰＲＯＢＡ×ＡＦ２×ＦＬＡＭＥ１） …（７６）
ただし、ＡＦ２ ：火炎核の反応面積（固定値）［ｍ²］、
このようにして、第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１、第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１、第１主燃焼期間ＢＵＲＮ２ｒｅｆ１の算出を終了したら図３３（Ａ）のステップ２２７に戻り燃焼質量割合の第１近似式Ｘ１’を算出する。この第１近似式の算出については図３７のフローにより説明する。

ここで、図３７のフローを説明する前に、どのようにして燃焼質量割合の近似式を求めたのかを先に説明する。

燃焼質量割合Ｘ［％］は一般にクランク角θ［ｄｅｇ］を変数とする指数関数として次式により記述されることが多い。

Ｘ＝１−ｅｘｐ［−ａ・｛（θ−θｓ）／θＢ｝＾（ｍ＋１）］
…（補４）
ただし、θｓ：熱発生開始クランク角、
θＢ：燃焼期間相当クランク角、
ａ、ｍ：プロファイルを決める定数、
（補４）式の定数ａ、ｍは燃焼室形状や点火プラグ位置、及び燃焼室内のガス流動特性等、エンジン機種によって一義的に定まる値である。従って、定数ａ、ｍは実験等により予め決めておくことができる。（補４）式はＷｉｅｂｅ関数といわれる。

本実施形態は、所定のパラメータに基づいてこれらの定数ａ、ｍを算出し、これら算出した定数ａ、ｍを上記（補４）式に代入した式を燃焼質量割合の近似式として求める。具体的に説明すると、（補４）式を次のように変形する。

Ｘ＝１−ｅｘｐ［−｛ａ／θＢ＾（ｍ＋１）｝・（θ−θｓ）＾（ｍ＋１）］
…（補５）
この（補５）式において定数ｋを、
ｋ≡ａ／｛θＢ＾（ｍ＋１）｝ …（補６）
の式で定義すると、（補４）式に代えて次式が得られる。

Ｘ＝１−ｅｘｐ［−ｋ・（θ−θｓ）＾（ｍ＋１）］…（補７）
さて、（補４）式のＷｉｅｂｅ関数を本実施形態に適用することを考える。

まず、熱発生開始クランク角θｓ［ｄｅｇＢＴＤＣ］は本実施形態では次式で与えられる。

θｓ＝θｉｇｎ［ｄｅｇＢＴＤＣ］−ＩＧＮＤＥＡＤ［ｄｅｇ］
…（補８）
いま、燃焼質量割合が２％のときのクランク角をθ2、燃焼質量割合が６０％のときのクランク角をθ60とすれば、これらを上記（補４）式に代入して次の２つの式が成立する。

０．０２＝１−［−ａ・｛（θ2−θｓ）／θＢ｝＾（ｍ＋１）］
…（補９）
０．６＝１−［−ａ・｛（θ60−θｓ）／θＢ｝＾（ｍ＋１）］
…（補１０）
（補９）、（補１０）式をそれぞれ変形して、
−ａ・｛（θ2−θｓ）／θＢ｝＾（ｍ＋１）＝ｌｎ（１−０．０２）
＝ｌｎ（０．９８）
…（補１１）
−ａ・｛（θ60−θｓ）／θＢ｝＾（ｍ＋１）＝ｌｎ（１−０．６）
＝ｌｎ（０．４）
…（補１２）
（補１１）、（補１２）式の両辺をそれぞれ割り算して、
｛（θ2−θｓ）／（θ60−θｓ）｝＾（ｍ＋１）
＝ｌｎ（０．９８）／ｌｎ（０．４）
…（補１３）
（補１３）式を定数ｍについて整理すると次式が得られる。

ｍ＝ｌｎ｛ｌｎ（０．９８）／ｌｎ（０．４）｝
／ｌｎ｛（θ2−θｓ）／（θ60−θｓ）｝−１
…（補１４）
ここで、（補１４）式の「ｌｎ」は自然対数を表している。

また、上記（補１１）式を変形すると次式が得られる。

ａ／｛θＢ＾（ｍ＋１）｝＝−ｌｎ（０．９８）／｛（θ2−θｓ）＾（ｍ＋１）｝
…（補１５）
（補１５）、（補６）式より、上記の定数ｋはこの定数ｍを用いて次式により表される。

ｋ＝ａ／｛θＢ＾（ｍ＋１）｝
＝−ｌｎ（０．９８）／｛（θ2−θｓ）＾（ｍ＋１）｝
…（補１６）
以上をまとめると、上記（補８）式により得られる熱発生開始クランク角θｓを上記（補１４）式に代入すれば、定数ｍが求まり、その求めた定数ｍを（補１６）式に代入すれば定数ｋが求まる。このようにして、定数ｍと定数ｋが求まると、この求めた定数ｍ、ｋを上記（補７）式に代入した式を燃焼質量割合の近似式として定めることができる。

図３７（図３３（Ａ）のステップ２２７のサブルーチン）は燃焼質量割合の第１近似式Ｘ１’を算出するためのものである。

図３７においてステップ３０１では、図３３（Ａ）のステップ２２２で得ている点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］、図３３（Ａ）のステップ２２５で算出している第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１［ｄｅｇ］を読み込む。

ステップ３０２では次式により熱発生開始クランク角θｓ［ｄｅｇＢＴＤＣ］を算出する。

θｓ＝θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１ …（７７）
（７７）式は上記（補８）式と同じである。

ステップ３０３ではこの熱発生開始クランク角θｓを用いて次式により定数ｍを算出する。

ｍ＝ｌｎ｛ｌｎ（０．９８）／ｌｎ（０．４）｝
／ｌｎ｛（θ2−θｓ）／（θ60−θｓ）｝−１ …（７８）
ただし、θ2 ：燃焼質量割合が２％のときのクランク角、
θ60：燃焼質量割合が６０％のときのクランク角、
（７８）式は上記（補１４）式と同じである。ここで、（７８）式のθ2、θ60は基本点火時期ＭＢＴＣＡＬより求める。または適合により定めておく。

ステップ３０４ではこのようにして算出した定数ｍを用いて次式によりもう一つの定数ｋを算出する。

ｋ＝ａ／｛θＢ＾（ｍ＋１）｝
＝−ｌｎ（０．９８）／｛（θ2−θｓ）＾（ｍ＋１）｝ …（７９）
（７９）式は上記（補１６）式と同じである。

ステップ３０５ではこの定数ｋと上記の定数ｍとを用いて次式により燃焼質量割合の第１近似式Ｘ１’［％］を次式により算出する。

Ｘ１’＝１−ｅｘｐ［−ｋ・（θ−θｓ）＾（ｍ＋１）］ …（８０）
（８０）式は上記（補７）式と同じである。

このようにして燃焼質量割合の第１近似式Ｘ１’の算出を終了したら図３３（Ａ）のステップ２２８に戻り、第１図示正仕事推定値Ｗｐｏｓ１を算出する。この第１図示正仕事推定値Ｗｐｏｓ１の算出については図３８、図３９のフローにより説明する。

図３８、図３９（図３３（Ａ）のステップ２２８のサブルーチン）は第１正仕事推定値Ｗｐｏｓ１を算出するためのものである。

図３８においてステップ３１１では、図３３（Ａ）のステップ２２２で得ている点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］、図３３（Ａ）のステップ２２５で算出している第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１［ｄｅｇ］、燃料量ＱＩＮＪ［ｋｇ］、 図５のステップ１２で算出している燃焼室５の吸気弁閉時期における容積ＶＩＶＣ［ｍ³］、図５のステップ１３で算出している燃焼室５の吸気弁閉時期における温度ＴＩＮＩ［Ｋ］、図５のステップ１４で算出している燃焼室５の吸気弁閉時期における圧力ＰＩＮＩ［Ｐａ］、温度センサ４８により検出される燃焼室壁温Ｔｗａｌｌ［Ｋ］、目標当量比ＴＦＢＹＡを読み込む。

ステップ３１２では吸気弁閉時期ＩＶＣ［ｄｅｇＢＴＤＣ］をクランク角前回値θ(old)［ｄｅｇＢＴＤＣ］にセットする。

ステップ３１３ではクランク角前回値θ(old)を所定角度ｃｏｎｓｔ２［ｄｅｇ］だけ進めた値をクランク角今回値θ(new)としてセットする。ここで、所定角度ｃｏｎｓｔ２は後述する瞬間仕事ｄＷを算出する周期を定める値であり、ｃｏｎｓｔ２の値を大きくし過ぎると演算時間は短くて済むものの図示正仕事推定値の演算精度が低下し、この逆にｃｏｎｓｔ２の値を小さくし過ぎると演算精度は上がるものの演算時間が長引くことになるので、最適な値を適合により設定する。ここでは所定角度ｃｏｎｓｔ２は１ｄｅｇ程度で考えている。

ステップ３１４ではこのクランク角今回値θ(new)と排気弁開時期ＥＶＯとを比較する。クランク角今回値θ(new)が排気弁開時期ＥＶＯに達していなければステップ３１５に進み燃焼室５のクランク角今回値θ(new)における容積Ｖｃｙｌ(new)［ｍ³］を算出する。この算出にはクランク角今回値θ(new)から次式によりまずピストン変位［ｍ］を算出する。

ピストン変位＝（クランク長＋コンロッド長）−クランク長・cos（θ(new)）
−ｓｑｒｔ｛コンロッド長＾２−クランク長＾２・sin（（θ(new)＾２）｝
…（８１）
（８１）式のクランク長［ｍ］、コンロッド長［ｍ］はエンジンの仕様により定まる値である。（８１）式の「ｓｑｒｔ」は平方根を採ることを表している。ただし、ここには簡単のためクランクシャフトの回転中心がシリンダの中心軸からオフセットしていない場合の式を示している。

上記（８１）式により算出されるピストン変位を用いて次式により燃焼室５のクランク角今回値θ(new)における容積Ｖｃｙｌ(new)［ｍ³］を算出する。

Ｖｃｙｌ(new)＝（π／４）・Ｄ＾２・ピストン変位＋Ｖｃ …（８２）
ただし、Ｄ ：シリンダボア径［ｍ］、
Ｖｃ：隙間容積［ｍ³］、
ステップ３１６では燃焼室５のクランク角今回値θ(new)における圧力Ｐ(new)と温度Ｔ(new)を次式により算出する。

Ｐ(new)＝Ｐ(old)×（Ｖｃｙｌ(old)／Ｖｃｙｌ(new)）＾κ …（８３）
Ｔ(new)＝Ｔ(old)×（Ｖｃｙｌ(old)／Ｖｃｙｌ(new)）＾（κ−１）
…（８４）
（８３）、（８４）式は燃焼室５内でガスが断熱圧縮されるとしたときの圧力と温度を熱力学の法則に則って表した式である。

ここで、（８３）、（８４）式のＶｃｙｌ(old)は燃焼室５のクランク角前回値θ(old)における容積前回値、またＰ(old)とＴ(old)は燃焼室５のクランク角前回値θ(old)における圧力前回値と温度前回値で、後述する図３９のステップ３２９により得られている。Ｐ(old)の初期値としては図５のステップ１４で算出している燃焼室５の吸気弁閉時期における圧力ＰＩＮＩを、Ｔ(old)の初期値としては図５のステップ１３で算出している燃焼室５の吸気弁閉時期における温度ＴＩＮＩを、Ｖｃｙｌ(old)の初期値としては図５のステップ１２で算出している燃焼室５の吸気弁閉時期における容積ＶＩＶＣを入れておく。

ピストン６が上昇して燃焼室５内でガスが断熱圧縮されるときにはＶｃｙｌ(old)／Ｖｃｙｌ(new)＞１．０となり、従って、上記（８３）、（８４）式よりクランク角前回値における圧力前回値Ｐ(old)よりもクランク角今回値における圧力Ｐ(new)のほうが高くなり、かつクランク角前回値における温度前回値Ｔ(old)よりもクランク角今回値における温度Ｔ(new)のほうが高くなる。

ステップ３１７では図３３（Ａ）のステップ２２７で求めた燃焼質量割合の第１近似式Ｘ１’にクランク角今回値θ(new)を代入することによって、つまりクランク角今回値θ(new)での燃焼質量割合ＭＢ(new)を次式により算出する。

ＭＢ(new)＝１−ｅｘｐ［−ｋ・（θ(new)
−（θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１））＾（ｍ＋１）］
…（８５）
ステップ３１８ではこの燃焼質量割合ＭＢ(new)と所定値ｃｏｎｓｔ３を、またクランク角今回値θ(new)と熱発生開始クランク角（θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１）とを比較する。ここで、燃焼が終了するのは燃焼質量割合が１００％になるときであるので、所定値ｃｏｎｓｔ３としてはこの１００％に十分近い値を設定する。つまり、燃焼質量割合ＭＢ(new)が所定値ｃｏｎｓｔ３以上になると燃焼が終了していることを表す。

クランク角今回値θ(new)が熱発生開始クランク角（θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１）に達していないときにはまだ熱発生がないし、燃焼質量割合ＭＢ(new)が所定値ｃｏｎｓｔ３以上であるときには燃焼が終了しているので、これらの場合にはステップ３２０に進んでクランク角前回値からクランク角今回値までの発生熱量ｄＱｂｕｒｎ＝０とする。

一方、クランク角今回値θ(new)が熱発生開始クランク角（θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１）を過ぎておりかつ燃焼質量割合ＭＢ(new)が所定値ｃｏｎｓｔ３未満であるときには熱発生がありかつ燃焼が終了していないと判断しステップ３１８よりステップ３１９に進んでクランク角前回値からクランク角今回値までの発生熱量ｄＱｂｕｒｎ［ｋＪ］を次式により算出する。

ｄＱｂｕｒｎ＝（ＭＢ(new)−ＭＢ(old)）・ＱＩＮＪ・ＣＦ＃…（８６）
ただし、ＣＦ＃：燃料の低位発熱量［ｋＪ／ｋｇ］、
ここで、（８６）式のＭＢ(old)はクランク角前回値θ(old)における燃焼質量割合前回値で、後述する図３９のステップ３２９により得られている。

ステップ３２１ではクランク角前回値からクランク角今回値までの冷却損失ｄＱｌｏｓｓ［ｋＪ］を次式により算出する。この算出については図４０のフローにより説明する。

図４０（図３８のステップ３２１のサブルーチン）においてステップ３４１では熱伝達係数ｈ［ｋｃａｌ／（ｍ²・ｈ・Ｋ）］を次式により算出する。

ｈ＝１１０・Ｄ＾（−０．２）・Ｐ(old)＾０．８・Ｔ(old)＾（−０．５３）
・［Ｃ１・Ｃｍ＋Ｃ２・（Ｖｃｙｌ(old)・ＴＩＮＩ）
／（ＰＩＮＩ・ＶＩＶＣ）・（Ｐ(old)−ＰＭ(old)）］＾０．８
…（８７）
ただし、Ｄ ：シリンダボア径［ｍ］、
Ｃｍ ：平均ピストン速度［ｍ／ｓ］、
ＰＭ(old)：モータリング圧力の前回値［Ｐａ］、
ここで、（８７）式のＰＭ(old)は燃焼室５のクランク角前回値θ(old)におけるモータリング時の圧力前回値で、後述するステップ３４４により得られている。ＰＭ(old)の初期値としては圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬを入れておく。なお、上記（８７）式のＰ(old)、Ｔ(old)、Ｖｃｙｌ(old)はそれぞれ燃焼室５のクランク角前回値θ(old)における圧力前回値、温度前回値、容積前回値で、後述する図３９のステップ３２９により得られている。

上記（８７）式の２つの定数Ｃ１、Ｃ２は一般に用いられている定数で、例えば一方の定数Ｃ１は排気行程において６．１８、圧縮行程と膨張行程において２．２８である。また、他方の定数Ｃ２は３．２４×１０＾（−３）［ｍ／（ｓ・Ｋ）］である。

ステップ３４２では、この熱伝達係数ｈを用いて次式によりクランク角前回値からクランク角今回値までの冷却損失ｄＱｌｏｓｓ［ｋＪ］を次式により算出する。

ｄＱｌｏｓｓ＝ｈ・Ａ・（Ｔ(old)−Ｔｗａｌｌ）・Δｔ …（８８）
ただし、Ａ ：燃焼室の表面積［ｍ²］、
Δｔ：クランク角前回値からクランク角今回値までの演算間隔［ｈ］
（８８）式右辺の温度前回値Ｔ(old)は燃焼室壁温Ｔｗａｌｌより高いので、（８８）式左辺の冷却損失ｄＱｌｏｓｓは正の値である。

（８８）式のΔｔは次のようにして求めればよい。すなわち、上記の所定角度ｃｏｎｓｔ２［ｄｅｇ］がクランク角前回値からクランク角今回値までの演算間隔であるので、これをそのときのエンジン回転速度ＮＲＰＭを用いて時間［ｈ］に換算する。

ステップ３４３では、燃焼室５のクランク角今回値θ(new)におけるモータリング時の圧力ＰＭ(new)を上記（８１）式と同様の式である次式により算出する。

ＰＭ(new)＝ＰＭ(old)×（Ｖｃｙｌ(old)／Ｖｃｙｌ(new)）＾κ…（８９）
ステップ３４４ではこのようにして算出した燃焼室５のクランク角今回値θ(new)におけるモータリング時の圧力ＰＭ(new)を、次回に備えて燃焼室５のクランク角前回値θ(new)におけるモータリング時の圧力前回値ＰＭ(old)に移しておく。

このようにしてクランク角前回値からクランク角今回値までの冷却損失ｄＱｌｏｓｓの算出を終了したら図３８のステップ３２２に戻り、定容比熱Ｃｖ［Ｊ／（ｋｇ・Ｋ）］を算出する。これは燃焼室５のクランク角前回値における温度前回値Ｔ(old)からそれぞれ所定のテーブルを検索して既燃ガスの定圧比熱Ｃｐ１及び比熱比κ１を、また同じく燃焼室５のクランク角前回値における温度前回値Ｔ(old)からそれぞれ所定のテーブルを検索して未燃ガスの定圧比熱Ｃｐ２及び比熱比κ２を求める。

ここで、定圧比熱Ｃｐ１、Ｃｐ２を対応する比熱比でκ１、κ２で除算すれば、各定容比熱Ｃｖ１（＝Ｃｐ１／κ１）、Ｃｖ２（＝Ｃｐ２／κ２）を求めることができるので、この求めた既燃ガスの定容比熱Ｃｖ１と未燃ガスの定容比熱Ｃｖ２との２つを燃焼質量割合前回値ＭＢ(old)で加重平均し、その算出した加重平均値を燃焼室５内のガスの定容比熱Ｃｖとして設定する。すなわち、次式により定容比熱Ｃｖを算出する。

Ｃｖ＝（Ｃｐ１／κ１）×ＭＢ(old)＋（Ｃｐ２／κ２）×（１−ＭＢ(old)）
…（９０）
ただし、Ｃｐ１：既燃ガスの定圧比熱、
Ｃｐ２：未燃ガスの定圧比熱、
κ１ ：既燃ガスの比熱比、
κ２ ：未燃ガスの比熱比、
ステップ３２３では燃焼室５内の全ガス量ＧＡＳＡＬＬ［ｋｇ］を算出する。これは燃料量ＱＩＮＪ、目標当量比ＴＦＢＹＡから次式により算出すればよい。

ＧＡＳＡＬＬ＝ＱＩＮＪ×（１＋１／ＴＦＢＹＡ） …（９１）
図３９のステップ３２４では、燃焼室５のクランク角今回値θ(new)における温度Ｔ(new)及び圧力Ｐ(new)を次式により改めて算出する。

Ｔ(new)＝Ｔ(old)×（Ｖｃｙｌ(old)／Ｖｃｙｌ(new)）＾（κ−１）
＋（ｄＱｂｕｒｎ−ｄＱｌｏｓｓ）／（ＧＡＳＡＬＬ・Ｃｖ）
…（９２）
Ｐ(new)＝｛Ｐ(old)×（Ｖｃｙｌ(old))／Ｖｃｙｌ(new）＾κ｝
・Ｔ(new)／｛Ｔ(old)×（Ｖｃｙｌ(old)／Ｖｃｙｌ(new)）＾（κ−１）｝
…（９３）
ここで、（９２）式右辺の第１項は燃焼室５内でガスが断熱圧縮されるとしたときの温度上昇分であるため上記の（８４）式と同じであり、（９２）式右辺の第２項が燃焼による温度上昇分である。また、（９３）式右辺の｛Ｐ(old)×（Ｖｃｙｌ(old))／Ｖｃｙｌ(new）＾κ｝は燃焼室５内でガスが断熱圧縮されるとしたときの圧力上昇分であるため上記の（８３）式と同じであり、（９３）式右辺のＴ(new)／｛Ｔ(old)×（Ｖｃｙｌ(old)／Ｖｃｙｌ(new)）＾（κ−１）｝が燃焼による温度上昇に伴って圧力が上昇する程度である。

ステップ３２５ではクランク角前回値からクランク角今回値までの平均仕事を瞬間仕事ｄＷ［ｋＪ］として次式により算出する。

ｄＷ＝（Ｐ(new)＋Ｐ(old)）×（Ｖｃｙｌ(new)−Ｖｃｙｌ(old)）／２
…（９４）
ここで、ピストン６が圧縮上死点を過ぎて下死点ＴＤＣに向かって下降するときには容積Ｖｃｙｌ(new)が容積前回値Ｖｃｙｌ(old)より大きくなり、（９４）式左辺の値は正の値となるのに対して、圧縮上死点ＴＤＣに向かって上昇するときには容積Ｖｃｙｌ(new)が容積前回値Ｖｃｙｌ(old)より小さく（９４）式左辺の値が負の値となるので、（９４）式左辺の値は正負の値を採り得るが、ここでは、絶対値を採って瞬間仕事ｄＷが常に正の値で得られるようにしておく。

ステップ３２６ではクランク角今回値θ(new)と圧縮上死点ＴＤＣとを比較する。クランク角今回値θ(new)が圧縮上死点ＴＤＣに達していなければエンジンが行う仕事は負の仕事であると判断し、ステップ３２７に進んでクランク角今回値θ(new)での負の仕事Ｗｎ(new)を次式により算出する。

Ｗｎ(new)＝Ｗｎ(old)＋ｄＷ …（９５）
ここで、（９５）式右辺の瞬間仕事ｄＷは常に正の値で算出しているので、（９５）式左辺の負の仕事Ｗｎ(new)も正の値で算出されている。

一方、クランク角今回値θ(new)が圧縮上死点ＴＤＣに達した後はエンジンが行う仕事は正の仕事であると判断してステップ３２６よりステップ３２８に進みクランク角今回値θ(new)での正の仕事Ｗｐ(new)を次式により算出する。

Ｗｐ(new)＝Ｗｐ(old)＋ｄＷ …（９６）
上記（９５）、（９６）式右辺のＷｎ(old)、Ｗｐ(old)はそれぞれクランク角前回値θ(old)での負の仕事前回値、クランク角前回値θ(old)での正の仕事前回値で、後述するステップ３２９により得られている。

ステップ３２９ではクランク角今回値θ(new)、クランク角今回値θ(new)での燃焼質量割合ＭＢ(new)、燃焼室５のクランク角今回値θ(new)における圧力Ｐ(new)、燃焼室５のクランク角今回値θ(new)における温度Ｔ(new)、燃焼室５のクランク角今回値θ(new)における容積Ｖｃｙｌ(new)、クランク角今回値θ(new)での負の仕事Ｗｎ(new)、クランク角今回値θ(new)での正の仕事ｗｐ(new)を、次回に備えてそれぞれクランク角前回値θ(old)、クランク角前回値θ(old)での燃焼質量割合前回値ＭＢ(old)、燃焼室５のクランク角前回値θ(old)における圧力前回値Ｐ(old)、燃焼室５のクランク角前回値θ(old)における温度前回値Ｔ(old)、燃焼室５のクランク角前回値θ(old)における容積前回値Ｖｃｙｌ(old)、クランク角前回値θ(old)での負の仕事前回値Ｗｎ(old)、クランク角前回値θ(old)での正の仕事前回値Ｗｐ(old)に移しておく。

ここで、クランク角前回値θ(old)の初期値は図３８のステップ３１２で与えられている。また、図３８のステップ３１６で前述したように、Ｐ(old)の初期値としては図５のステップ１４で算出している燃焼室５の吸気弁閉時期における圧力ＰＩＮＩを、Ｔ(old)の初期値としては図５のステップ１３で算出している燃焼室５の吸気弁閉時期における温度ＴＩＮＩを、Ｖｃｙｌ(old)の初期値としては図５のステップ１２で算出している燃焼室５の吸気弁閉時期における容積ＶＩＶＣを入れておく。残った前回値であるＭＢ(old)、Ｗｎ(old)、Ｗｐ(old)の各初期値としては図３８のフローの処理開始時にゼロを入れておく。

ステップ３３０では、図３８のステップ３１３と同じに、クランク角前回値θ(old)を所定角度ｃｏｎｓｔ２だけ進めた値をクランク角今回値θ(new)としてセットした後に図３８のステップ３１４に戻り、図３８のステップ３１４〜３２３、図３９のステップ３２４〜３３０の操作を繰り返す。つまり、クランク角θが熱発生開始クランク角（θｉｇｎ−ＩＧＮＤＥＡＤｒｅｆ１）になると、発生熱量ｄＱｂｕｒｎと冷却損失ｄＱｌｏｓｓとが所定クランク角ｃｏｎｓｔ２毎に算出され、これら算出された発生熱量ｄＱｂｕｒｎ及び冷却損失ｄＱｌｏｓｓ並びに熱力学の法則に基づいて燃焼室５内の燃焼中の圧力Ｐと温度Ｔが、さらにはこのうち圧力Ｐに基づいて瞬時仕事ｄＷが所定クランク角ｃｏｎｓｔ２毎に算出される。この場合に、圧縮上死点ＴＤＣまではその瞬時仕事ｄＷは負の仕事であるとしてＷｎに、また圧縮上死点ＴＤＣ後になるとその瞬時仕事ｄＷは正の仕事であるとしてＷｐにそれぞれ所定クランク角ｃｏｎｓｔ２毎に積算されてゆく。

やがて、図３８のステップ３１４でクランク角今回値θ(new)が排気弁開時期ＥＶＯに達したときには燃焼室５内の燃焼ガスによる仕事の終了タイミングであると判断し、ステップ３１４よりステップ３３１に進み次式により第１図示正仕事推定値Ｗｐｏｓ１［ｋＪ］を算出し、これで図３８、図３９での処理を終了する。

Ｗｐｏｓ１＝Ｗｐ(old)−Ｗｎ(old) …（９７）
このようにして第１図示正仕事推定値Ｗｐｏｓ１の算出を終了したら図３３（Ａ）のステップ２２９に戻り、第１軸トルク推定値Ｔｒｑｅｓｔ１を算出する。この算出については図４１のフローにより説明する。

図４１（図３３（Ａ）のステップ２２９のサブルーチン）においてステップ３５１では、図３３（Ａ）のステップ２２８により算出されている第１図示正仕事推定値Ｗｐｏｓ１［ｋＪ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、エンジンの負荷、回転速度ＮＲＰＭ［ｒｐｍ］を読み込む。

ステプ３５２ではポンピングロスＷｎｅｇ１［ｋＪ］を次式により算出する。

Ｗｎｅｇ＝ＰＣＯＬ×Ｈx …（９８）
ただし、Ｈx：ピストンピンのＴＤＣからの距離の最大値と最小値の差［ｍ］
ステップ３５３では摩擦損失Ｗｎｅｇ２［ｋＪ］を算出する。これはエンジンの負荷と回転速度ＮＲＰＭから図４２を内容とするマップを検索することにより求めればよい。図４２のように摩擦損失Ｗｎｅｇ２は高負荷になるほど、また回転速度ＮＲＰＭが高くなるほど大きくなる値である。

ステップ３５４では第１図示正仕事推定値Ｗｐｏｓ１にポンピングロスＷｎｅｇ１と摩擦損失Ｗｎｅｇ２を加算して、つまり次式により第１軸トルク推定値Ｔｒｑｅｓｔ１［ｋＪ］を算出する。

Ｔｒｑｅｓｔ１＝Ｗｐｏｓ＋Ｗｎｅｇ１＋Ｗｎｅｇ２ …（９９）
これで図３３（Ａ）のフローの説明を全て終える。

次に、図３３（Ｂ）（図２のステップ２のサブルーチン）において上記の図３３（Ａ）と異なるのは点火時期の扱いだけである。すなわち、図３３（Ｂ）においてステップ２３１では図１３のステップ４３で算出されている基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］、点火時期リタード量前回値ＲＥＴ(old)［ｄｅｇ］を読み込む。ここで、点火時期リタード量前回値ＲＥＴ(old)は図２のステップ６で得られる値である（後述する）。

ステップ２３２では点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］の値を点火時期前回値θｉｇｎｚ［ｄｅｇＢＴＤＣ］に移すと共に、基本点火時期ＭＢＴＣＡＬから点火時期リタード量前回値ＲＥＴ(old)を差し引いた値（つまり基本点火時期ＭＢＴＣＡＬより点火時期リタード量前回値ＲＥＴ(old)だけリタード側の値）を点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］として算出する。なお、ＲＥＴ(old)の初期値としてはゼロを入れておく。

ステップ２３４では第２初期燃焼期間ＢＵＲＮ１ｒｅｆ２［ｄｅｇ］を、ステップ２３５では第２点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ２［ｄｅｇ］を、ステップ２３６では第２主燃焼期間ＢＵＲＮ２ｒｅｆ２［ｄｅｇ］を、ステップ２３７では燃焼質量割合の第２近似式Ｘ２’を、ステップ２３８では第２図示正仕事推定値Ｗｐｏｓ２［ｋＪ］を、ステップ２３９では第２軸トルク推定値Ｔｒｑｅｓｔ２［ｋＪ］をそれぞれ算出する。

ここで、第２初期燃焼期間ＢＵＲＮ１ｒｅｆ２の算出については図４３のフローにより、第２点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ２の算出については図４４のフローにより、第２主燃焼期間ＢＵＲＮ２ｒｅｆ２の算出については図４５のフローにより、燃焼質量割合の第２近似式Ｘ２’の算出については図４６のフローにより、第２図示正仕事推定値Ｗｐｏｓ２の算出については図４７、図４８のフローにより、第２軸トルク推定値Ｔｒｑｅｓｔ２の算出については図４９（図３３（Ｂ）のステップ２３９のサブルーチン）のフローにより説明する。

この場合に、図４３（図３３（Ｂ）のステップ２３４のサブルーチン）、図４４（図３３（Ｂ）、図４５（図３３（Ｂ）のステップ２３６のサブルーチン）、図４６（図３３（Ｂ）のステップ２３７のサブルーチン）、図４７、図４８（図３３（Ｂ）のステップ２３８のサブルーチン）において、上記の図３４、図３５、図３６、図３７、図３８、図３９と異なるのは、図３４のステップ２４１における点火時期前回値θｉｇｎｚ、図３５のステップ２６１における点火時期θｉｇｎ、図３６のステップ２８１における点火時期θｉｇｎ、点火時期前回値θｉｇｎｚ、図３７のステップ３０１における点火時期θｉｇｎ、図３８のステップ３１１における点火時期θｉｇｎがいずれも基本点火時期ＭＢＴＣＡＬであるのに対して、図４３のステップ３６１における点火時期前回値θｉｇｎｚ、図４４のステップ３８１における点火時期θｉｇｎ、図４５のステップ４０１における点火時期θｉｇｎ、点火時期前回値θｉｇｎｚ、図４６のステップ４２１における点火時期θｉｇｎ、図４７のステップ４３１における点火時期θｉｇｎがいずれも基本点火時期ＭＢＴＣＡＬより点火時期リタード量前回値ＲＥＴ(old)だけリタードさせた値である点にある。他は同様であるので、個々の説明は省略する。

このように、図３３（Ｂ）において初期燃焼期間、点火無駄時間相当クランク角、主燃焼期間、燃焼質量割合の近似式、図示正仕事推定値、軸トルク推定値を、図３３（Ａ）とは別に算出するのは、基本点火時期ＭＢＴＣＡＬより点火時期リタード量前回値ＲＥＴ(old)だけリタードさせた点火時期で点火したときの初期燃焼期間、点火無駄時間相当クランク角、主燃焼期間、燃焼質量割合の近似式、図示正仕事推定値、軸トルク推定値が、基本点火時期ＭＢＴＣＡＬで点火したときの初期燃焼期間、点火無駄時間相当クランク角、主燃焼期間、燃焼質量割合の近似式、図示正仕事推定値、軸トルク推定値とは異なる値になるためである。

このようにして第１軸トルク推定値Ｔｒｑｅｓｔ１、第２軸トルク推定値Ｔｒｑｅｓｔ２の算出を共に終了したら、図２に戻りステップ３で要求軸トルクＴｒｑｔｇｔを算出する。この要求軸トルクＴｒｑｔｇｔを算出するには、アクセルセンサ４２により検出されるアクセル開度（アクセルペダルの踏み込み量）ＡＰＯと回転速度ＮＲＰＭとから図５０を内容とするマップを検索することにより求めればよい。図５０に示したように要求軸トルクＴｒｑｔｇｔはアクセル開度ＡＰＯが大きくなるほどまた回転速度ＮＲＰＭが高くなるほど大きくなる値である。

ステップ４では点火時期リタード量ＲＥＴ(new)を算出する。この点火時期リタード量ＲＥＴ(new)の算出については図５１のフローにより説明する。

図５１（図２のステップ４のサブルーチン）においてステップ４７１では図２のステップ２により算出されている第１軸トルク推定値Ｔｒｑｅｓｔ１［ｋＪ］、図２のステップ３により算出されている要求軸トルクＴｒｑｔｇｔ［ｋＪ］を読み込む。

ステップ４７２では第１軸トルク推定値Ｔｒｑｅｓｔ１と要求軸トルクＴｒｑｔｇｔとの差に基づいて、つまり次式により点火時期リタード量ＲＥＴ(new)［ｄｅｇ］を算出する。

ＲＥＴ(new)＝（Ｔｒｑｅｓｔ１−Ｔｒｑｔｇｔ）×ｃｏｎｓｔ２０１
…（１００）
ただし、ｃｏｎｓｔ２０１：所定値、
この点火時期リタード量ＲＥＴ(new)は、基本点火時期ＭＢＴＣＡＬでの軸トルク推定値である第１軸トルク推定値Ｔｒｑｅｓｔ１が要求軸トルクＴｒｑｔｇｔより大きい、つまり軸トルクが出過ぎている場合に、点火時期を基本点火時期ＭＢＴＣＡＬよりリタードさせて要求軸トルクが得られるようにするためのものである。従って、（１００）式によれば、第１軸トルク推定値Ｔｒｑｅｓｔ１のほうが要求軸トルクＴｒｑｔｇｔより小さい場合があり得るが、この場合は考えない。

上記の所定値ｃｏｎｓｔ２０１は軸トルクを点火時期リタード量へと換算するための値で、予め適合する。

ステップ４７３では基本点火時期ＭＢＴＣＡＬよりこの点火時期リタード量ＲＥＴ(new)だけリタードさせたと仮定したときの軸トルク推定値を第３軸トルク推定値Ｔｒｑｅｓｔ３として算出する。

この第３軸トルク推定値Ｔｒｑｅｓｔ３の算出については図５２のフローにより説明する。図５２において図３３（Ｂ）と同一の部分には同一のステップ番号をつけている。

ここで、第３軸トルク推定値Ｔｒｑｅｓｔ３の算出方法そのものは図３３（Ｂ）で説明した第２軸トルク推定値Ｔｒｑｅｓｔ２の算出方法と同じである。図５２において図３３（Ｂ）の場合と異なるのは、図３３（Ｂ）では点火時期リタード量前回値ＲＥＴ(old)を用いたのに対して、図５２では図５１のステップ４７２で算出した点火時期リタード量ＲＥＴ(new)を用いる点だけである。すなわち、ステップ４８１では図１３のステップ４３で算出されている基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］、図５１のステップ４７２で算出されている点火時期リタード量ＲＥＴ(new)［ｄｅｇ］を読み込み、ステップ４８２で点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］の値を点火時期前回値θｉｇｎｚ［ｄｅｇＢＴＤＣ］に移すと共に、基本点火時期ＭＢＴＣＡＬから点火時期リタード量ＲＥＴ(new)を差し引いた値（つまり基本点火時期ＭＢＴＣＡＬより点火時期リタード量ＲＥＴ(new)だけリタード側の値）を点火時期θｉｇｎ［ｄｅｇＢＴＤＣ］として算出している。

図５２のステップ２３３〜２３９では図３３（Ｂ）のステップ２３３〜２３９と同じに第２軸トルク推定値Ｔｒｑｅｓｔ２を算出し、図５２のステップ４８３でこの算出した第２軸トルク推定値Ｔｒｑｅｓｔ２を第３軸トルク推定値Ｔｒｑｅｓｔ３に移す。

このようにして第３軸トルク推定値Ｔｒｑｅｓｔ３の算出を終了したら図５１のステップ４７４に戻りこの第３軸トルク推定値Ｔｒｑｅｓｔ３から要求軸トルクＴｒｑｔｇｔを差し引いて軸トルク差ｄＴｒｑ（＝Ｔｒｑｅｓｔ３−Ｔｒｑｔｇｔ）を算出し、この軸トルク差ｄＴｒｑと所定値ｃｏｎｓｔ２０２をステップ４７５で比較する。

ここで、所定値ｃｏｎｓｔ２０２は要求軸トルクの許容範囲を定める値である。軸トルク差ｄＴｒｑが所定値ｃｏｎｓｔ２０２以上あれば第３軸トルク推定値Ｔｒｑｅｓｔ３が要求軸トルクＴｒｑｔｇｔの許容範囲内にない、つまりそのときの第３軸トルク推定値Ｔｒｑｅｓｔ３によれば軸トルクが出過ぎていると判断し、ステップ４７６に進んで点火時期リタード量ＲＥＴ(new)を点火時期リタード量前回値ＲＥＴ(old)に移したあと、ステップ４７７で点火時期リタード量ＲＥＴ(new)［ｄｅｇ］を所定値だけ大きくする。つまり、次式により点火時期リタード量ＲＥＴ(new)を更新する。

ＲＥＴ(new)＝ＲＥＴ(old)＋（Ｔｒｑｅｓｔ３−Ｔｒｑｔｇｔ）
×ｃｏｎｓｔ２０１
…（１０１）
ただし、ＲＥＴ(new) ：更新後の点火時期リタード量、
ＲＥＴ(old) ：更新前の点火時期リタード量、
ｃｏｎｓｔ２０１：大きくする側への更新の程度、
この点火時期リタード量ＲＥＴ(new)の大きくする側への更新は、基本点火時期ＭＢＴＣＡＬより所定の点火時期リタード量（上記（１００）式のリタード量）だけリタードさせたと仮定したときの点火時期での軸トルク推定値である第３軸トルク推定値Ｔｒｑｅｓｔ１も要求軸トルクＴｒｑｔｇｔより大きい、つまり軸トルクが出過ぎている場合に、点火時期リタード量をさらに大きくして要求軸トルクが得られるようにするためのものである。従って、（１０１）式によれば、第３軸トルク推定値Ｔｒｑｅｓｔ３のほうが要求軸トルクＴｒｑｔｇｔより小さい場合があり得るが、この場合は考えない。

（１０１）式の所定値ｃｏｎｓｔ２０１は上記（１００）式と同じである。

この後はステップ４７３に戻ってステップ４７３〜４７７の操作を繰り返す。この操作により第３軸トルク推定値Ｔｒｑｅｓｔ３が減少して軸トルク差ｄＴｒｑが所定値ｃｏｎｓｔ２０２未満になるとステップ４７５より「ＥＮＤ」へと進んで図５１の処理を終了する。この場合に、ステップ４７３〜４７７の操作は実際に点火を行うことなく実行するものであり、ほぼ一瞬で終了する。

一方、ステップ４７５で軸トルク差ｄＴｒｑが所定値ｃｏｎｓｔ２０２未満であるとき、つまり当初より第３軸トルク推定値Ｔｒｑｅｓｔ３が要求軸トルクＴｒｑｔｇｔの許容範囲内に収まっているときにもステップ４７５より「ＥＮＤ」へと進んで図５１の処理を終了する。

このようにして点火時期リタード量ＲＥＴ(new)の算出を終了したら図２のステップ５に戻り、この点火時期リタード量ＲＥＴ(new)だけ基本点火時期ＭＢＴＣＡＬよりリタードさせて、つまり次式により点火時期指令値ＱＡＤＶ［ｄｅｇＢＴＤＣ］を設定する。

ＱＡＤＶ＝ＭＢＴＣＡＬ−ＲＥＴ(new) …（１０２）
ステップ６では点火時期リタード量ＲＥＴ(new)を点火時期リタード量前回値ＲＥＴ(old)に移す。この点火時期リタード量前回値ＲＥＴ(old)の値が前述した図３３（Ｂ）のステップ２３１、２３２で用いられる。

上記（１０２）式により設定した点火時期指令値ＱＡＤＶは、ステップ７で点火レジスタに移され、実際のクランク角がこの点火時期指令値ＱＡＤＶと一致したタイミングでエンジンコントローラ３１より一次電流を遮断する点火信号が点火コイル１３に出力される。

次に、図５３はエンジントルク補正量としての要求空気増加率Ｒｑａｉを算出するためのもので、燃料の噴射毎に実行する。

ここで、要求空気増加率Ｒｑａｉを導入したのは次の理由からである。すなわち、図５１で前述したように第１軸トルク推定値Ｔｒｑｅｓｔ１が要求軸トルクＴｒｑｔｇｔより大きい場合に点火時期がリタードされるので（図５１のステップ４７２参照）、基本点火時期ＭＢＴＣＡＬよりリタードされた点火時期で点火したときに発生する軸トルク（つまり第２軸トルク推定値Ｔｒｑｅｓｔ２）は、基本点火時期ＭＢＴＣＡＬで点火したときに発生する軸トルク（つまり第１軸トルク推定値Ｔｒｑｅｓｔ１）より低下する。この点火時期リタードに伴う軸トルクの低下を避けるため、トルク補償を行うようにしたものである。

まずステップ４９１では、図２のステップ２で算出されている第１軸トルク推定値Ｔｒｑｅｓｔ１［ｋＪ］および第２軸トルク推定値Ｔｒｑｅｓｔ２［ｋＪ］を読み込む。

ステップ４９２では、第１軸トルク推定値Ｔｒｑｅｓｔ１からの軸トルク不足量Ｔｓ［ｋＪ］を次式により算出する。

Ｔｓ＝Ｔｒｑｅｓｔ１−Ｔｒｑｅｓｔ２ …（１０３）
ステップ４９３ではこの軸トルク不足量Ｔｓと所定値ｃｏｎｓｔ２０３を比較する。ここで、所定値ｃｏｎｓｔ２０３は軸トルク不足量の許容範囲を定めるもので、予め設定しておく。軸トルク不足量Ｔｓが所定値ｃｏｎｓｔ２０３以上あるときにはステップ４９３よりステップ４９４に進み要求空気増加率Ｒｑａｉ(new)を所定値ｃｏｎｓｔ２０４だけ大きくする。つまり次式により要求空気増加率Ｒｑａｉ(new)を更新する。

Ｒｑａｉ(new)＝Ｒｑａｉ(old)＋ｃｏｎｓｔ２０４ …（１０４）
ただし、Ｒｑａｉ(new)：更新後の要求空気増加率、
Ｒｑａｉ(old)：更新前の要求空気増加率、
ｃｏｎｓｔ２０４：大きくする側への更新量（正の値）、
ここで、Ｒｑａｉ(old)の初期値としては１．０を入れておく。エンジンの運転後に初めて軸トルク不足量Ｔｓが所定値ｃｏｎｓｔ２０３以上になったときにはＲｑａｉ(old)＝１．０であることより要求空気増加率Ｒｑａｉ(new)は１．０を超える値となる。

ステップ４９５では次回に備えてＲｑａｉ(new)の値を前回値を表すＲｑａｉ(old)に移して今回の処理を終了する。

本実施形態では、燃料量と空気量については特開平９−２８７５１３号公報に記載されているのと同じ制御を行っている。この燃料量と空気量の制御について説明すると、アクセルセンサ４２により検出されるアクセル開度ＡＰＳとエンジン回転速度ＮＲＰＭとに基づいて理論空燃比を得るための基準目標吸入空気量ｔＴＰと、理論空燃比を含みリーン空燃比やリッチ空燃比での運転を可能とするための目標当量比ｔＤＭＬとを算出し、これら基準目標吸入空気量ｔＴＰと目標当量比ｔＤＭＬから目標空燃比を得るための目標吸入空気量ｔＴＰ’（＝ｔＴＰ／ｔＤＭＬ）を算出し、この目標吸入空気量ｔＴＰ’に基づいて目標絞り弁開度ｔＴＰＳを算出し、この目標絞り弁開度ｔＴＰＳが実現されるようにスロットルモータ２４で絞り弁２３の開度を制御することで空気量の制御を行いつつ、エンジン回転速度とエアフローメータ３２により検出される実吸入空気量とに基づいて算出される基本燃料噴射パルス幅ＴＰ［ｍｓ］を前記目標当量比ｔＤＭＬ及び無効パルス幅ＴＳで補正して得られる燃料噴射パルス幅ＴＩ［ｍｓ］を燃料インジェクタ２１に与えることで燃料量の制御を行っている。

こうした燃料量と空気量の制御を行うものを前提として本実施形態では、目標吸入空気量（ｔＴＰやｔＴＰ’）に上記図５３のステップ４９４（またはステップ４９６）で更新される要求空気増加率Ｒｑａｉ(new)を乗算した値を改めて目標吸入空気量（ｔＴＰやｔＴＰ’）として空気量の制御を行う。図５３のステップ４９４で更新される要求空気増加率Ｒｑａｉ(new)は１．０を超える値であるので、この１．０を超える値である要求空気増加率Ｒｑａｉ(new)が目標吸入空気量に乗算されると、乗算後の目標空気量は乗算前より大きくなり、この乗算後の目標空気量が得られるように絞り弁２３の開度が乗算前より大きくされる。これにより吸入空気量が増え、この空気量の増加を介して燃料量が増量され、同じ点火時期でもエンジンの発生する軸トルクが増加する。そして、要求空気増加率Ｒｑａｉ(new)による目標空気量の増加と並行して図２のステップ２により第２軸トルク推定値Ｔｒｑｅｓｔ２の算出が行われ、要求空気増加率Ｒｑａｉ(new)による目標空気量の増加で軸トルクが増えれば、算出される第２軸トルク推定値Ｔｒｑｅｓｔ２も大きくなり、図５３のステップ４９２で算出される軸トルク不足量Ｔｓが小さくなる。それでも、軸トルク不足量Ｔｓが所定値ｃｏｎｓｔ２０３以上であれば図５３のステップ４９４で要求空気増加率Ｒｑａｉ(new)をさらに大きくなる側に更新する。これにより、エンジンで発生する軸トルクが増え軸トルク不足量Ｔｓが減るので、図５３のステップ４９４を繰り返せばやがて軸トルク不足量Ｔｓが所定値ｃｏｎｓｔ２０３未満に収まる。このときにはそのときの要求空気増加率Ｒｑａｉ(new)の値を維持するためステップ４９３よりステップ４９６に進んで要求空気増加率前回値Ｒｑａｉ(old)を要求増加率今回値Ｒｑａｉ(new)に移して今回の処理を終了する。

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

本実施形態（請求項６に記載の発明）によれば、エンジンの運転状態を検出し、この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合（Ｘ１’、Ｘ２’）を推定し、この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定し、この推定した軸トルク（Ｔｒｑｅｓｔ１、Ｔｒｑｅｓｔ２）に基づいてトルク制御（吸入空気量の増加補正や点火時期のリタード補正）を行うように構成したので、圧力センサにより燃焼室内の圧力を直接検出することなくエンジンの発生する軸トルクを推定できる。

本実施形態（請求項１、６に記載の発明）によれば、軸トルク推定手段が、燃焼中のクランク角における燃焼質量割合の近似式（Ｘ１’、Ｘ２’）を算出する近似式算出手段（図３３（Ａ）のステップ２２７、図３３（Ｂ）のステップ２３７参照）と、この燃焼質量割合の近似式（Ｘ１’、Ｘ２’）に基づいて図示正仕事推定値（Ｗｐｏｓ１、Ｗｐｏｓ２）を算出する図示正仕事推定値算出手段（図３３（Ａ）のステップ２２８、図３３（Ｂ）のステップ２３８参照）と、この図示正仕事推定値（Ｗｐｏｓ１、Ｗｐｏｓ２）にポンピングロスＷｎｅｇ１と摩擦損失Ｗｎｅｇ２を加算した値を軸トルク推定値（Ｔｒｑｅｓｔ１、Ｔｒｑｅｓｔ２）として算出する軸トルク推定値算出手段（図３３（Ａ）のステップ２２９、図３３（Ｂ）のステップ２３９参照）とを含むので、図示正仕事推定値（Ｗｐｏｓ１、Ｗｐｏｓ２）を精度よく算出することができる。

本実施形態（請求項２、７に記載の発明）によれば、エンジンの運転状態は熱発生開始時期θｓであり、燃焼質量割合の近似式（Ｘ１’、Ｘ２’）はこの熱発生開始時期θｓと２つの定数（ｍ、ｋ）とを用いた指数関数の式であるので（図３７のステップ３０５、図４６のステップ４２５参照）、燃焼開始からの燃焼質量割合の変化をほぼ厳密に与えることができる。

本実施形態（請求項３、８に記載の発明）によれば基本点火時期ＭＢＴＣＡＬ（ＭＢＴの得られる点火時期）で点火するときと、基本点火時期ＭＢＴＣＡＬより点火時期リタード量前回値ＲＥＴ(old)だけリタード側の値（ＭＢＴの得られる点火時期よりリタードさせた点火時期）で点火したときとで燃焼質量割合の近似式を相違させるので（図３３（Ａ）のステップ２２２、２２７、図３７のステップ３０５、図３３（Ｂ）のステップ２３２、２３７、図４６のステップ４２５参照）、第１軸トルク推定値Ｔｒｑｅｓｔ１（ＭＢＴの得られる点火時期で点火したときの軸トルク推定値）と、第２軸トルク推定値Ｔｒｑｅｓｔ２（ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときの軸トルク推定値）とをいずれも精度よく算出できる。

本実施形態（請求項１１に記載の発明）によれば、軸トルク推定値算出手段が、ＭＢＴの得られる点火時期で点火したときの軸トルク推定値を第１軸トルク推定値Ｔｒｑｅｓｔ１として算出する第１軸トルク推定値算出手段（図３３（Ａ）のステップ２２９参照）と、ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときの軸トルク推定値を第２軸トルク推定値Ｔｒｑｅｓｔ２として算出する第２軸トルク推定値算出手段（図３３（Ｂ）のステップ２３９参照）とを含み、トルク制御手段が、この第２軸トルク推定値Ｔｒｑｅｓｔ２が第１軸トルク推定値Ｔｒｑｅｓｔ１と一致するように吸入空気量を増加する空気量増加手段（図５２参照）であるので、点火時期がＭＢＴの得られる点火時期よりリタードされ（図２のステップ４、５参照）、ドライバの意図しない軸トルクの低下が発生したときであっても、第１軸トルク推定値Ｔｒｑｅｓｔ１（ＭＢＴの得られる点火時期で点火したときの軸トルク）と同じ軸トルクが得られる。

第１軸トルク推定値Ｔｒｑｅｓｔ１（ＭＢＴの得られる点火時期で点火したときの軸トルク推定値）が要求軸トルクＴｒｑｔｇｔを超えているときには軸トルクが出過ぎているのであり、この場合に本実施形態（請求項１２に記載の発明）によれば、トルク制御手段が、第１軸トルク推定値Ｔｒｑｅｓｔ１が要求軸トルクＴｒｑｔｇｔを超えるとき、この第１軸トルク推定値Ｔｒｑｅｓｔ１が要求軸トルクＴｒｑｔｇｔと一致するように所定の点火時期リタード量ＲＥＴ(new)だけ点火時期をリタード側に補正する点火時期補正手段（図５１のステップ４７１、４７２参照）を含むので、第１軸トルク推定値Ｔｒｑｅｓｔ１が要求軸トルクＴｒｑｔｇｔを超えるときであってもドライバの意図通りの軸トルクを発生させることができる。

前記所定の点火時期リタード量ＲＥＴ(new)だけ点火時期をリタードさせたときの軸トルクが要求軸トルクＴｒｑｔｇｔを超えているときにはまだ軸トルクの出過ぎであり、この場合に、本実施形態（請求項１３に記載の発明）によれば、軸トルク推定値算出手段が、基本点火時期ＭＢＴＣＡＬ（ＭＢＴの得られる点火時期）より前記所定の点火時期リタード量ＲＥＴ(new)だけ点火時期をリタードさせたと仮定したときの軸トルク推定値を第３軸トルク推定値Ｔｒｑｅｓｔ３として算出する第３軸トルク推定値算出手段（図５１のステップ４７３、図５２参照）を含み、トルク制御手段が、この第３軸トルク推定値Ｔｒｑｅｓｔ３が要求軸トルクＴｒｑｔｇｔを超えるとき、この第３軸トルク推定値Ｔｒｑｅｓｔ３が要求軸トルクＴｒｑｔｇｔと一致するように点火時期をリタード側に補正する点火時期補正手段（図５１のステップ４７３〜４７７参照）を含むので、第３軸トルク推定値Ｔｒｑｅｓｔ３がまだ要求軸トルクＴｒｑｔｇｔを超えるときであってもドライバの意図通りの軸トルクを発生させることができる。

実施形態では、図５１に示したように、第１軸トルク推定値Ｔｒｑｅｓｔ１を要求軸トルクＴｒｑｔｇｔへと下げるための点火時期リタード量だけリタードさせたと仮定したときの軸トルク推定値を求め、この軸トルク推定値がまだ要求軸トルクＴｒｑｔｇｔより大きいときには点火時期リタード量を大きくなる側に更新する場合で説明したが、簡単には、図５１においてステップ４７３〜４７７の操作を省略した態様も考え得る。

実施形態では、エンジンの運転状態は熱発生開始時期θｓであり、燃焼質量割合の近似式がこの熱発生開始時期θｓと２つの定数（ｍ、ｋ）とを用いた指数関数の式である場合で説明したが、エンジンの運転状態は燃焼開始から所定クランク角までの燃焼期間であり、この燃焼期間を複数に分割して算出する燃焼期間算出手段を備えさせ、この複数の燃焼期間に基づいて燃焼質量割合の近似式を設定することもできる（第２実施形態）。

この第２実施形態について説明すると、いま燃焼質量割合を求めるためのクランク角として、圧縮上死点ＴＤＣを基準のゼロとしこれより遅角側をプラス、これより進角側をマイナスに採ったクランク角Θ［ｄｅｇＡＴＤＣ］を考える。このとき、基本点火時期ＭＢＴＣＡＬで点火したときの燃焼質量割合ＢＲは基本点火時期ＭＢＴＣＡＬ、第１点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ１、第１初期燃焼期間ＢＵＲＮ１ｒｅｆ１、第１主燃焼期間ＢＵＲＮ２ｒｅｆ１を用いて次のような一次式で近似できる。

第１燃焼遅れ期間；
ＢＲ＝０ …（１０５）
第１初期燃焼期間；
ＢＲ＝ＳＳ１×（Θ＋ＭＢＴＣＡＬ−ＩＧＮＤＥＡＤｒｅｆ１）
…（１０６）
第１主燃焼期間 ；
ＢＲ＝０．０２＋ＳＳ２×（Θ＋ＭＢＴＣＡＬ−ＩＧＮＤＥＡＤｒｅｆ１
−ＢＵＲＮ１ｒｅｆ１）
…（１０７）
ただし、ＳＳ１：０．０２／ＢＵＲＮ１ｒｅｆ１、
ＳＳ２：０．５８／ＢＵＲＮ２ｒｅｆ１、
同様にして、基本点火時期ＭＢＴＣＡＬより点火時期リタード量前回値ＲＥＴ(old)だけリタードさせた点火時期で点火したときの燃焼質量割合ＢＲは基本点火時期ＭＢＴＣＡＬ、第２点火無駄時間相当クランク角ＩＧＮＤＥＡＤｒｅｆ２、第２初期燃焼期間ＢＵＲＮ１ｒｅｆ２、第２主燃焼期間ＢＵＲＮ２ｒｅｆ２を用いて次のような一次式で近似できる。

第２燃焼遅れ期間；
ＢＲ＝０ …（１０８）
第２初期燃焼期間；
ＢＲ＝ＳＳ１×（Θ＋ＭＢＴＣＡＬ−ＩＧＮＤＥＡＤｒｅｆ２）
…（１０９）
第２主燃焼期間 ；
ＢＲ＝０．０２＋ＳＳ２×（Θ＋ＭＢＴＣＡＬ−ＩＧＮＤＥＡＤｒｅｆ２
−ＢＵＲＮ１ｒｅｆ２）
…（１１０）
ただし、ＳＳ１：０．０２／ＢＵＲＮ１ｒｅｆ２、
ＳＳ２：０．５８／ＢＵＲＮ２ｒｅｆ２、
図５４は、基本点火時期ＭＢＴＣＡＬで点火したときの燃焼質量割合と、基本点火時期ＭＢＴＣＡＬより点火時期リタード量前回値ＲＥＴ(old)だけリタードさせた点火時期で点火したときの燃焼質量割合の一次式近似波形をモデル的に示したものである。

このように第２実施形態（請求項４、９に記載の発明）によれば、エンジンの運転状態が燃焼開始から所定クランク角までの燃焼期間であり、この燃焼期間を初期燃焼期間と主燃焼期間の２つ（複数）に分割して算出する燃焼期間算出手段（図３３（Ａ）のステップ２２４〜２２６、図３３（Ｂ）のステップ２３４〜２３６参照）を備え、この複数の燃焼期間（ＢＵＲＮ１ｒｅｆ１、ＢＵＲＮ１ｒｅｆ２、ＢＵＲＮ２ｒｅｆ１、ＢＵＲＮ２ｒｅｆ２、ＩＧＮＤＥＡＤｒｅｆ１、ＩＧＮＤＥＡＤｒｅｆ２）に基づいて燃焼質量割合の近似式（上記（１０５）〜（１１０）式参照）を設定するので、燃焼質量割合の変化を簡単な一次式で近似することができ、燃焼質量割合の近似式の設定が容易になる。

実施形態では、エンジンの運転状態を検出する運転状態検出手段と、この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合を推定する推定手段と、この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定する軸トルク推定手段と、この推定した軸トルクに基づいてトルク制御を行うトルク制御手段とを備える場合で説明したが（請求項６に記載の発明）、エンジンの運転状態を検出する運転状態検出手段と、この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合を推定する推定手段と、この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定する軸トルク推定手段とで軸トルク演算装置を構成することができる。

請求項１、６に記載の発明において、運転状態検出手段の機能は図３７のステップ３０２、図４６のステップ４２２または図３３（Ａ）のステップ２２４〜２２６、図３３（Ｂ）のステップ２３４〜２３６により、推定手段の機能は図３３（Ａ）のステップ２２７、図３３（Ｂ）のステップ２３７により、軸トルク推定手段の機能は図３３（Ａ）のステップ２２７〜２２９、図３３（Ｂ）のステップ２３７〜２３９により、トルク制御手段の機能は図５１または図５３によりそれぞれ果たされている。

一実施形態のエンジンの制御システム図。 点火時期制御を説明するためのフローチャート。 燃焼室の圧力変化図。 燃焼質量割合の変化を説明する特性図。 物理量の算出を説明するためのフローチャート。 エンジンのクランクシャフトとコネクティングロッドの位置関係を説明するダイアグラム。 水温補正係数の特性図。 当量比補正係数の特性図。 基準クランク角の特性図。 初期燃焼期間の算出を説明するためのフローチャート。 温度上昇率の特性図。 主燃焼期間の算出を説明するためのフローチャート。 基本点火時期の算出を説明するためのフローチャート。 内部不活性ガス率の算出を説明するためのフローチャート。 内部不活性ガス量の算出を説明するためのフローチャート。 ＥＶＣ時不活性ガス量の算出を説明するためのフローチャート。 オーバーラップ中吹き返し不活性ガス量の算出を説明するためのフローチャート。 過給判定フラグ、チョーク判定フラグの設定を説明するためのフローチャート。 過給無しかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 過給無しかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 過給有りかつチョーク無し時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 過給有りかつチョーク有り時のオーバーラップ中吹き返し不活性ガス流量の算出を説明するためのフローチャート。 排気弁閉時期における燃焼室容積の特性図。 不活性ガスのガス定数の特性図。 オーバーラップ中の積算有効面積の特性図。 オーバーラップ中の積算有効面積の説明図。 不活性ガスの比熱比の特性図。 混合気の比熱比の特性図。 点火無駄時間相当クランク角の算出を説明するためのフローチャート。 点火無駄時間基本値の特性図。 温度補正係数の特性図。 空気過剰率補正係数の特性図。 第１軸トルク推定値の算出を説明するためのフローチャート。 第２軸トルク推定値の算出を説明するためのフローチャート。 第１初期燃焼期間の算出を説明するためのフローチャート。 第１点火無駄時間相当クランク角の算出を説明するためのフローチャート。 第１主燃焼期間の算出を説明するためのフローチャート。 燃焼質量割合の第１近似式の算出を説明するためのフローチャート。 第１図示正仕事推定値の算出を説明するためのフローチャート。 第１図示正仕事推定値の算出を説明するためのフローチャート。 冷却損失の算出を説明するためのフローチャート。 第１軸トルク推定値の算出を説明するためのフローチャート。 摩擦損失の特性図。 第２初期燃焼期間の算出を説明するためのフローチャート。 第２点火無駄時間相当クランク角の算出を説明するためのフローチャート。 第２主燃焼期間の算出を説明するためのフローチャート。 燃焼質量割合の第２近似式の算出を説明するためのフローチャート。 第２図示正仕事推定値の算出を説明するためのフローチャート。 第２図示正仕事推定値の算出を説明するためのフローチャート。 第２初期燃焼期間の算出を説明するためのフローチャート。 要求軸トルクの特性図。 点火時期リタード量の算出を説明するためのフローチャート。 第３軸トルク推定値の算出を説明するためのフローチャート。 要求空気増加率の算出を説明するためのフローチャート。 燃焼質量割合の変化を示す波形図。

符号の説明

１ エンジン
５ 燃焼室
１１ 点火装置（火花点火手段）
２１ 燃料インジェクタ
３１ エンジンコントローラ

Claims (13)

1. エンジンの運転状態を検出する運転状態検出手段と、
   この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合を推定する推定手段と、
   この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定する軸トルク推定手段と
   を備え、
   前記軸トルク推定手段は、
   この燃焼中のクランク角における燃焼質量割合の近似式を算出する近似式算出手段と、
   この燃焼質量割合の近似式に基づいて図示正仕事推定値を算出する図示正仕事推定値算出手段と、
   この図示正仕事推定値にポンピングロスと摩擦損失を加算した値を軸トルク推定値として算出する軸トルク推定値算出手段と
   を含むことを特徴する軸トルク演算装置。
2. 前記エンジンの運転状態は熱発生開始時期であり、前記燃焼質量割合の近似式は２つの定数を用いた熱発生開始時期の指数関数の式であることを特徴する請求項１に記載の軸トルク演算装置。
3. ＭＢＴの得られる点火時期で点火するときと、ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときとで前記燃焼質量割合の近似式を相違させることを特徴する請求項２に記載の軸トルク演算装置。
4. 前記エンジンの運転状態は燃焼開始から所定クランク角までの燃焼期間であり、この燃焼期間を複数に分割して算出する燃焼期間算出手段を備え、この複数の燃焼期間に基づいて前記燃焼質量割合の近似式を設定することを特徴する請求項１に記載の軸トルク演算装置。
5. ＭＢＴの得られる点火時期で点火したときと、ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときとで前記燃焼質量割合の近似式を別々に設定することを特徴する請求項４に記載の軸トルク演算装置。
6. エンジンの運転状態を検出する運転状態検出手段と、
   この検出した運転状態に基づいて燃焼室内のガスの燃焼中のクランク角における燃焼質量割合を推定する推定手段と、
   この推定した燃焼質量割合に基づいてエンジンの発生する軸トルクを推定する軸トルク推定手段と、
   この推定した軸トルクに基づいてトルク制御を行うトルク制御手段と
   を備え、
   前記軸トルク推定手段は、
   この燃焼中のクランク角における燃焼質量割合の近似式を算出する近似式算出手段と、
   この燃焼質量割合の近似式に基づいて図示正仕事推定値を算出する図示正仕事推定値算出手段と、
   この図示正仕事推定値にポンピングロスと摩擦損失を加算した値を軸トルク推定値として算出する軸トルク推定値算出手段と
   を含むことを特徴とするエンジンの制御装置。
7. 前記エンジンの運転状態は熱発生開始時期であり、前記燃焼質量割合の近似式は２つの定数を用いた熱発生開始時期の指数関数の式であることを特徴とする請求項６に記載のエンジンの制御装置。
8. ＭＢＴの得られる点火時期で点火するときと、ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときとで前記燃焼質量割合の近似式を相違させることを特徴とする請求項７に記載のエンジンの制御装置。
9. 前記エンジンの運転状態は燃焼開始から所定クランク角までの燃焼期間であり、この燃焼期間を複数に分割して算出する燃焼期間算出手段を備え、この複数の燃焼期間に基づいて前記燃焼質量割合の近似式を設定することを特徴とする請求項６に記載のエンジンの制御装置。
10. ＭＢＴの得られる点火時期で点火したときと、ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときとで前記燃焼質量割合の近似式を別々に設定することを特徴とする請求項９に記載のエンジンの制御装置。
11. 前記軸トルク推定値算出手段は、
    ＭＢＴの得られる点火時期で点火したときの軸トルク推定値を第１軸トルク推定値として算出する第１軸トルク推定値算出手段と、
    ＭＢＴの得られる点火時期よりリタードさせた点火時期で点火したときの軸トルク推定値を第２軸トルク推定値として算出する第２軸トルク推定値算出手段と
    を含み、
    前記トルク制御手段は、この第２軸トルク推定値が前記要求要求軸トルクと一致するように吸入空気量を増加する空気量増加手段であることを特徴とする請求項６に記載のエンジンの制御装置。
12. 前記トルク制御手段は、前記第１軸トルク推定値が要求軸トルクを超えるとき、この第１軸トルク推定値が要求軸トルクと一致するように所定の点火時期リタード量だけ点火時期をリタード側に補正する点火時期補正手段を含むことを特徴とする請求項１１に記載のエンジンの制御装置。
13. 前記軸トルク推定値算出手段は、ＭＢＴの得られる点火時期より前記所定の点火時期リタード量だけ点火時期をリタードさせたと仮定したときの軸トルク推定値を第３軸トルク推定値として算出する第３軸トルク推定値算出手段を含み、前記トルク制御手段は、この第３軸トルク推定値が前記要求軸トルクを超えるとき、この第３軸トルク推定値が要求軸トルクと一致するように点火時期をリタード側に補正する点火時期補正手段を含むことを特徴とする請求項１２に記載のエンジンの制御装置。

Images