(57)【特許請求の範囲】
【請求項１】
吸気弁のリフト特性とエンジン回転速度と吸気弁閉時期からタンブル形成角までのクランク角区間とに基づいて燃焼室内のガス流動としてのタンブル強度を推定するタンブル強度推定処理手順と、
この燃焼室内のガス流動に基づいて燃焼室内ガスの乱れ強さを算出する乱れ強さ算出処理手順と、
この燃焼室内ガスの乱れ強さに基づいて燃焼室内ガスの乱流状態での燃焼速度である乱流燃焼速度を算出する乱流燃焼速度算出処理手順と、
この乱流燃焼速度に基づいてＭＢＴの得られる基本点火時期を算出する基本点火時期算出処理手順と
を含むことを特徴とするエンジンの点火時期制御方法。
【請求項２】
吸気弁のバルブリフト量及び作動角を連続的に可変制御する可変バルブ機構及び吸気弁のバルブタイミングを進遅角する可変バルブタイミング機構を備え、
運転条件に応じて吸気弁の閉時期とバルブリフト量とをこれらの機構を用いて制御する制御処理手順を含むことを特徴とする請求項１に記載のエンジンの点火時期制御方法。
【請求項３】
前記乱流燃焼速度を少なくとも吸気弁のリフト特性に基づいて算出することを特徴とする請求項１に記載のエンジンの点火時期制御方法。
【請求項４】
前記乱流燃焼速度を少なくとも吸気弁の閉時期に基づいて算出することを特徴とする請求項１に記載のエンジンの点火時期制御方法。
【請求項５】
エンジンの圧縮工程開始時期を検出する圧縮工程開始時期検出装置またはエンジンの圧縮工程開始時期を推定する圧縮工程開始時期推定装置を備え、
前記乱流燃焼速度を少なくともこの圧縮工程開始時期に基づいて算出することを特徴とする請求項１に記載のエンジンの点火時期制御方法。
【請求項６】
エンジンの圧縮工程開始時期を検出する圧縮工程開始時期検出装置またはエンジンの圧縮工程開始時期を推定する圧縮工程開始時期推定装置を備え、
前記乱流燃焼速度をこの圧縮工程開始時期に基づいても算出することを特徴とする請求項３に記載のエンジンの点火時期制御方法。
【請求項７】
前記乱流燃焼速度を燃焼室内のガス量に基づいても算出することを特徴とする請求項３から６までのいずれか一つに記載のエンジンの点火時期制御方法。
【請求項８】
燃焼室内ガスの層流状態での燃焼速度である層流燃焼速度を算出する層流燃焼速度算出処理手順を含み、
前記乱流燃焼速度をこの層流燃焼速度に基づいても算出することを特徴とする請求項３から７までのいずれか一つに記載のエンジンの点火時期制御方法。
【請求項９】
燃焼室内の吸気流動を変化させ得る可変吸気流動装置を備え、
前記乱流燃焼速度をこの可変吸気流動装置の状態に応じても算出することを特徴とする請求項３から８までのいずれか一つに記載のエンジンの点火時期制御方法。
【請求項１０】
吸気弁のリフト特性とエンジン回転速度と吸気弁閉時期からタンブル形成角までのクランク角区間とに基づいて燃焼室内のガス流動としてのタンブル強度を推定するタンブル強度推定手段と、
この燃焼室内のガス流動に基づいて燃焼室内ガスの乱れ強さを算出する乱れ強さ算出手段と、
この燃焼室内ガスの乱れ強さに基づいて燃焼室内ガスの乱流状態での燃焼速度である乱流燃焼速度を算出する乱流燃焼速度算出手段と、
この乱流燃焼速度に基づいてＭＢＴの得られる基本点火時期を算出する基本点火時期算出手段と
を備えることを特徴とするエンジンの点火時期制御装置。
【請求項１１】
吸気弁のバルブリフト量及び作動角を連続的に可変制御する可変バルブ機構及び吸気弁のバルブタイミングを進遅角する可変バルブタイミング機構を備え、
運転条件に応じて吸気弁の閉時期とバルブリフト量とをこれらの機構を用いて制御する制御手段を含むことを特徴とする請求項１０に記載のエンジンの点火時期制御装置。
【請求項１２】
前記乱流燃焼速度を少なくとも吸気弁のリフト特性に基づいて算出すことを特徴とする請求項１０に記載のエンジンの点火時期制御装置。
【請求項１３】
前記乱流燃焼速度を少なくとも吸気弁の閉時期に基づいて算出することを特徴とする請求項１０に記載のエンジンの点火時期制御装置。
【請求項１４】
エンジンの圧縮工程開始時期を検出する圧縮工程開始時期検出装置またはエンジンの圧縮工程開始時期を推定する圧縮工程開始時期推定装置を備え、
前記乱流燃焼速度を少なくともこの圧縮工程開始時期に基づいて算出することを特徴とする請求項１０に記載のエンジンの点火時期制御装置。
【請求項１５】
エンジンの圧縮工程開始時期を検出する圧縮工程開始時期検出装置またはエンジンの圧縮工程開始時期を推定する圧縮工程開始時期推定装置を備え、
前記乱流燃焼速度をこの圧縮工程開始時期に基づいても算出することを特徴とする請求項１２に記載のエンジンの点火時期制御装置。
【請求項１６】
前記乱流燃焼速度を燃焼室内のガス量に基づいても算出することを特徴とする請求項１２から１５までのいずれか一つに記載のエンジンの点火時期制御装置。
【請求項１７】
燃焼室内ガスの層流状態での燃焼速度である層流燃焼速度を算出する層流燃焼速度算出手段を備え、
前記乱流燃焼速度をこの層流燃焼速度に基づいても算出することを特徴とする請求項１２から１６までのいずれか一つに記載のエンジンの点火時期制御装置。
【請求項１８】
燃焼室内の吸気流動を変化させ得る可変吸気流動装置を備え、
前記乱流燃焼速度をこの可変吸気流動装置の状態に応じても算出することを特徴とする請求項１２から１７までのいずれか一つに記載のエンジンの点火時期制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、エンジン（内燃機関）の点火時期制御方法及びエンジンの点火時期制御装置に関する。
【背景技術】
【０００２】
燃焼室ガスの乱れ強さを、吸気流速に基づいて算出し、この燃焼室内ガスの乱れ強さとエンジン回転速度とに基づいて乱流燃焼速度を算出し、この乱流燃焼速度に基づいてＭＢＴの得られる点火時期を算出するものがある（特許文献１参照）。
【０００３】
ここで、燃焼室内ガスの乱れ強さとは、燃焼室内のガス流動に伴う乱れのうちの１つの特性値、また、乱流燃焼速度とは、燃焼室内ガスの乱流状態での燃焼速度のことである。
【特許文献１】特開２００４−３０１０４５公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
ところで、燃焼室内ガスの乱れ強さによって燃焼室内における乱流燃焼速度が変化し、燃焼室内ガスの乱れ強さには、燃焼室内のガス流動が大きく影響する。従って、本来なら、この燃焼室内のガス流動そのものに基づいて燃焼室内ガスの乱れ強さを推定すべきである。
【０００５】
しかしながら、燃焼室内のガス流動そのものに基づいて燃焼室ガスの乱れ強さを推定するものは今まで開示されていない。
【０００６】
一方、上記特許文献１の技術は、燃焼室内ガスの乱れ強さと、吸気流速との間に相関があるとするものであるが、吸気流速は必ずしも燃焼室内のガス流動を正確に表すものでない。というのも、スワール強度、タンブル強度といった燃焼室内のガス流動は、後述する（３）式のように、理論的には燃焼室内のさまざまな場所に分布する小渦の角運動量（ｒ×ｍ・ｖ；ただし、ｒはシリンダの中心から小渦までの距離、ｖは小渦のガス流速、ｍは小渦のガス質量）の総和で表されるので、吸気流速だけでは小渦の角運動量を正確に表せないからである。このため、吸気流速から推定される燃焼室内のガス流動が実際の燃焼室内のガス流動よりずれる分だけ、燃焼室内ガスの乱れ強さの算出精度が落ち、そのため、算出される乱流燃焼速度の実際の乱流燃焼速度からの誤差が大きくなり、ＭＢＴの得られる点火時期の推定精度が悪化してしまう。
【０００７】
そこで本発明は、燃焼室内ガスの乱れ強さの算出方法及び算出手段を改良した点火時期制御方法及び点火時期制御装置を提供することを目的とする。
【課題を解決するための手段】
【０００８】
本発明は、吸気弁のリフト特性とエンジン回転速度と吸気弁閉時期からタンブル形成角までのクランク角区間とに基づいて燃焼室内のガス流動としてのタンブル強度を推定し、この燃焼室内のガス流動に基づいて燃焼室内ガスの乱れ強さを算出し、この燃焼室内ガスの乱れ強さに基づいて乱流燃焼速度を算出し、この乱流燃焼速度に基づいてＭＢＴの得られる基本点火時期を算出するように構成する。
【発明の効果】
【０００９】
タンブル強度といった燃焼室内のガス流動は、燃焼室内のさまざまな場所に分布する小渦の角運動量の総和で表されるので、例えばシリンダの中心からのｉ番目の小渦までの距離をｒi、ｉ番目の小渦のガス質量をｍi、ｉ番目の小渦のガス流速をｖiとすると、後述する（３）式によりタンブル強度Ｉｔ（燃焼室内のガス流動）が表されるのであるが、燃焼室内の全ての小渦についてこれらの値ｍi、ｒi、ｖiをエンジンコントローラにおいてオンラインで算出することは現在のところ不可能である。
【００１０】
そこで、本発明者は、タンブルが圧縮工程中に小渦に変換されることにより、燃焼室内ガスの乱れが強くなるので、タンブル強度が強いほど燃焼室内ガスの乱れが強くなると考え、タンブル強度（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さを推定することとした。すなわち、本発明によれば、タンブルといった燃焼室内のガス流動に基づいて燃焼室内ガスの乱れ強さ（ｕ）を算出し、この燃焼室内ガスの乱れ強さ（ｕ）に基づいて、乱流燃焼速度（ＳＴ）を算出し、この乱流燃焼速度（ＳＴ）に基づいてＭＢＴの得られる基本点火時期（ＭＢＴＣＡＬ）を算出するので、燃焼室内のガス流動の影響を受けて乱流燃焼速度（ＳＴ）が変動しても、その変動する乱流燃焼速度（ＳＴ）を精度よく算出できることから、ＭＢＴの得られる基本点火時期（ＭＢＴＣＡＬ）の予測精度を向上させることができる。
【発明を実施するための最良の形態】
【００１１】
以下、図面に基づき本発明の実施形態について説明する。
【００１２】
図１は、エンジンの点火時期制御方法の実施に直接使用するエンジンの点火時期制御装置の概略構成を示している。
【００１３】
空気は吸気コレクタ２に蓄えられた後、吸気マニホールド３を介して各気筒の燃焼室５に導入される。燃料は各気筒の吸気ポート４に配置された燃料インジェクタ２１より噴射供給される。空気中に噴射された燃料は気化しつつ空気と混合してガス（混合気）を作り、燃焼室５に流入する。この混合気は吸気弁１５が閉じることで燃焼室５内に閉じこめられ、ピストン６の上昇によって圧縮される。
【００１４】
この圧縮混合気に対して高圧火花により点火を行うため、パワートランジスタ内蔵の点火コイルを各気筒に配した電子配電システムの点火装置１１を備える。すなわち、点火装置１１は、バッテリからの電気エネルギーを蓄える点火コイル１３と、点火コイル１３の一次側への通電、遮断を行うパワートランジスタ（図示しない）と、燃焼室５の天井に設けられ点火コイル１３の一次電流の遮断によって点火コイル１３の二次側に発生する高電圧を受けて、火花放電を行う点火プラグ１４とからなっている。
【００１５】
圧縮上死点より少し手前で点火プラグ１４により火花が飛ばされ圧縮混合気に着火されると、火炎が広がりやがて爆発的に燃焼し、この燃焼によるガス圧がピストン６を押し下げる仕事を行う。この仕事はクランクシャフト７の回転力として取り出される。燃焼後のガス（排気）は排気弁１６が開いたときに排気通路８へと排出される。
【００１６】
排気通路８には三元触媒９、１０を備える。三元触媒９、１０は排気の空燃比が理論空燃比を中心とした狭い範囲（ウインドウ）にあるとき、排気中に含まれるＨＣ、ＣＯ、ＮＯｘといった有害三成分を同時に効率よく除去できる。空燃比は吸入空気量と燃料量の比であるので、エンジンの１サイクル（４サイクルエンジンではクランク角で７２０°区間）当たりに燃焼室５に導入される吸入空気量と、燃料インジェクタ２１からの燃料噴射量との比が理論空燃比となるように、エンジンコントローラ３１ではエアフローメータ３２からの吸入空気流量の信号とクランク角センサ（３３、３４）からの信号に基づいて燃料インジェクタ２１からの燃料噴射量を定めると共に、三元触媒９の上流に設けたＯ2センサ３５からの信号に基づいて空燃比をフィードバック制御している。
【００１７】
吸気コレクタ２の上流には絞り弁２３がスロットルモータ２４により駆動される、いわゆる電子制御スロットル２２を備える。運転者が要求するトルクはアクセルペダル４１の踏み込み量（アクセル開度）に現れるので、エンジンコントローラ３１ではアクセルセンサ４２からの信号に基づいて目標トルクを定め、この目標トルクを実現するための目標空気量を定め、この目標空気量が得られるようにスロットルモータ２４を介して絞り弁２３の開度を制御する。
【００１８】
吸気弁１５、排気弁１６は、クランクシャフト７を動力源として、各々吸気側カムシャフト２５及び排気側カムシャフト２６に設けられたカムの動作により開閉駆動される。吸気側には、吸気弁１５のバルブリフト量及び作動角を連続的に可変制御する多節リンク状の機構で構成される可変バルブ機構（以下、「ＶＥＬ機構」という。）２８を備える。このＶＥＬ機構２８には吸気弁１５のバルブリフト量及び作動角を検出する作動角センサ（図示しない）が併設されている。
【００１９】
同じく吸気側には、クランクシャフト７と吸気側カムシャフト２５との回転位相差を連続的に可変制御して、吸気弁１５のバルブタイミングを進遅角する可変バルブタイミング機構（以下、「ＶＴＣ機構」という。）２７を備える。また、吸気側カムシャフト２５の他端には吸気側カムシャフト２５の回転位置を検出するためのカム角センサ３４が併設されている。
【００２０】
これらＶＥＬ機構２８及びＶＴＣ機構２７（可変動弁装置）の具体的な構成は特開２００３−３８７２号公報により公知であるので、その詳しい説明は省略する。
【００２１】
ＶＥＬ機構２８、ＶＴＣ機構２７の各アクチュエータに指令して、吸気弁１５のリフト特性（バルブタイミング（開閉時期）や吸気弁１５のバルブリフト量）を変えると燃焼室５に残留する不活性ガスの量が変化する。燃焼室５内の不活性ガスの量が増えるほどポンピングロスが減って燃費がよくなるので、運転条件によりどのくらいの不活性ガスが燃焼室５内に残留したらよいかを目標吸気弁閉時期や目標バルブリフト量にして予め定めており、エンジンコントローラ３１ではそのときの運転条件（エンジンの負荷と回転速度）より目標吸気弁閉時期と目標バルブリフト量とを定め、それら目標値が得られるようにＶＴＣ機構２７及びＶＥＬ機構２８の各アクチュエータを介して吸気弁１５の閉時期とバルブリフト量とを制御する。
【００２２】
吸気温度センサ４３からの吸気温度の信号、吸気圧力センサ４４からの吸気圧力の信号、排気温度センサ４５からの排気温度の信号も入力されるエンジンコントローラ３１では、パワートランジスタ１３を介して点火プラグ１４の一次側電流の遮断時期である点火時期を制御する。
【００２３】
さて、ＭＢＴ（最大トルクの得られる最小進角値）で混合気に点火した場合に、混合気の燃焼圧力が最大値Ｐｍａｘとなるクランク角θｐｍａｘを基準クランク角とすると、基準クランク角は燃焼方式によらずほぼ一定である。また、燃焼室内における燃焼解析によれば、燃焼室５に供給された燃料に対する燃料質量の比率を表す燃焼質量割合は点火時に０％であり、完全燃焼によって１００％に達する。そして、基準クランク角における燃焼質量割合は一定で約６０％であるとして、そのときの運転条件（エンジンの負荷と回転速度）より燃焼速度を求め、この燃焼速度に基づいて燃焼期間ＢＵＲＮを算出し、基準クランク角よりこの燃焼期間ＢＵＲＮと着火遅れ時間相当角ＩＧＮＤＥＡＤとの合計のクランク角区間だけ進角側のクランク角位置を基本点火時期ＭＢＴＣＡＬとして算出する点火時期制御方法を提案している（特開２００４−３３２６４７号公報参照）。
【００２４】
燃焼解析に基づくこうした点火時期制御方法（この点火時期制御方法を、以下「先行点火時期制御方法」という。）は吸気弁１５のバルブリフト量及び吸気弁１５の開閉タイミングが変化しないエンジン（コンベンショナルエンジン）を対象に本出願人が開発してきた経過があるので、現在でもコンベンショナルエンジンに適用する限り先行点火時期制御方法でなんら問題ないのであるが、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンに対しても、この先行点火時期制御方法を適用したとき、上記の基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期と合わないことが判明している。その理由は、ＶＥＬ機構２８及びＶＴＣ機構２７が非作動状態であるときに点火時期制御に用いる各種の値を適合していれば、当然ながら、ＶＥＬ機構２８及びＶＴＣ機構２７が非作動状態であるときに基本点火時期ＭＢＴＣＡＬはＭＢＴの得られる点火時期と一致する。
【００２５】
この状態つまり同じ運転条件においてＶＥＬ機構２８やＶＴＣ機構２７を作動させてリフト特性（吸気弁１５のバルブリフト量や吸気弁閉時期）を変化させると、燃焼室５内のガス流動であるタンブルやスワールの各強度が変化し、これに伴って燃焼室内ガスの乱流状態での燃焼速度である乱流燃焼速度が変化する。すると、ＭＢＴの得られる点火時期も変化する。しかしながら、ＶＥＬ機構２８やＶＴＣ機構２７を作動させたからといって運転条件が同じであるため、基本点火時期ＭＢＴＣＡＬは変化しない。従って、ＶＥＬ機構２８やＶＴＣ機構２７を作動させた途端に基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期からずれて燃費が悪くなる。
【００２６】
例えば、ＶＥＬ機構２８の作動で吸気弁１５のバルブリフト量が大きい状態から小さい状態に切換えられたとすると、切換前よりタンブルやスワールが形成されにくくなり乱流燃焼速度が遅くなって燃焼期間が長引く。従って、同じ運転条件で考えると、バルブリフト量が小さい状態であるときには、バルブリフト量が大きい状態のときよりＭＢＴの得られる点火時期が進角側に移動するのであるが、バルブリフト量が小さい状態のときにも、バルブリフト量が大きい状態のときと同じ燃焼期間ＢＵＲＮを算出したのでは、実際より短い燃焼期間を算出してしまうことになり、基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期よりも遅すぎることとなる。
【００２７】
ＶＴＣ機構２７の作動で吸気弁閉時期が遅い状態から早い状態へと進角される場合も同様である。すなわち、吸気弁閉時期が遅い状態から早い状態へと進角されると、進角される前よりタンブルやスワールが形成されにくいために乱流燃焼速度が遅くなって燃焼期間が長くなる。従って、同じ運転条件で考えると、吸気弁閉時期が早い状態であるときには、吸気弁閉時期が遅い状態のときよりＭＢＴの得られる点火時期が進角側に移動するのであるが、吸気弁閉時期が早いときにも吸気弁閉時期が遅い状態のときと同じ燃焼期間ＢＵＲＮを算出したのでは、実際より短い燃焼期間を算出してしまうことになり、基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期より遅すぎてしまう。
【００２８】
そこで本実施形態では、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象としていても、物理モデルに従った制御構造をもち、可能な限り適合の不要な制御とするため、コンベンショナルエンジンに適用している先行点火時期制御方法に対して〈１〉燃焼速度の算出方法、〈２〉基準クランク角の設定方法、〈３〉着火遅れ時間の算出方法を変更する。
【００２９】
ここで、上記〈１〉だけでなく、〈２〉と〈３〉をも追加しているのは、制御方法を再検討した結果である。まず、上記〈２〉をも追加しているのは次の理由による。すなわち、コンベンショナルエンジンに適用している先行点火時期制御方法に対してこれまで正しいとされてきた、基準クランク角と燃焼質量割合が６０％のときのクランク角位置との関係を、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンについて実験してみると、基準クランク角と燃焼質量割合が６０％のときのクランク角位置との間に大きなずれ（クランク角差）があり、そのずれがそのまま点火時期算出値（ＭＢＴＣＡＬ）の推定誤差になってしまうため、基準クランク角位置の設定方法を改める必要が生じたためである。
【００３０】
上記〈３〉をも追加しているのは次の理由による。〈３〉の着火遅れ時間とは、点火タイミングを起点として燃焼質量割合が０％である間の時間のことであるが、この着火遅れ時間の推定誤差が大きい。その理由は、燃焼解析装置によって燃焼質量割合が０％であるクランク角位置を計測したときの計測誤差が大きいために、この計測した値を正にして着火遅れ時間ＤＥＡＤＴＩＭＥの適合を行うと、燃焼期間ＢＵＲＮの推定誤差にこの計測誤差も含まれ誤差が大きくなってしまうためである。従って、燃焼解析装置による計測誤差が入ってこないように、着火遅れ時間の算出方法を改めるようにしたものである。
【００３１】
以下、項を分けて説明する。
〈１〉燃焼速度の算出方法
コンベンショナルエンジンに適用している先行点火時期制御方法では、燃焼速度Ｓｂを燃焼室内ガスの層流状態での燃焼速度である層流燃焼速度ＳＬと、燃焼室内ガスの乱流状態での燃焼速度である乱流燃焼速度ＳＴの和、つまり次の（補１）式により燃焼速度Ｓｂを求めている。
【００３２】
Ｓｂ＝ＳＬ＋ＳＴ                                        …（補１）
（補１）式の層流燃焼速度ＳＬ、乱流燃焼速度ＳＴはそれぞれ次の（補２）式、（補３）式により与えている。
【００３３】
ＳＬ＝ＳＬ０×（Ｔ／２９８）＾ｎ×（ｐ／１０１．３２５）＾ｄ
×（１−ｂ×ＭＲＥＳＲ＾ｋ）              …（補２）
ＳＴ＝Ｓｔ・Ｎｅ                                        …（補３）
ただし、ＳＬ０        ：標準状態での層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ            ：未燃ガス温度［Ｋ］、
ｐ            ：燃焼室内圧力［ｋＰａ］、
ＭＲＥＳＲ    ：内部不活性ガス率［％］、
Ｓｔ          ：乱流燃焼速度係数、
ｎ、ｄ、ｂ、ｋ：係数、
（補３）式の乱流燃焼速度ＳＴは、燃焼室５内におけるガスの乱れ強さによって決まる。また、（補３）式のように、乱流燃焼速度Ｓｔをエンジン回転速度Ｎｅの関数としたのは、燃焼室５内におけるガスの乱れ強さは吸気弁１５のバルブリフト量、開閉タイミングが同じ場合、ピストンスピード（エンジン回転速度）に比例するため、乱流燃焼速度Ｓｔの推定にはエンジン回転速度Ｎｅの感度を持たせれば十分であると考えたものである。
【００３４】
しかしながら、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンでは、コンベンショナルエンジンと相違して、吸気弁１５のバルブリフト量や開閉タイミングが運転条件によって大きく変わり、それらバルブリフト量や開閉タイミングの変化による燃焼室５内におけるガスの乱れ強さの変化を無視できないので、本実施形態では、上記（補３）式に代え、次のようにして乱流燃焼速度を算出する。すなわち、ＶＥＬ機構２８の作動に伴う吸気弁１５のバルブリフト量の変化と、ＶＴＣ機構２７の作動に伴う吸気弁閉時期の変化とに対応するため、特にタンブルの影響を燃焼速度の推定式に加える。
【００３５】
まず、燃焼速度Ｓｂ［ｍ／ｓｅｃ］は次の（１）式により表される。（１）式右辺第２項は、乱流燃焼速度ＳＴを与える式で、このように乱流燃焼速度ＳＴを、ｕ／ＳＬを変数としてこの変数の指数関数で与えるとする点は公知である。
【００３６】
Ｓｂ＝ＳＬ＋ｂ・（ｕ／ＳＬ）＾ａ                        …（１）
ただし、ＳＬ    ：層流燃焼速度、
ｕ    ：乱れ強さ、
ａ，ｂ：適合係数、
（１）式の、燃焼室内のガス流動に伴う乱れのうちの１つの特性値である乱れ強さｕ［ｍ／ｓｅｃ］は次の（２）式により表わされる。これは、燃焼室内ガスの乱れ強さｕは回転速度Ｎｅと比例の関係を持つことが、また、その比例定数は、スワール強度、タンブル強度と相関があることが知られているので、これを式に表したものである。
【００３７】
ｕ＝ｃ・Ｉｔ・Ｎｅ                                      …（２）
ただし、Ｉｔ：タンブル強度［無名数］
Ｎｅ：エンジン回転速度［ｒｐｍ］、
ｃ  ：適合係数、
ここで、本実施形態では（２）式にスワール強度は入れていない。これは、今回対象としているエンジンにおいては、ピストン６の中心軸に対して周方向に旋回する流れであるスワールは、ピストン６の中心軸に直交する軸に対して周方向に旋回する流れであるタンブルに比べ、乱流燃焼速度に与える影響が少ないと判断し、今回は考慮しないためである。
【００３８】
ただし、本発明ではスワール強度は対象外というのではない。対象とするエンジンにおいてタンブルよりスワールのほうが乱流燃焼速度に与える影響が大きいときには、スワール強度Ｉｓを、後述するタンブル強度Ｉｔの推定方法と同様の推定方法により推定すればよい。
【００３９】
次に、（２）式のタンブル強度Ｉｔ（燃焼室内のガス流動）の推定方法を説明する。まず、図２（ａ）〜図２（ｅ）は燃焼室５内におけるタンブルの生成から消滅までをピストン６の動きに合わせて図解したものである。順に説明すると、図２（ａ）は吸気弁開時期ＩＶＯ、つまり吸気弁１５が開いた瞬間で、燃焼室５内においてさまざまな方向にガスが流れ込んでいる。図２（ｂ）はピストン６が下降している吸気上死点後４５ｄｅｇ付近（４５ｄｅｇＡＴＤＣ）を示し、この４５ｄｅｇＡＴＤＣ付近からタンブルが形成され始める。図２（ｃ）は吸気弁閉時期ＩＶＣ付近を示し、吸気弁閉時期ＩＶＣ付近までタンブル強度が増加していく。したがって、吸気弁閉時期ＩＶＣが遅いほどタンブルが強くなる。
【００４０】
図２（ｄ）は吸気上死点後３１５ｄｅｇ付近まで、つまりピストン６が上昇する圧縮行程中を示し、圧縮工程中に徐徐にタンブルが小渦（乱れ）に変換されていく。図２（ｅ）は吸気上死点後３１５ｄｅｇ付近を示し、この３１５ｄｅｇＡＴＤＣ付近でタンブルは消え全て小渦に変換される。
【００４１】
一方、図２（ｆ）と図２（ｇ）はＶＥＬ機構２８の作動、非作動に伴う吸気弁１５のバルブリフト量の相違でタンブルの形成がどのように違うのかを示している。このうち、左側に示す図２（ｆ）は吸気弁１５のバルブリフト量が小さい場合で、バルブリフト量が小さいときには、吸気がシリンダ壁側とシリンダ中心側の両方に分かれて流れ、お互いに流れを打ち消しあうため、タンブルが形成されにくい。これに対して右側に示す図２（ｇ）は吸気弁１５のバルブリフト量が大きい場合で、バルブリフト量が大きいときには、ほとんどの吸気がシリンダ中心側へと流れるため、タンブルを形成しやすい。
【００４２】
さて、タンブル強度Ｉｔ［無名数］は次式により小渦の有する角運動量の総和として算出されることがＳＡＥペーパーに記載されている（ＳＡＥ９８１０４８）。
【００４３】
Ｉｔ＝ΣＭiｒi×ｖi                                      …（３）
ただし、ｍi：ｉ番目の小渦におけるガスの質量、
ｒi：シリンダの中心からｉ番目の小渦までの距離、
ｖi：ｉ番目の小渦でのガス流速、
しかしながら、燃焼室内の全ての小渦についてこれらの値ｍi、ｒi、ｖiをエンジンコントローラ３１においてオンラインで算出することは現在のところ不可能であるので、上記（３）式に代えて次の（４）式でタンブル強度Ｉｔ［無名数］を近似する。
【００４４】
Ｉｔ＝ａ・（ｍ＾ｂ・Ｎｅ＾ｃ）                          …（４）
ただし、ｍ    ：燃焼室内ガス質量［ｇ］、
Ｎｅ  ：エンジン回転速度［ｒｐｍ］、
ａ    ：係数、
ｂ、ｃ：適合係数、
（４）式のｍ＾ｂはタンブル強度に対する燃焼室内ガス質量ｍの、またＮｅ＾ｃはタンブル強度に対する回転速度Ｎｅの各影響を反映させたものである。ここで、燃焼室内ガス質量ｍは吸気弁閉時期ＩＶＣに燃焼室５に入っているガス質量のことである。
【００４５】
さらに、タンブルは上記図２（ａ）〜図２（ｅ）で示したように、吸気弁閉時期ＩＶＣが４５ｄｅｇＡＴＤＣ付近より遅いほど強くなり、また、図２（ｆ）、図２（ｇ）で示したように吸気弁１５のバルブリフト量が大きいほうが強くなると考え、（４）式の係数ａを新たに次の（５）式で表す。
【００４６】
ａ＝ｋＶＥＬ・（ＩＶＣ−θ０ｔ）                      …（５）
ただし、ｋＶＥＬ：タンブル強度のバルブリフト量補正係数［１／ｄｅｇ］、
ＩＶＣ  ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θ０ｔ  ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
従って、（５）式を（４）式に代入した次の（６）式によりタンブル強度Ｉｔを近似する。
【００４７】
Ｉｔ＝ｋＶＥＬ・（ＩＶＣ−θ０ｔ）・ｍ＾ｂ・Ｎｅ＾ｃ    …（６）
ただし、ｋＶＥＬ：タンブル強度のバルブリフト量補正係数［１／ｄｅｇ］、
ＩＶＣ  ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θ０ｔ  ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｂ、ｃ  ：適合係数、
（６）式の吸気弁閉時期ＩＶＣの単位は吸気上死点を起点として遅角側に計測する値［ｄｅｇＡＴＤＣ］であるので、吸気弁閉時期ＩＶＣが吸気上死点より遅れるほどＩＶＣの値が大きくなり、従って（６）式よりタンブル強度Ｉｔが大きくなる。また、（６）式のバルブリフト量補正係数ｋＶＥＬはバルブリフト量の関数であり、バルブリフト量が大きいときには、バルブリフト量が小さいときより大きくなり、従って（６）式よりタンブル強度Ｉｔが大きくなる。
【００４８】
（６）式によりタンブル強度Ｉｔを近似した新しい考え方が妥当かどうかを確認するため、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンについて、ＶＥＬ機構２８は非作動状態つまり吸気弁１５のバルブリフト量を一定として、ＶＴＣ機構２７のみを作動し吸気弁閉時期ＩＶＣを変化させる実験を行ったところ、図３に示すようにタンブル強度と吸気弁閉時期ＩＶＣとの関係を示す結果が得られた。図３より吸気弁１５のバルブリフト量一定の条件で、吸気弁閉時期ＩＶＣが大きくなる（つまりＩＶＣが吸気上死点より遅れる）ほどタンブル強度Ｉｔが強くなることが確かめられた。
【００４９】
上記（６）式のタンブル形成開始角θ０ｔは一定値（例えば４５ｄｅｇＡＴＤＣ）である。タンブル形成開始角θ０ｔはエンジン仕様に依存し、ＶＥＬ機構２８、ＶＴＣ機構２７の有無には依存しない。また、図２より、吸気弁閉時期ＩＶＣはタンブル形成開始角θ０ｔより必ず遅角側の値である。従って、上記（６）式の（ＩＶＣ−θ０ｔ）の値は必ず正の値になると考えている。
【００５０】
上記（６）式のタンブル強度のバルブリフト量補正係数ｋＶＥＬはバルブリフト量Ｌｉｆｔ［ｍ］（あるいはバルブ作動角）に応じて定めている。ここでは簡単のため、例えば、ＶＥＬ機構２８の作動、非作動により吸気弁１５のバルブリフト量が大小の２段に切換可能であり、通常の運転条件（運転条件１とする）ではＶＥＬ機構２８を非作動状態とし、運転条件１より所定の運転条件（運転条件２とする）に移行したときにはＶＥＬ機構２８を作動させ、図４下段に示したように実線のバルブリフト特性から破線のバルブリフト特性へとバルブリフト量を小さくしているとして説明すると、ＶＥＬ機構２８の非作動状態でバルブリフト量が第１バルブリフト量Ｌｉｆｔ１［ｍ］となり、ＶＥＬ機構２８の作動状態でバルブリフト量が第２バルブリフト量Ｌｉｆｔ２［ｍ］となる（Ｌｉｆｔ２＜Ｌｉｆｔ１）。従って、図４上段に示したように、ＶＥＬ機構２８の非作動時（あるいは運転条件１のとき）には第１バルブリフト量Ｌｉｆｔ１に対応する第１バルブリフト量補正係数ｋＶＥＬ１を、これに対してＶＥＬ機構２８の作動時（あるいは運転条件２のとき）には第２バルブリフト量Ｌｉｆｔ２に対応する第２バルブリフト量補正係数ｋＶＥＬ２をバルブリフト量補正係数ｋＶＥＬとして設定する。このように、ＶＥＬ機構２８の作動で吸気弁１５のバルブリフト量が第１バルブリフト量Ｌｉｆｔ１より第２バルブリフト量Ｌｉｆｔ２へと小さくなったときに、バルブリフト量補正係数ｋＶＥＬを第１バルブリフト量補正係数ｋＶＥＬ１より第２バルブリフト量補正係数ｋＶＥＬ２へと小さくするのは、バルブリフト量が小さいときのほうがタンブルが弱くなる（従って、タンブル強度Ｉｔが小さくなる）からである。
【００５１】
このように、（６）式によれば、タンブル強度Ｉｔを、新たに、吸気弁１５のバルブリフト量、吸気弁閉時期ＩＶＣ、燃焼室内ガス質量ｍ、エンジン回転速度Ｎｅの関数としている。
〈２〉基準クランク角の設定方法
コンベンショナルエンジンに適用している先行点火時期制御方法では、次の（補４）式のように、燃焼圧が最大となるクランク角θｐｍａｘを、燃焼質量割合が６０％となるときのクランク角位置に設定すると共に、この燃焼圧が最大となるクランク角θｐｍａｘを基準クランク角としている。なお、燃焼質量割合がｘ％となるときのクランク角位置を、以下、「θｍｂ ｘ％」で表す。例えば、燃焼質量割合が６０％となるときのクランク角位置は「θｍｂ６０％」である。また、θｍｂ ｘ％の起点は圧縮上死点とする。
【００５２】
θｍｂ６０％＝θｐｍａｘ                                …（補４）
また、燃焼圧が最大となるクランク角θｐｍａｘを次の（補５）式のようにエンジン回転速度Ｎｅの関数で与えている。
【００５３】
θｐｍａｘ＝ｈ・（Ｎｅ）＾ｉ                            …（補５）
ただし、ｈ、ｉ：適合係数、
一方、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象として実験してみたところ、図５に示したように、回転速度Ｎｅが一定の条件のもとで、燃焼圧が最大となるクランク角θｐｍａｘとθｍｂ ｘ％との関係を表す実験結果が得られた。ただし、図５にはθｍｂ ｘ％として燃焼質量割合が０％となるときのクランク角位置であるθｍｂ０％、燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％、燃焼質量割合が６０％となるときのクランク角位置であるθｍｂ６０％の３つの場合だけを示しており、図５において左はθｐｍａｘとθｍｂ０％の、中央はθｐｍａｘとθｍｂ１０％の、右はθｐｍａｘとθｍｂ６０％の各関係を整理したものである。実際には燃焼質量割合が０％、１０％、６０％以外の値（２０％、３０％、４０％、５０％）となるときのクランク角位置での実験結果も得ている（図７参照）。
【００５４】
燃焼圧が最大となるクランク角θｐｍａｘと、燃焼質量割合が６０％となるときのクランク角位置であるθｍｂ６０％とでは、図５右側に示したように最大誤差が２．６ｄｅｇもあり、θｐｍａｘとθｍｂ０％との間の誤差を示す図５左側と同様であることがわかる。
【００５５】
次に、図６は、吸入空気量（エンジン負荷）が一定の条件のもとでの、燃焼圧が最大となるクランク角θｐｍａｘとエンジン回転速度Ｎｅとの関係を示したものである。図示のように、傾向としては回転速度Ｎｅが大きくなるほどθｐｍａｘが進角側にずれているが、回転速度Ｎｅに対するθｐｍａｘのバラツキが大きく、燃焼圧が最大となるクランク角θｐｍａｘは上記（補５）式のようにエンジン回転速度Ｎｅだけでは表しきれないことがわかる。
【００５６】
図７は燃焼質量割合を横軸に、燃焼質量割合がｘ％となるときのクランク角位置であるθｍｂ ｘ％をθｐｍａｘで直線近似した場合のＲの２乗値を縦軸に採り、図５に示したθｍｂ０％時、θｍｂ１０％時、θｍｂ６０％時の３つの実験結果及び図示しない図５と同様のθｍｂ２０％時、θｍｂ３０％時、θｍｂ４０％時、θｍｂ５０％時の残りの実験結果とを整理したものである。縦軸のＲの２乗値は、この値が１．０に近づくほどθｐｍａｘとθｍｂ ｘ％（ｘ＝０、１０、２０、３０、４０、５０、６０）との間に相関があることを表すので、図７によれば、燃焼圧が最大となるクランク角θｐｍａｘはθｍｂ１０％付近で一番相関があることがわかる。
【００５７】
なお、図７にはＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンでもエンジン仕様の異なる２種類のエンジンの特性を示している。２つのエンジンの違いは主に燃焼室形状にあり、第２エンジンは第１エンジンよりも燃焼室が平べったい形状のものである。第１エンジンでは燃焼圧が最大となるクランク角θｐｍａｘは燃焼質量割合が１０％となるときのクランク角位置から５０％となるときのクランク角位置との間で広く相関があるのに対して、第２エンジンになると燃焼圧が最大となるクランク角θｐｍａｘは燃焼質量割合が１０％となるときのクランク角位置付近でだけ相関がある。従って、２つのエンジンに共通して燃焼圧が最大となるクランク角θｐｍａｘが高い相関を有するのは燃焼質量割合が１０％となるときのクランク角位置付近にあるときだけである。このように、エンジン機種が相違しても、燃焼圧が最大となるクランク角θｐｍａｘはθｍｂ１０％となるときのクランク角位置付近で一番相関がある、という結果が初めて得られた。
【００５８】
こうした実験結果を受けて本実施形態では、燃焼圧が最大となるクランク角θｐｍａｘに基づいて、燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％を算出し、このθｍｂ１０％を基準クランク角として設定する、つまり次の（７）式によりθｍｂ１０％［ｄｅｇＡＴＤＣ］を算出する。
【００５９】
θｍｂ１０％＝ａ・θｐｍａｘ−ｂ                        …（７）
ただし、ａ：適合係数［無名数］、
ｂ：適合クランク角［ｄｅｇ］、
（７）式は燃焼圧が最大となるクランク角θｐｍａｘをｂ（正の値）のクランク角だけ進角側にシフトしてθｍｂ１０％を求めるとする式である。この結果、図８に示したように本実施形態での基準クランク角（＝θｍｂ１０％）はコンベンショナルエンジンに適用している先行点火時期制御方法における基準クランク角（＝θｍｂ６０％）よりも進角側にくることとなる。
【００６０】
そして、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンでは、燃焼圧が最大となるクランク角θｐｍａｘは、吸気弁閉時期ＩＶＣが遅いほど、また吸気弁１５のバルブリフト量が大きいほど圧縮上死点より遅れると考え、また、燃焼室内ガス質量（負荷）によって影響されると判断し、さらにＦＴＡ（Ｆａｕｌｔ  Ｔｒｅｅ  Ａｎａｌｙｓｉｓ）より燃焼室内のガス流速（タンブル、スワール）の影響をもつと判断し、次の（８）式により燃焼圧が最大となるクランク角θｐｍａｘ［ｄｅｇＡＴＤＣ］を算出することとした。従って、（８）式と上記（補５）式とを比較すれば、ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｄを新たに導入したものである。
【００６１】
θｐｍａｘ＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｄ・Ｎｅ＾ｅ
…（８）
ただし、ｋ１ＶＥＬ：流速係数［無名数］、
ＩＶＣ    ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θｔ０    ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｍ        ：燃焼室内ガス質量［ｇ］、
Ｎｅ      ：エンジン回転速度［ｒｐｍ］、
ｄ、ｅ    ：適合定数［無名数］、
（８）式に示す燃焼圧が最大となるクランク角θｐｍａｘに対するこうした新しい考え方が妥当かどうかを確認するため、吸気弁１５のバルブリフト量及び吸気弁閉時期ＩＶＣを一定として、燃焼圧が最大となるクランク角θｐｍａｘと吸入空気量（燃焼室内ガス質量ｍの代表値）、エンジン回転速度Ｎｅ（ガス流速の代表値）の関係を実験したところ、図９に示す結果が得られた。図９によれば、吸気弁１５のバルブリフト量及び吸気弁閉時期ＩＶＣが一定の条件で、吸入空気量が同じであればエンジン回転速度Ｎｅが大きいときのほうがエンジン回転速度Ｎｅが小さいときより燃焼圧が最大となるクランク角θｐｍａｘが大きくなる（圧縮上死点よりの遅角量が大きくなる）こと、またエンジン回転速度Ｎｅが同じであれば吸入空気量が多くなるほど燃焼圧が最大となるクランク角θｐｍａｘが大きくなる（圧縮上死点よりの遅角量が大きくなる）ことが確かめられた。
【００６２】
なお、図９の縦軸の単位［ｄｅｇＡＴＤＣ］は圧縮上死点を起点とするクランク角であり、吸気上死点を起点とするものでない。
【００６３】
上記（８）式のタンブル形成開始角θ０ｔについては上記（６）式のところで前述したところと同じである。すなわち、タンブル形成開始角θ０ｔは一定値（例えば４５ｄｅｇＡＴＤＣ）である。タンブル形成開始角θ０ｔはエンジンの仕様に依存し、ＶＥＬ機構２８、ＶＴＣ機構２７の有無には依存しない。図２より、吸気弁閉時期ＩＶＣはタンブル形成開始角θ０ｔより必ず遅角側の値である。従って、（８）式の（ＩＶＣ−θ０ｔ）の値は必ず正の値になる。
【００６４】
上記（８）式の流速係数ｋ１ＶＥＬは吸気弁１５のバルブリフト量Ｌｉｆｔ（あるいはバルブ作動角）に応じて定めている。ここでも簡単のため、前述したようにＶＥＬ機構２８の作動、非作動により吸気弁１５のバルブリフト量が大小の２段に切換可能であり、運転条件１ではＶＥＬ機構２８を非作動状態とし、運転条件１より運転条件２に移行したときにはＶＥＬ機構２８を作動させ、図１０下段に示したように実線のバルブリフト特性から破線のバルブリフト特性へとバルブリフト量を小さくしているとして説明すると、図１０上段に示したように、ＶＥＬ機構２８の非作動時（あるいは運転条件１のとき）には第１バルブリフト量Ｌｉｆｔ１に対応する第１流速係数ｋ１ＶＥＬ１を、これに対してＶＥＬ機構２８の作動時（あるいは運転条件２のとき）には第２バルブリフト量Ｌｉｆｔ２に対応する第２流速係数ｋ１ＶＥＬ２を流速係数ｋ１ＶＥＬとして設定する。このように、バルブリフト量が小さくなったときに流速係数ｋ１ＶＥＬを第１流速係数ｋ１ＶＥＬ１より第２流速係数ｋ１ＶＥＬ２へと小さくするのは、バルブリフト量が小さいときのほうがタンブルが弱くなる（従って、θｐｍａｘがより遅角側に移動する）からである。
【００６５】
また、上記（８）式のｍ＾ｄはθｐｍａｘに対する筒内ガス質量ｍの、Ｎｅ＾ｅはθｐｍａｘに対する回転速度Ｎｅの各影響を反映させている。
【００６６】
このように、（８）式によれば、燃焼圧が最大となるクランク角θｐｍａｘを、エンジン回転速度Ｎｅに加えて、燃焼室内ガス質量ｍ、吸気弁１５のバルブリフト量及び吸気弁閉時期ＩＶＣの関数でもあるとして新たに構成している。
〈３〉着火遅れ時間の算出方法
コンベンショナルエンジンに適用している先行点火時期制御方法では、点火タイミングを起点として燃焼質量割合が０％となるときのクランク角位置θｍｂ０％である間の区間を着火遅れ時間ＤＥＡＤＴＩＭＥとして算出している。
【００６７】
しかしながら、燃焼質量割合が０％となるときのクランク角位置であるθｍｂ０％は冷却損失の算出誤差やガス漏れの影響を受けやすく、正確に計測することが難しい。そのため、このθｍｂ０％の値を正として着火遅れ時間ＤＥＡＤＴＩＭＥの適合を行うと燃焼速度の推定誤差に着火遅れ時間の計測誤差が含まれてしまう。また、このような理由により、一般的にもθｍｂ０％を正にして燃焼シミュレーションの適合を行うこともまれである。
【００６８】
これに対して本実施形態では、点火タイミングを起点として燃焼質量割合が２％となるときのクランク角位置であるθｍｂ２％までの時間を着火遅れ時間τとみなす。このため、本実施形態ではθｍｂ２％よりθｍｂ１０％までのクランク角区間が燃焼期間となるが、これはコンベンショナルエンジンに適用している先行点火時期制御方法では主燃焼期間ＢＵＲＮ２に相当しているので、本実施形態では、燃焼期間の推定を主燃焼期間のみとすることとなる。
【００６９】
また、この着火遅れ時間τは純粋な意味での着火遅れ時間ではなく、点火タイミングから燃料が２％燃えるまでの期間である。そのため、コンベンショナルエンジンに適用している先行点火時期制御方法において提案しているアレニウスの式（自己着火の式）を用いることができないため、本実施形態では、着火遅れ時間τ［ｓｅｃ］を新たに次の（９）式により算出することとする。この（９）式は点火直後の燃焼速度を算出する式を用いて算出するものである。
【００７０】
τ＝（Ｄkernel−Ｄ０）
／２［（Ｔａｄ／Ｔ）・ＳＬ＋｛（２／３）・ｋ｝＾（１／２）］
…（９）
ただし、Ｄkernel：θｍｂ２％時の火炎直径［ｍ］、
Ｄ０    ：点火タイミング直後の火炎直径［ｍ］、
Ｔａｄ  ：火炎温度［Ｋ］、
Ｔ      ：未燃ガス温度［Ｋ］、
ＳＬ    ：層流燃焼速度［ｍ／ｓｅｃ］、
ｋ      ：燃焼室内ガスの運動エネルギー、
（９）式の火炎直径Ｄkernelは燃焼質量割合が２％となったときの火炎の直径であるので、次の（１０）式により新たに算出する。
【００７１】
Ｄkernel＝｛（Ｖｃｙｌ・［１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）］
・６）／π｝＾（１／３）              …（１０）
ただし、Ｖｃｙｌ：燃焼室内容積［ｍ3］、
ｘｂ    ：燃焼質量割合（＝０．０２）、
κ      ：ポリトロープ指数（＝１．３４）、
（１０）式は次のようにして導いたものである。火炎の形状を球と仮定すると、火炎直径Ｄkernelを与える式は次の（１１ａ）式である。
【００７２】
Ｄkernel＝｛（Ｖｂ・６）／π｝＾（１／３）              …（１１ａ）
ただし、Ｖｂ：火炎体積、
熱力学の公式より体積燃焼割合の式は次の（１１ｂ）式で表される（公知）。
【００７３】
体積燃焼割合＝１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）
…（１１ｂ）
この（１１ｂ）式と、
Ｖｂ＝Ｖｃｙｌ・体積燃焼割合                            …（１１ｃ）
の関係式とを利用すれば次の（１１ｄ）式が得られる。
【００７４】
Ｖｂ＝Ｖｃｙｌ・［１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）］
…（１１ｄ）
この（１１ｄ）式を（１１ａ）式に代入すれば上記（１０）式が得られる。
【００７５】
上記（９）式の点火タイミング直後の火炎の直径Ｄ０は１ｍｍとする。火炎温度Ｔａｄと未燃ガス温度Ｔの比であるＴａｄ／Ｔは適合定数とする。層流燃焼速度ＳＬの算出方法はコンベンショナルエンジンに適用している先行点火時期制御方法と同じとする。
【００７６】
上記（９）式の燃焼室内ガスの運動エネルギーｋは次式により算出する。
【００７７】
ｋ＝（１／２）ｍ（ｆ・Ｎｅ）＾２                        …（１２）
ただし、ｍ  ：燃焼室内ガス質量［ｇ］、
Ｎｅ：エンジン回転速度［ｒｐｍ］、
ｆ  ：適合定数、
これで、項分け説明を終了する。
【００７８】
次に、上記〈１〉〜〈３〉の変更を行った後の新しい点火時期制御方法（この新しい点火時期制御方法を以下「変更後点火時期制御方法」という。）を、以下にまとめる。ただし、適合係数については改めて取り直している。
【００７９】
ＭＢＴの得られる点火時期（この点火時期を「基本点火時期」という。）ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を次式により算出する。
【００８０】
ＭＢＴＣＡＬ＝−｜（τ＋ＢＴ）・Ｎｅ・６−θｍｂ１０％｜…（１３）
ただし、ＢＴ        ：燃焼時間［ｓｅｃ］、
τ          ：着火遅れ時間［ｓｅｃ］、
θｍｂ１０％：基準クランク角［ｄｅｇＡＴＤＣ］、
これを図解したのが図１１である。変更後火時期制御方法では、点火タイミングからθｍｂ２％までの時間である着火遅れ時間τ［ｓｅｃ］と、θｍｂ２％からθｍｂ１０％までの時間である燃焼時間ＢＴ［ｓｅｃ］とを加算し、この加算値［ｓｅｃ］にＮｅ［ｒｐｍ］・６を乗算することによってクランク角区間［ｄｅｇ］に変換する。この換算されたクランク角区間（（τ＋ＢＴ）・Ｎｅ・６）は点火タイミングからθｍｂ１０％（＝θｐｍａｘ）までのクランク角区間である。従って、基準クランク角（＝θｐｍａｘ）であるθｍｂ１０％［ｄｅｇＡＴＤＣ］よりこの換算されたクランク角区間だけ進角側の値を基本点火時期ＭＢＴＣＡＬとして算出する。
【００８１】
なお、（１３）式において絶対値をとりマイナスの符号を付けているのは、θｍｂ１０％の単位は圧縮上死点より遅角側に計測するクランク角［ｄｅｇＡＴＤＣ］であるのに対して、ＭＢＴＣＡＬの単位は圧縮上死点より進角側に計測するクランク角［ｄｅｇＢＴＤＣ］であるため、θｍｂ１０％の単位をＭＢＴＣＡＬの単位へと変換してやる必要があるためである。
【００８２】
ここで、基準クランク角であるθｍｂ１０％［ｄｅｇＡＴＤＣ］は、次の（１４）式により算出する。
【００８３】
θｍｂ１０％＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｃ１・Ｎｅ＾ｃ２＋ｃ３
…（１４）
ただし、ｋ１ＶＥＬ：流速係数［無名数］、
ＩＶＣ    ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θｔ０    ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｍ        ：燃焼室内ガス質量［ｇ］、
Ｎｅ      ：エンジン回転速度［ｒｐｍ］、
ｃ１〜ｃ３：適合係数、
前述したように、吸気弁閉時期ＩＶＣとタンブル形成開始角θｔ０の単位は吸気上死点を起点としており、一方、θｍｂ１０％は圧縮上死点を起点としている。従って、吸気上死点を起点とする単位より圧縮上死点を起点とする単位への変換は適合係数ｃ３により行うこととなる。例えば４気筒エンジンでは吸気上死点と圧縮上死点の間に１８０ｄｅｇのずれがあるので、適合係数ｃ３に１８０ｄｅｇを入れてやればよい。
【００８４】
（１３）式の燃焼時間ＢＴ［ｓｅｃ］は次式により算出する。
【００８５】
ＢＴ＝０．０８・Ｖｃｙｌ／｛（（３・ｘｂ＋１）＾（１／κ））・Ａｂ・Ｓｂ｝
＝０．０７０７・Ｖｃｙｌ／（Ａｂ・Ｓｂ）            …（１５）
ただし、Ｖｃｙｌ：燃焼室内容積［ｍ3］、
ｘｂ    ：燃焼時間算出時の平均燃焼質量割合（＝０．０６）、
κ      ：ポリトロープ指数、
Ａｂ    ：火炎表面積［ｍ2］、
Ｓｂ    ：燃焼速度［ｍ／ｓｅｃ］、
（１５）式はコンベンショナルエンジンに適用している先行点火時期制御方法における式と基本的に同じ式である。すなわち、（１５）式においてＶｃｙｌ／｛（３・ｘｂ＋１）＾（１／κ）｝は既燃ガスの質量割合であり、燃焼時間ＢＴはこの既燃ガスの質量割合に比例し、燃焼速度Ｓｂに反比例するとする式である。（１５）式の火炎表面積Ａｂはエンジンの仕様から算出することができる。
【００８６】
なお、（１５）式の燃焼質量割合ｘｂとしては、θｍｂ２％からθｍｂ１０％までのクランク角区間（燃焼期間）における値を用いる必要がある。この場合に、燃焼質量割合はθｍｂ２％からθｍｂ１０％までのクランク角区間で２％から１０％へと変化するので、ここでは、θｍｂ２％とθｍｂ１０％とを平均した値であるθｍｂ６％時の燃焼質量割合の値つまり６％を用いる。
【００８７】
上記（１５）式の燃焼速度Ｓｂ［ｍ／ｓｅｃ］は次式により算出する。
【００８８】
Ｓｂ＝ＳＬ＋［｛ｋＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｃ４・Ｎｅ＾ｃ５｝
／ＳＬ］＾ｃ６                                    …（１６）
ただし、ＳＬ      ：層流燃焼速度［ｍ／ｓｅｃ］、
ｋＶＥＬ  ：タンブル強度のバルブリフト量補正係数［１／ｄｅｇ］、
ＩＶＣ    ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θ０ｔ    ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｃ４、ｃ５：適合係数、
ｃ６      ：適合係数、
（１６）式右辺第２項が乱流燃焼速度［ｍ／ｓｅｃ］を与える式であり、乱流燃焼速度を（１６）式右辺第２項で近似している点が新しい。つまり、本実施形態では、燃焼室内ガスの乱れ強さｕ［ｍ／ｓｅｃ］を次の（１６−１）式により算出し、この燃焼室内ガスの乱れ強さｕに基づいて次の（１６−２）式により乱流燃焼速度ＳＴ［ｍ／ｓｅｃ］を算出している。
【００８９】
ｕ＝ｋＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｃ４・Ｎｅ＾ｃ５  …（１６−１）
ＳＴ＝［ｕ／ＳＬ］＾ｃ６                                …（１６−２）
これに対して、（１６）式の層流燃焼速度ＳＬ［ｍ／ｓｅｃ］としては、コンベンショナルエンジンに適用している先行点火時期制御方法と同じでよく、従って次の（１７）式により算出する。
【００９０】
ＳＬ＝ＳＬ０・（Ｔ／２９８）＾ａ・（ｐ／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ＾ｃ）            …（１７）
ただし、ＳＬ０    ：標準状態での層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ        ：未燃ガス温度［Ｋ］、
ｐ        ：燃焼室内圧力［ｋＰａ］、
ＭＲＥＳＲ：内部不活性ガス率［％］、
ａ、ｂ、ｃ：係数、
一方、上記（１３）式の着火遅れ時間τ［ｓｅｃ］は次の（１８）式により新たに算出する。
【００９１】
τ＝［｛｛Ｖｃｙｌ・［１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）］・６｝            ／π｝＾（１／３）−０．００１］
／｛ｃ７／ＳＬ＋ｃ８・Ｎｅ・ｍ＾（１／２）｝
＝｛０．４９０・Ｖｃｙｌ＾（１／３）｝
／｛ｃ７・ＳＬ＋ｃ８・Ｎｅ・ｍ＾（１／２）｝    …（１８）
ただし、ｘｂ      ：着火遅れ時間算出時の燃焼質量割合（＝０．０２）、
κ        ：ポリトロープ指数（＝１．３４）、
ＳＬ      ：層流燃焼速度［ｍ／ｓｅｃ］、
ｍ        ：燃焼室内ガス質量［ｇ］、
Ｎｅ      ：エンジン回転速度［ｒｐｍ］、
Ｖｃｙｌ  ：燃焼室内容積［ｍ3］、
ｃ７、ｃ８：適合定数、
なお、（１８）式の燃焼質量割合ｘｂとしては、θｍｂ０％からθｍｂ２％までのクランク角区間（着火遅れ期間）における値を用いる必要がある。この場合に、燃焼質量割合はθｍｂ０％からθｍｂ２％までのクランク角区間で０％より２％まで変化するのであるが、ここでは、θｍｂ０％とθｍｂ２％とを平均した値であるθｍｂ１％時の燃焼質量割合の値を用いるのではなく、着火遅れ期間の終期であるθｍｂ２％時の燃焼質量割合の値つまり２％を用いることとする。
【００９２】
上記（１８）式の層流燃焼速度ＳＬ［ｍ／ｓｅｃ］としては、次の（１９）式により算出する。
【００９３】
ＳＬ＝ＳＬ０・（Ｔ／２９８）＾ａ・（ｐ／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ＾ｃ）            …（１９）
ただし、ＳＬ０    ：標準状態での層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ        ：未燃ガス温度［Ｋ］、
ｐ        ：燃焼室内圧力［ｋＰａ］、
ＭＲＥＳＲ：内部不活性ガス率［％］、
ａ、ｂ、ｃ：係数、
ここで、着火遅れ時間τを算出するのに用いる（１９）式の層流燃焼速度ＳＬと、燃焼時間ＢＴを算出するのに用いる上記（１７）式の層流燃焼速度ＳＬとでは、式そのものは変わらず、また標準状態での層流燃焼速度ＳＬ０、内部不活性ガス率ＭＲＥＳＲの値も変わらないが、未燃ガス温度Ｔ、燃焼室内圧力ｐとして代入する値が（１９）式と（１７）式とで、後述するように相違することとなる。
【００９４】
上記（１７）式及び（１９）式の標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂ、ｃ［無名数］は次の（２０）式〜（２３）式により算出する。
【００９５】
ＳＬ０＝（０．２６３２−０．８４７２／（φ−１．１３）＾２）
…（２０）
ａ＝２．１８−０．８０・（φ−１）                      …（２１）
ｂ＝−０．１６＋０．２２・（φ−１）                    …（２２）
ｃ＝１                                                  …（２３）
ただし、φ：当量比［無名数］、
ここで、（２０）式〜（２３）式についてはＳＡＥペーパーにより公知である（ＳＡＥ１９９９１０１７５参照）。
【００９６】
これで変更後点火時期制御方法をまとめたものの説明を終了する。
【００９７】
次に、上記（１１ｂ）式に示した体積燃焼割合の求め方と、上記（１７）式、（１９）式に示した未燃ガス温度の算出方法とについて補足説明を行う。
【００９８】
まず、上記（１１ｂ）式に示した体積燃焼割合の求め方を簡単に説明する。
【００９９】
既燃ガスの質量燃焼割合をｘ、体積燃焼割合をｙとして、燃焼室内の燃焼が定容燃焼同様に、
ｙ＝ｆ（ｘ）
＝１＋（ｘ−１）／（１＋ｘ（ｋ１−１）＾（１／κ））
…（２４）
ただし、ｋ１≒４〜５、
の式により書けると仮定すると、（２４）式のｋ１に４を代入して計算することにより、上記（１１ｂ）式が容易に得られる。
【０１００】
上記（１７）式、（１９）式に示した未燃ガス温度Ｔの算出方法をまとめて説明すると、断熱変化と仮定し、熱力学の次の（２５）式を用いて燃焼室５のθｍｂ ｘ％時の未燃ガス温度Ｔｍｂ ｘ％を算出する。
【０１０１】
Ｔｍｂ ｘ％＝Ｔｉｖｃ・（Ｐｍｂ ｘ％／Ｐｉｖｃ）＾｛（κ−１）／κ｝                                                                      …（２５）
ただし、Ｔｍｂ  ｘ％：θｍｂ ｘ％時の未燃ガス温度［Ｋ］、
Ｔｉｖｃ    ：ＩＶＣ時の燃焼室内温度［Ｋ］、
Ｐｍｂ  ｘ％：θｍｂ ｘ％時の燃焼室内圧力［ｋＰａ］、
Ｐｉｖｃ    ：ＩＶＣ時の燃焼室内圧力［ｋＰａ］、
κ          ：比熱比（固定値で１．３〜１．４）、
ここで、（２５）式の燃焼室５のθｍｂ ｘ％時の未燃ガス圧力Ｐｍｂ ｘ％の算出方法を説明すると、この燃焼室５のθｍｂ ｘ％時の未燃ガス圧力Ｐｍｂ ｘ％は次の（２６）式により与えられる。
【０１０２】
Ｐｍｂ ｘ％＝Ｐｉｖｃ＋Δｐ                              …（２６）
ただし、Δｐ：燃焼室内未燃ガス圧力の変化量、
ここで、燃焼室内未燃ガス圧力の変化量Δｐは、容積変化による圧力変化分Δｐｖと、燃焼による圧力変化分Δｐｃとに分けて算出する。つまり、次の（２７）式により燃焼室内未燃ガス圧力の変化量Δｐを算出する。
【０１０３】
Δｐ＝Δｐｖ＋Δｐｃ                                    …（２７）
このうち、（２７）式の容積変化による圧力変化分Δｐｖは熱力学の式である次の（２８）式により算出することができる。
【０１０４】
Δｐｖ＝Ｐｉｖｃ・｛（Ｖｉｖｃ／Ｖｍｂ ｘ％）＾κ−１｝ …（２８）
ただし、Ｖｉｖｃ   ：ＩＶＣ時の燃焼室内容積、
Ｖｍｂ ｘ％：θｍｂ ｘ％時の燃焼室内容積、
ここで、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃはピストン６の位置から求めることができる。燃焼室５のθｍｂ ｘ％時における未燃ガス容積Ｖｍｂ  ｘ％については、図１２に示したように、θｍｂ ｘ％とＶｍｂ ｘ％の関係を表す特性を予め作成しておき、θｍｂ ｘ％から当該特性を用いて、燃焼室５のθｍｂ ｘ％時における未燃ガス容積Ｖｍｂ ｘ％を求めればよい。図１２は要するにクランク角に対する燃焼室容積の特性である。
【０１０５】
一方、上記（２７）式の燃焼による圧力変化分Δｐｃは、次の（２９）式により簡易に算出すればよい。
【０１０６】
Δｐｃ＝（ｘ／１００）×Δｐｔｏｔａｌ                  …（２９）
ただし、Δｐｔｏｔａｌ：燃焼によるトータルの圧力上昇分、
ここで、燃焼によるトータルの圧力上昇分Δｐｔｏｔａｌについては、図１３に示したように、燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］（または燃料噴射量）とΔｐｔｏｔａｌの関係を表す特性を予め作成しておき、燃料噴射パルス幅Ｔｉ（または燃料噴射量）から当該特性を用いて燃焼によるトータルの圧力上昇分Δｐｔｏｔａｌを求める。
【０１０７】
これで補足説明を終える。
【０１０８】
次に、エンジンコントローラ３１で実行される上記変更後点火時期制御方法における基本点火時期ＭＢＴＣＡＬの算出方法を図１４、図１６、図１７、図１８、図１９のフローチャートを参照しながら詳述する。
【０１０９】
図１４は基本点火時期ＭＢＴＣＡＬの算出に必要な各種の物理量を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。
【０１１０】
まずステップ１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、バルブリフト量Ｌｉｆｔ［ｍ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、内部不活性ガス率ＭＲＥＳＲ［％］、燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］、クランク角センサ（３３、３４）により検出されるエンジン回転速度Ｎｅ［ｒｐｍ］、総ガス質量ＭＧＡＳ［ｇ］を読み込む。
【０１１１】
ここで、吸気弁閉時期ＩＶＣはＶＴＣ機構２７に与える指令値から既知である。あるいはカム角センサ３４により実際の吸気弁閉時期を検出してもかまわない。
【０１１２】
吸気弁１５のバルブリフト量ＬｉｆｔはＶＥＬ機構２８に与える指令値から既知である。つまり、図４下段、図１０下段に示したように、ＶＥＬ機構２８の非作動時には第１バルブリフト量Ｌｉｆｔ１が、これに対してＶＥＬ機構２８の作動時には第２バルブリフト量Ｌｉｆｔ２がバルブリフト量Ｌｉｆｔである。
【０１１３】
燃料噴射パルス幅Ｔｉは図示しない燃料噴射パルス幅の算出ルーチンにおいて算出されている。例えば、シーケンシャル噴射時に燃料インジェクタ２１に与える燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］の算出式は次のようなものである。
【０１１４】
Ｔｉ＝Ｔｐ×Ｔｆｂｙａ×（α＋αｍ−１）×２＋Ｔｓ      …（３０）
ただし、Ｔｐ      ：基本噴射パルス幅［ｍｓｅｃ］、
Ｔｆｂｙａ：［無名数］、
α        ：空燃比フィードバック補正係数［無名数］、
αｍ      ：空燃比学習値［無名数］、
Ｔｓ      ：無効パルス幅［ｍｓｅｃ］、
（３０）式の目標当量比Ｔｆｂｙａは無名数であり、理論空燃比を１４．７とすると、次の（３１）式により表される値である。
【０１１５】
Ｔｆｂｙａ＝１４．７／目標空燃比                        …（３１）
例えば（３１）式より目標空燃比が理論空燃比（１４．７）のときＴｆｂｙａ＝１．０となり、目標空燃比が例えば２２．０といったリーン側の値であるとき、Ｔｆｂｙａは１．０未満の正の値となる。
【０１１６】
クランク角センサはクランクシャフト７のポジションを検出するポジションセンサ３３と、吸気用カムシャフト２５ポジションを検出するフェーズセンサ（＝カム角センサ）３４とからなり、これら２つのセンサ３３、３４からの信号に基づいてエンジン回転速度Ｎｅ［ｒｐｍ］が算出されている。
【０１１７】
内部不活性ガス率ＭＲＥＳＦＲは燃焼室５内に残留する不活性ガス量を燃焼室５内の総ガス量で除した値で、その算出については総ガス質量ＭＧＡＳと共に図１６により後述する。
【０１１８】
ステップ２では、燃焼室５の吸気弁閉時期ＩＶＣにおける容積（圧縮開始時期の容積）Ｖｉｖｃ［ｍ3］を算出する。燃焼室５の吸気弁閉時期における容積Ｖｉｖｃは、ピストン６のストローク位置によって決まる。ピストン６のストローク位置はエンジンのクランク角位置によって決まる。
【０１１９】
図１５を参照して、エンジンのクランクシャフト７１の回転中心７２がシリンダの中心軸７３からオフセットしている場合を考える。コネクティングロッド７４、コネクティングロッド７４とクランクシャフト７１との結節点７５、コネクティングロッド７４とピストンをつなぐピストンピン７６が図に示す関係にあるとする。このときの、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃは次の（４１）式〜（４５）式で表すことができる。
【０１２０】
Ｖｉｖｃ＝ｆ１（θivc）＝Ｖｃ＋（π／４）Ｄ2・Ｈivc
…（４１）
Ｖｃ＝（π／４）Ｄ2・Ｈx／（ε−１）                  …（４２）
Ｈivc＝｛（ＣＮＤ＋ＳＴ2／２）−（ＣＲoff−ＰＩＳoff）2｝1/2
−｛（ＳＴ／２）・cos（θivc＋θoff）｝
＋（ＣＮＤ2−Ｘ2）1/2                            …（４３）
Ｘ  ＝（ＳＴ／２）・sin（θivc＋θoff）−ＣＲoff＋ＰＩＳoff
…（４４）
θoff＝arcsin｛（ＣＲoff−ＰＩＳoff）／（ＣＮＤ・（ＳＴ／２））｝
…（４５）
ただし、Ｖｃ     ：隙間容積［ｍ3］、
ε       ：圧縮比、
Ｄ       ：シリンダボア径［ｍ］、
ＳＴ     ：ピストンの全ストローク［ｍ］、
Ｈivc    ：吸気弁閉時期におけるピストンピン７６の
ＴＤＣからの距離［ｍ］、
Ｈx      ：ピストンピン７６のＴＤＣからの距離の最大値と最小値の                           差[ｍ］、
ＣＮＤ   ：コネクティングロッド７４の長さ［ｍ］、
ＣＲoff  ：結節点７５のシリンダ中心軸７３からのオフセット距離
［ｍ］、
ＰＩＳoff：クランクシャフト回転中心７２のシリンダ中心軸７３から                             のオフセット距離［ｍ］、
θivc    ：吸気弁閉時期のクランク角［ｄｅｇＡＴＤＣ］、
θoff    ：ピストンピン７６とクランクシャフト回転中心７２とを結                           ぶ線がＴＤＣにおいて垂直線となす角度［ｄｅｇ］、
Ｘ       ：結節点７５とピストンピン７６との水平距離［ｍ］、
吸気弁閉時期のクランク角θivcは前述のように、エンジンコントローラ３１からＶＴＣ機構２７への指令信号によって決まるので、既知である。（４１）式〜（４５）式にこのときのクランク角θivc（＝ＩＶＣ）を代入すれば、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃを算出することができる。したがって、実用上は燃焼室５の吸気弁閉時期における容積Ｖｉｖｃは吸気弁閉時期ＩＶＣをパラメータとするテーブルで設定したものを用いる。
【０１２１】
ステップ３では、燃焼室５の吸気弁閉時期ＩＶＣにおける温度（圧縮開始時期温度）Ｔｉｖｃ［Ｋ］を算出する。燃焼室５に流入するガスの温度は、燃焼室５に流入する新気と燃焼室５に残留する不活性ガスとが混じったガスの温度であり、燃焼室５に流入する新気の温度は吸気コレクタ２内の新気温度ＴＣＯＬに等しく、また燃焼室５内に残留する不活性ガスの温度は排気ポート部近傍の排気温度ＴＥＸＨで近似できるので、燃焼室５の吸気弁閉時期ＩＶＣにおける温度Ｔｉｖｃは吸気弁閉時期ＩＶＣになったタイミングでの、吸気コレクタ２内の新気温度ＴＣＯＬ、排気温度ＴＥＸＨ、燃焼室５内に残留する不活性ガスの割合である内部不活性ガス率ＭＲＥＳＲから次の（４６）式により求めることができる。
【０１２２】
Ｔｉｖｃ＝ＴＥＸＨ×ＭＲＥＳＲ＋ＴＣＯＬ×（１−ＭＲＥＳＲ）
…（４６）
ステップ４では燃焼室５の吸気弁閉時期ＩＶＣにおける圧力（圧縮開始時期圧力）Ｐｉｖｃ［ｋＰａ］を算出する。すなわち、吸気弁閉時期ＩＶＣになったタイミングでのコレクタ内圧力ＰＣＯＬを吸気弁閉時期ＩＶＣにおける圧力Ｐｉｖｃとして取り込む。
【０１２３】
ステップ５では、燃料噴射パルス幅Ｔｉ（または燃料噴射量）から図１３を内容とするテーブルを検索することにより、トータルの燃焼による圧力上昇分Δｐｔｏｔａｌ［Ｐａ］を算出する。
【０１２４】
ステップ６では、バルブリフト量Ｌｉｆｔから図４上段を内容とするテーブルを検索することにより、タンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］を算出する。すなわち、ＶＥＬ機構２８の非作動時には第１バルブリフト量補正係数ｋＶＥＬ１を、これに対してＶＥＬ機構２８の作動時には第２バルブリフト量補正係数ｋＶＥＬ２を、バルブリフト量補正係数ｋＶＥＬとして算出する。
【０１２５】
ステップ７では、同じくバルブリフト量Ｌｉｆｔから図１０上段を内容とするテーブルを検索することにより、流速係数ｋ１ＶＥＬ［無名数］を算出する。すなわち、ＶＥＬ機構２８の非作動時には第１流速係数ｋ１ＶＥＬ１を、これに対してＶＥＬ機構２８の作動時には第２流速係数ｋ１ＶＥＬ２を、流速係数ｋ１ＶＥＬとして算出する。
【０１２６】
ステップ８では、このうちの流速係数ｋ１ＶＥＬ、吸気弁閉時期ＩＶＣ、総ガス質量ＭＧＡＳ、エンジン回転速度Ｎｅを用いて、次の（４７）式により燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％［ｄｅｇＡＴＤＣ］を算出する。
【０１２７】
θｍｂ１０％＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ１・Ｎｅ＾ｃ２
＋ｃ３                    …（４７）
ただし、ｋ１ＶＥＬ：流速係数［無名数］、
ＩＶＣ    ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θｔ０    ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ＭＧＡＳ  ：総ガス質量［ｇ］、
Ｎｅ      ：エンジン回転速度［ｒｐｍ］、
ｃ１〜ｃ３：適合係数、
（４７）式は、上記（１４）式において、燃焼室内ガス質量ｍとして総ガス質量ＭＧＡＳを用いたもので、基本的に上記（１４）式と変わらない。
【０１２８】
ここで、燃焼圧が最大となるクランク角θｐｍａｘは（４７）式右辺第１項、つまり次の（４７−１）式により与えていることになる。
【０１２９】
θｐｍａｘ＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ１・Ｎｅ＾ｃ２
…（４７−１）
ステップ９では燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％から単純にクランク角で４ｄｅｇを差し引いた値を燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％［ｄｅｇＡＴＤＣ］として求める。これはθｍｂ２％よりθｍｂ１０％までのクランク角区間で燃焼質量割合が直線的に変化するとみなしてθｍｂ６％を求めるようにしたものある。
【０１３０】
ステップ１０では、このθｍｂ６％から図１２を内容とするテーブルを検索することにより、燃焼室５のθｍｂ６％時における容積Ｖｍｂ６％を算出する。
【０１３１】
ステップ１１、１２はステップ９、１０と同様である。 ステップ１１では燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％から単純にクランク角で８ｄｅｇを差し引いた値を燃焼質量割合が２％となるときのクランク角位置であるθｍｂ２％［ｄｅｇＡＴＤＣ］として求める。ステップ１２ではこのθｍｂ２％から図１２を内容とするテーブルを検索することにより、燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％を算出する。
【０１３２】
図１６は燃焼室５内の内部不活性ガス率ＭＲＥＳＲ［％］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。このフローは上記図１４のフローに先立って実行する。
【０１３３】
ステップ２１ではエアフローメータ３２の出力と目標当量比Ｔｆｂｙａを読み込む。ステップ２２ではエアフロメータ３２の出力に基づいて、燃焼室５に流入する新気量（シリンダ新気量）ＭＡＣＹＬ［ｇ］を算出する。このシリンダ新気量ＭＡＣＹＬの算出方法については公知の方法を用いればよい（特開２００１−５００９１号公報参照）。
【０１３４】
ステップ２３では、燃焼室５内の内部不活性ガス量ＭＲＥＳ［ｇ］を算出する。この内部不活性ガス量ＭＲＥＳの算出についても公知の方法を用いればよい（特開２００５−１７１８５６号公報参照）。
【０１３５】
ステップ２４では、この内部不活性ガス量ＭＲＥＳ、シリンダ新気量ＭＡＣＹＬ、目標当量比Ｔｆｂｙａから次の（４８）式により燃焼室５の総ガス質量ＭＧＡＳ［ｇ］を算出する。
【０１３６】
ＭＧＡＳ＝ＭＡＣＹＬ×(１＋Ｔｆｂｙａ／１４．７)＋ＭＲＥＳ
…（４８）
（４８）式右辺の括弧内の「１」は新気分、「Ｔｆｂｙａ／１４．７」は燃料分である。
【０１３７】
ステップ２５では、この総ガス質量ＭＧＡＳ、内部不活性ガス量ＭＲＥＳを用いて、次の（４９）式により内部不活性ガス率ＭＲＥＳＲ（燃焼室５内の総ガス質量に対する内部不活性ガス量の割合）［％］を算出する。
【０１３８】
ＭＲＥＳＲ＝ＭＲＥＳ／ＭＧＡＳ                          …（４９）
図１７は燃焼時間ＢＴ［ｓｅｃ］を算出するためのもの、図１８は着火遅れ時間τ［ｓｅｃ］を算出するためのもので、それぞれ一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１７、図１８は図１４に続けて実行する。図１７、図１８はどちらを先に実行してもかまわない。
【０１３９】
ここでは図１７から先に説明すると、ステップ３１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、図１４のステップ２で算出されている燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ［ｍ3］、図１４のステップ３で算出されている燃焼室５の吸気弁閉時期における温度Ｔｉｖｃ［Ｋ］、図１４のステップ４で算出されている燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ［Ｐａ］、図１４のステップ５で算出されているトータルの燃焼による圧力上昇分Δｐｃｔｏｔａｌ［ｋＰａ］、図１４のステップ６で算出されているタンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］、図１４のステップ１０で算出されている燃焼室５のθｍｂ６％時における容積Ｖｍｂ６％［ｄｅｇＡＴＤＣ］、図１６のステップ２４で算出されている総ガス質量ＭＧＡＳ［ｇ］、図１６のステップ２５で算出されている内部不活性ガス率ＭＲＥＳＲ［％］、燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］、目標当量比Ｔｆｂｙａ［無名数］、エンジン回転速度Ｎｅ［ｒｐｍ］を読み込む。
【０１４０】
ステップ３２では、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ、燃焼室５のθｍｂ６％時における容積Ｖｍｂ６％から、次の（５０）式により、容積変化による圧力変化分Δｐｖ［ｋＰ］を算出する。
【０１４１】
Δｐｖ＝Ｐｉｖｃ・｛（Ｖｉｖｃ／Ｖｍｂ６％）＾κ−１｝  …（５０）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（５０）式は上記（２８）式と同じものである。ここでの容積変化による圧力上昇分Δｐｖはθｍｂ２％よりθｍｂ１０％までのクランク角区間（燃焼期間）における圧力上昇分であるので、θｍｂ２％とθｍｂ１０％の平均値であるθｍｂ６％を用いている。
【０１４２】
ステップ３３では、燃焼によるトータルの圧力上昇分Δｐｃｔｏｔａｌから、次の（５１）式により、燃焼による圧力上昇分Δｐｃ［ｋＰ］を算出する。
【０１４３】
Δｐｃ＝０．０６×Δｐｔｏｔａｌ                        …（５１）
（５１）式は上記（２９）式において、燃焼質量割合ｘとして６％を用いたものである。
【０１４４】
ステップ３４では、これら燃焼による圧力上昇分Δｐｃ、容積変化による圧力上昇分Δｐｖを燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃに加算して、つまり次の（５２）式により燃焼室５のθｍｂ６％時における未燃ガス圧力Ｐｍｂ６％［ｋＰ］を算出する。
【０１４５】
Ｐｍｂ６％＝Ｐｉｖｃ＋Δｐｖ＋Δｐｃ                    …（５２）
ステップ３５では、この燃焼室５のθｍｂ６％時における未燃ガス圧力Ｐｍｂ６％、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における温度Ｔｉｖｃを用いて、次の（５３）式により燃焼室５のθｍｂ６％時における未燃ガス温度Ｔｍｂ６％［Ｋ］を算出する。
【０１４６】
Ｔｍｂ６％＝Ｔｉｖｃ・（Ｐｍｂ６％／Ｐｉｖｃ）＾｛（κ−１）／κ｝                                                                      …（５３）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（５３）式は上記（２５）式において燃焼質量割合ｘが６％のときの値である。
【０１４７】
ステップ３６では、目標当量比Ｔｆｂｙａから次の（５４）式〜（５６）式により標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂ［無名数］を求め、ステップ３７で、これら標準状態での層流燃焼速度ＳＬ０、係数ａ、ｂと、燃焼室５のθｍｂ６％時における未燃ガス温度Ｔｍｂ６％、燃焼室５のθｍｂ６％時における未燃ガス圧力Ｐｍｂ６％、内部不活性ガス率ＭＲＥＳＲを用いて、次の（５７）式によりθｍｂ６％時の層流燃焼速度ＳＬ［ｍ／ｓ］を算出する。
【０１４８】
ＳＬ０＝（０．２６３２−０．８４７２／（Ｔｆｂｙａ−１．１３）＾２）
…（５４）
ａ＝２．１８−０．８０・（Ｔｆｂｙａ−１）              …（５５）
ｂ＝−０．１６＋０．２２・（Ｔｆｂｙａ−１）            …（５６）
ＳＬ＝ＳＬ０・（Ｔｍｂ６％／２９８）＾ａ・（Ｐｍｂ６％／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ）                …（５７）
（５４）式〜（５６）式は、上記（２０）式〜（２２）式において、当量比φとして目標当量比Ｔｆｂｙａを用いたもので、基本的に（２０）式〜（２２）式と変わらない。同様に、（５７）式は、上記（１７）式において、未燃ガス温度Ｔ、燃焼室内圧力ｐとしてＴｍｂ６％、Ｐｍｂ６％を用いたもので、基本的に（１７）式と変わらない。
【０１４９】
ステップ３８では、このθｍｂ６％時の層流燃焼速度ＳＬ、タンブル強度のバルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣ、総ガス質量ＭＧＡＳ、エンジン回転速度Ｎｅを用いて、次の（５８）式により、燃焼室内のガス流動としてのタンブル強度Ｉｔ［無名数］を算出する。
【０１５０】
Ｉｔ＝ｋＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ４・Ｎｅ＾（ｃ５−１）
…（５８）
ただし、θｔ０    ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｃ４、ｃ５：適合係数、
ＶＴＣ機構２７の作動で吸気弁閉時期ＩＶＣが遅角側に遅らされるほど、つまりタンブル開始角を起点としてＩＶＣまでのクランク角区間（ＩＶＣ−θｔ０）が大きくなるほど実際のタンブル強度が大きくなるが、（５８）式よれば、吸気弁閉時期ＩＶＣが遅れるとき（ＩＶＣの値が大きくなるとき）タンブル強度Ｉｔも大きくなるのであり、算出値としてのＩｔは実際のタンブル強度と良く一致することとなる。
【０１５１】
また、ＶＥＬ機構２８の作動で吸気弁のバルブリフト量が大きくなるほど実際のタンブル強度が大きくなるが、図４上段によればバルブリフト量Ｌｉｆｔが大きくなるときバルブリフト量補正係数ｋＶＥＬが大きくなり、（５８）式よれば、バルブリフト量補正係数ｋＶＥＬが大きくなるときタンブル強度Ｉｔも大きくなるのであり、算出値としてのＩｔは実際のタンブル強度と良く一致することとなる。
【０１５２】
さらに、（５８）式よれば、総ガス質量ＭＧＡＳが大きくなるほどタンブル強度Ｉｔが大きくなる。
【０１５３】
ステップ３９では、このタンブル強度Ｉｔと回転速度Ｎｅから次の（５９）式により燃焼室内ガスの乱れ強さｕ［ｍ／ｓｅｃ］を算出する。
【０１５４】
ｕ＝Ｉｔ・Ｎｅ                                          …（５９）
これは、タンブルが圧縮工程中に小渦に変換されることにより、燃焼室内ガスの乱れが強くなるので、タンブル強度Ｉｔが強いほど燃焼室内ガスの乱れが強くなると考え、タンブル強度Ｉｔ（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さｕを推定するものである。
【０１５５】
ステップ４０では、この燃焼室内ガスの乱れ強さｕから次の（６０）式により乱流燃焼速度ＳＴ［ｍ／ｓｅｃ］を算出する。
【０１５６】
ＳＴ＝［ｕ／ＳＬ］＾ｃ６                                …（６０）
ただし、ｃ６：適合係数、
ステップ４１では、このようにして求めた乱流燃焼速度ＳＴと、ステップ３７で算出している、燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％時の層流燃焼速度ＳＬとから、次の（６１）式により燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％時の燃焼速度Ｓｂ［ｍ／ｓ］を算出する。
【０１５７】
Ｓｂ＝ＳＬ＋ＳＴ                                        …（６１）
なお、ステップ３８〜４１に代えて、燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％時の燃焼速度Ｓｂ［ｍ／ｓ］を次の一つの式により算出してもかまわない。
【０１５８】
Ｓｂ＝ＳＬ＋［｛ｋＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ４・Ｎｅ＾ｃ５｝
／ＳＬ］＾ｃ６                                    …（６２）
ただし、θｔ０          ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｃ４、ｃ５、ｃ６：適合係数、
（６２）式は、上記（１６）式において、燃焼室内ガス質量ｍとして総ガス質量ＭＧＡＳを用いたもので、基本的に（１６）式と変わらない。
【０１５９】
ステップ４２では、この燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％時の燃焼速度Ｓｂ、燃焼室５のθｍｂ６％時におけ容積Ｖｍｂ６％を用いて、次の（６３）式により燃焼時間ＢＴ［ｓ］を算出する。
【０１６０】
ＢＴ＝０．０７０７・Ｖｍｂ６％／（Ａｂ・Ｓｂ）          …（６３）
ただし、Ａｂ：火炎表面積［ｍ2］、
（６３）式は、上記（１５）式において、燃焼室内容積ＶｃｙｌとしてＶｍｂ６％を用いたもので、基本的に（１５）式と変わらない。
【０１６１】
図１８は着火遅れ時間τ［ｓｅｃ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。
【０１６２】
ステップ５１〜５７は図１６のステップ３１〜３７と同様である。すなわち、ステップ５１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、図１４のステップ２で算出されている燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ［ｍ3］、図１４のステップ３で算出されている燃焼室５の吸気弁閉時期における温度Ｔｉｖｃ［Ｋ］、図１４のステップ４で算出されている燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ［Ｐａ］、図１４のステップ５で算出されているトータルの燃焼による圧力上昇分Δｐｃｔｏｔａｌ［Ｐａ］、図１４のステップ６で算出されているタンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］、図１４のステップ１２で算出されている燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％［ｄｅｇＡＴＤＣ］、図１６のステップ２４で算出されている総ガス質量ＭＧＡＳ［ｇ］、図１６のステップ２５で算出されている内部不活性ガス率ＭＲＥＳＲ［％］、燃料噴射パルス幅Ｔｉ［ｍｓ］、目標当量比Ｔｆｂｙａ［無名数］、エンジン回転速度Ｎｅ［ｒｐｍ］を読み込む。
【０１６３】
ステップ５２では、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ、燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％から、次の（６４）式により、容積変化による圧力変化分Δｐｖ［ｋＰ］を算出する。
【０１６４】
Δｐｖ＝Ｐｉｖｃ・｛（Ｖｉｖｃ／Ｖｍｂ２％）＾κ−１｝  …（６４）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（６４）式は上記（２８）式と同じものである。ここでの容積変化による圧力上昇分Δｐｖはθｍｂ０％よりθｍｂ２％までのクランク角区間（着火遅れ期間）における圧力上昇分であるが、θｍｂ０％とθｍｂ２％の平均値であるθｍｂ１％は用いず、着火遅れ期間の終期の値であるθｍｂ２％を用いている。
【０１６５】
ステップ５３では、燃焼によるトータルの圧力上昇分Δｐｃｔｏｔａｌから、次の（６５）式により、燃焼による圧力上昇分Δｐｃ［ｋＰ］を算出する。
【０１６６】
Δｐｃ＝０．０２×Δｐｔｏｔａｌ                        …（６５）
（６５）式は上記の（２９）式において、燃焼質量割合ｘとして２％を用いたものである。
【０１６７】
ステップ５４では、これら燃焼による圧力上昇分Δｐｃ、容積変化による圧力上昇分Δｐｖを燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃに加算して、つまり次の（６６）式により燃焼室５のθｍｂ２％時における未燃ガス圧力Ｐｍｂ２％［ｋＰ］を算出する。
【０１６８】
Ｐｍｂ２％＝Ｐｉｖｃ＋Δｐｖ＋Δｐｃ                    …（６６）
ステップ５５では、この燃焼室５のθｍｂ２％時における未燃ガス圧力Ｐｍｂ２％、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における温度Ｔｉｖｃを用いて、次の（６７）式により燃焼室５のθｍｂ２％時における未燃ガス温度Ｔｍｂ２％［Ｋ］を算出する。
【０１６９】
Ｔｍｂ２％＝Ｔｉｖｃ・（Ｐｍｂ２％／Ｐｉｖｃ）＾｛（κ−１）／κ｝                                                                      …（６７）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（６７）式は上記（２５）式において燃焼質量割合ｘが２％のときの値である。
【０１７０】
ステップ５６では、目標当量比Ｔｆｂｙａから次の（６８）式〜（７０）式により標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂを求め、ステップ５７で、これら標準状態での層流燃焼速度ＳＬ０、係数ａ、ｂと、燃焼室５のθｍｂ２％時における未燃ガス温度Ｔｍｂ２％、燃焼室５のθｍｂ２％時における未燃ガス圧力Ｐｍｂ２％、内部不活性ガス率ＭＲＥＳＲを用いて、次の（７１）式によりθｍｂ２％時の層流燃焼速度ＳＬ［ｍ／ｓ］を算出する。
【０１７１】
ＳＬ０＝（０．２６３２−０．８４７２／（Ｔｆｂｙａ−１．１３）＾２）
…（６８）
ａ＝２．１８−０．８０・（Ｔｆｂｙａ−１）              …（６９）
ｂ＝−０．１６＋０．２２・（Ｔｆｂｙａ−１）            …（７０）
ＳＬ＝ＳＬ０・（Ｔｍｂ２％／２９８）＾ａ・（Ｐｍｂ２％／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ）                …（７１）
（６８）式〜（７０）式は、上記（２０）式〜（２２）式において、当量比φとして目標当量比Ｔｆｂｙａを用いたもので、基本的に（２０）式〜（２２）式と変わらない。同様に、（７１）式は、上記（１９）式において、未燃ガス温度Ｔ、燃焼室内圧力ｐとしてＴｍｂ２％、Ｐｍｂ２％を用いたもので、基本的に（１９）式と変わらない。
【０１７２】
ステップ５８では、このθｍｂ２％時の層流燃焼速度ＳＬ、燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％、エンジン回転速度Ｎｅ、総ガス質量ＭＧＡＳを用いて、次の（７２）式により着火遅れ時間τ［ｓｅｃ］を算出する。
【０１７３】
τ＝｛０．４９０・Ｖｍｂ２％＾（１／３）｝
／｛ｃ７・ＳＬ＋ｃ８・Ｎｅ・ＭＧＡＳ＾（１／２）｝
…（７２）
ただし、ｃ７、ｃ８：適合係数、
（７２）式は、上記（１８）式において、燃焼室内容積Ｖｃｙｌ、燃焼室内ガス質量ｍとして、それぞれ燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％、総ガス質量ＭＧＡＳを用いたもので、基本的に（１８）式と変わらない。
【０１７４】
図１９は基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１７、図１８のうち遅く実行されるフローに続けて実行する。
【０１７５】
ステップ６１では、図１４のステップ８で算出されている燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％［ｄｅｇＡＴＤＣ］、図１７のステップ４２で算出されている燃焼時間ＢＴ［ｓｅｃ］、図１８のステップ５８で算出されている着火遅れ時間τ［ｓｅｃ］、エンジン回転速度Ｎｅ［ｒｐｍ］を読み込む。
【０１７６】
ステップ６２では、着火遅れ時間τ、燃焼期間ＢＴ、燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％、エンジン回転速度Ｎｅを用いて、次の（７３）式により基本点火時期ＭＢＴＣＡＬ［ｄｅｇＡＴＤＣ］を算出し、ステップ６３で次の（７４）式によりその絶対値を採り、マイナスの符号を付けたものを改めて基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］として算出する。
【０１７７】
ＭＢＴＣＡＬ＝（τ＋ＢＴ）・Ｎｅ・６−θｍｂ１０％      …（７３）
ＭＢＴＣＡＬ＝−｜ＭＢＴＣＡＬ｜                        …（７４）
（７３）式、（７４）式は、上記（１３）式と同じ式である。
【０１７８】
このようにして算出した基本点火時期ＭＢＴＣＡＬは、点火時期指令値として点火レジスタに移され、実際のクランク角がこの点火時期指令値と一致したタイミングでエンジンコントローラ３１より一次電流を遮断する点火信号が点火コイル１３に出力される。
【０１７９】
ここで、本実施形態の作用効果を説明する。
【０１８０】
タンブル強度Ｉｔ（燃焼室内のガス流動）は、燃焼室内のさまざまな場所に分布する小渦の角運動量の総和で表されるので、例えばシリンダ（燃焼室５）の中心からｉ番目の小渦までの距離をｒi、ｉ番目の小渦のガス質量をｍi、ｉ番目の小渦のガス流速をｖiとすると、上記（３）式によりタンブル強度Ｉｔ（燃焼室内のガス流動）が表されるのであるが、燃焼室内の全ての小渦についてこれらの値ｍi、ｒi、ｖiをエンジンコントローラ３１においてオンラインで算出することは現在のところ不可能である。
【０１８１】
そこで、本発明者は、タンブルが圧縮工程中に小渦に変換されることにより、燃焼室内ガスの乱れｕが強くなるため、タンブル強度Ｉｔが強いほど燃焼室内ガスの乱れが強くなると考え、タンブル強度Ｉｔ（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さｕを推定するものとした。すなわち、本実施形態（請求項１、１０に記載の発明）によれば、タンブル強度Ｉｔ（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さｕを算出し（図１７のステップ３８、３９参照）、この燃焼室内ガスの乱れ強さｕに基づいて乱流燃焼速度ＳＴを算出し（図１７のステップ４０参照）、この乱流燃焼速度ＳＴに基づいてＭＢＴの得られる基本点火時期を算出する（図１７のステップ４１、４２、図１９のステップ６１〜６３参照）ので、タンブル強度Ｉｔ（燃焼室内のガス流動）の影響を受けて乱流燃焼速度ＳＴが変動しても、その変動する乱流燃焼速度ＳＴを精度よく算出できることから、基本点火時期ＭＢＴＣＡＬ（ＭＢＴの得られる基本点火時期）の予測精度を向上させることができる。
【０１８２】
上記のように、燃焼室内の全ての小渦についてｍi、ｒi、ｖiをエンジンコントローラ３１においてオンラインで算出することは現在のところ不可能であるので、タンブル強度Ｉｔを求めるには何らかの近似が必要である。そこで、本発明者は燃焼室内のガス流動を新たに吸気弁１５のリフト特性とエンジン回転速度Ｎｅとに基づいて推定するものとした。ここで、吸気弁１５のリフト特性とは、吸気弁のバルブリフト量（あるいは吸気弁のバルブ作動角）または吸気弁の開閉時期のことである。
【０１８３】
すなわち、本実施形態（請求項２、１１に記載の発明）によれば、タンブル強度Ｉｔ（燃焼室内のガス流動）を吸気弁１５のバルブリフト量Ｌｉｆｔ及び吸気弁閉時期ＩＶＣ（リフト特性）とエンジン回転速度Ｎｅとに基づいて推定するので（図１４のステップ６、図１７のステップ３８参照）、簡易な構成でタンブル強度Ｉｔ（燃焼室内のガス流動）を精度よく算出できることになった。
【０１８４】
ＶＥＬ機構２８を備えるときには、ＶＥＬ機構２８の作動状態で吸気弁１５のリフト量が小さいときと、ＶＥＬ機構２８の非作動状態で吸気弁１５のリフト量が大きいときでは燃焼室内のガス流動の向きが変わるため（図２（ｆ）と図２（ｇ）を参照）、燃焼室内ガスの乱れ強さｕも変わってしまうのであるが、本実施形態（請求項３、１２に記載の発明）によれば、乱流燃焼速度ＳＴを少なくとも吸気弁１５のバルブリフト量Ｌｉｆｔ（リフト特性）に基づいて算出するので（図１４のステップ６、図１７のステップ３８、３９、４０参照）、ＶＥＬ機構２８の作動、非作動に伴う吸気弁１５のバルブリフト量Ｌｉｆｔの大小に関係なく乱流燃焼速度ＳＴの推定精度を向上できる。
【０１８５】
ＶＴＣ機構２７の働きにより、吸気弁閉時期ＩＶＣが遅いほど燃焼室内ガスの乱れを形成する期間が長くなり、燃焼室内ガスの乱れ強さｕが強くなることに対応し、本実施形態（請求項４、１３に記載の発明）によれば、乱流燃焼速度ＳＴを少なくとも吸気弁閉時期ＩＶＣに基づいて算出するので（図１７のステップ３８、３９、４０参照）、吸気弁閉時期ＩＶＣの進角、遅角に関係なく乱流燃焼速度ＳＴの推定精度を向上できる。
【０１８６】
燃焼室５の総ガス質量ＭＧＡＳ（燃焼室内のガス量）が多いほどガスの持つ運動エネルギーが多くなり燃焼室内ガスの乱れ強さｕも強くなり燃焼が早くなることに対応し、本実施形態（請求項７、１６に記載の発明）によれば、乱流燃焼速度ＳＴを燃焼室５の総ガス質量ＭＧＡＳ（燃焼室内のガス量）に基づいても算出するので（図１７のステップ３８、３９、４０参照）、燃焼室５の総ガス質量ＭＧＡＳ（燃焼室内のガス量）の多少に関係なく乱流燃焼速度ＳＴの推定精度を向上できる。
【０１８７】
燃焼室内ガスの乱れが少ない状態では層流燃焼速度ＳＬが乱流燃焼速度ＳＴに与える影響が大きいことに対応し、本実施形態（請求項第８、１７の発明によれば、層流燃焼速度ＳＴを算出する層流燃焼速度算出処理手順を含み（図１７のステップ３７参照）、乱流燃焼速度ＳＴをこの層流燃焼速度ＳＬに基づいても算出するので（図１７のステップ４０参照）、燃焼室内ガスの乱れが少ない状態においても乱流燃焼速度ＳＴの推定精度を向上できる。
【０１８８】
ところで、実施形態にいう吸気弁閉時期ＩＶＣは、実際に吸気弁１５が閉じる瞬間のことである。実際には、吸気管内圧力の脈動による影響で、吸気弁閉時期ＩＶＣよりも早く燃焼室内圧力が上昇し始めると、吸気弁は開いているのにそれ以上空気が燃焼室５へと入らないことがある。こうした場合には、燃焼室内のガス流動（燃焼室内ガスの乱れ）に影響するのは、吸気弁閉時期ＩＶＣではなく、吸気弁閉時期ＩＶＣの前に燃焼室内圧力が実際に上昇し始める時期である。ここで、吸気弁閉時期ＩＶＣの前に燃焼室内圧力が実際に上昇し始める時期を「圧縮行程開始時期」として新たに定義すれば、吸気弁閉時期ＩＶＣと圧縮行程開始時期とは異なる値となる。
【０１８９】
そこで、吸気弁閉時期ＩＶＣに代え、この新たに定義した圧縮行程開始時期に基づいて乱流燃焼速度ＳＴを算出することで、吸気管内圧力の脈動による影響により、吸気弁閉時期ＩＶＣよりも早く燃焼室内圧力が上昇し始め、吸気弁１５は開いているのにそれ以上空気が燃焼室内に入らない場合においても、乱流燃焼速度ＳＴの推定精度を向上できることとなる（請求項５、６、１４、１５に記載の発明）。
【０１９０】
圧縮行程開始時期を求める方法には、次のように検出する方法と推定する方法の２つが考えられる。
【０１９１】
検出する方法：燃焼室に臨んで燃焼室内の圧力を検出するセンサ（圧縮行程開始時期検出装置）を設けておき、このセンサにより検出される燃料室内圧力に基づいて、吸気弁閉時期付近で、燃焼室内圧力が上昇を開始するときのクランク角を圧縮行程開始時期として計測（検出）する。
【０１９２】
推定方法：吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］と圧縮行程開始時期［ｄｅｇＡＴＤＣ］の間のクランク角度差ＤＣＡ［ｄｅｇ］はエンジン回転速度Ｎｅの感度を持っているので、次の式により圧縮行程開始時期を算出（推定）することができる。
【０１９３】
圧縮行程開始時期＝ＩＶＣ−ＤＣＡ                          …（７５）
（７５）式のクランク角度差ＤＣＡは、そのときのエンジン回転速度Ｎｅから図２０を内容とするテーブルを検索することにより求める。
【０１９４】
実施形態に対して、さらに、スワールコントロールバルブ、タンブルコントロールバルブのように、燃焼室５への吸気流動を変化させ得る可変吸気流動装置を備えさせることが考えられる。例えばスワールコントロールバルブを作動させた状態では、スワールコントロールバルブを作動させていない状態より、燃焼室内ガスの乱流燃焼速度ＳＴが大きくなるので、スワールコントロールバルブを作動させるときには、バルブの作動に伴って乱流燃焼速度ＳＴが大きくなる分だけ、乱流燃焼速度ＳＴが大きくなる側に補正してやればよい。すなわち、燃焼室内ガスの乱流燃焼速度ＳＴをスワールコントロールバルブの作動状態（作動状態、非作動状態）（可変吸気流動装置の状態）に応じても算出させる（請求項９、１８に記載の発明）。これにより、可変吸気流動装置を備えている場合においても、燃焼室内ガスの乱流燃焼速度ＳＴを精度よく算出できる。
【０１９５】
実施形態では、ＶＥＬ機構２８及びＶＴＣ機構２７を備える場合で説明したが、本発明は、ＶＥＬ機構２８、ＶＴＣ機構２７のいずれかを備える場合あるいはいずれも備えない場合にも適用がある。ＶＥＬ機構２８、ＶＴＣ機構２７のいずれも備えない場合には、バルブリフト量補正係数ｋＶＥＬ及び吸気弁閉時期ＶＴＣが固定値となるだけある。
【０１９６】
請求項１に記載の発明において、乱れ強さ算出処理手順は図１７のステップ３８、３９により、乱流燃焼速度算出処理手順は図１７のステップ４０により、基本点火時期算出処理手順は図１９のステップ６１、６２、６３によりそれぞれ果たされている。
【０１９７】
請求項１０に記載の発明において、乱れ強さ算出手段の機能は図１７のステップ３８、３９により、乱流燃焼速度算出手段の機能は図１７のステップ４０により、基本点火時期算出手段の機能は図１９のステップ６１、６２、６３によりそれぞれ果たされている。
【図面の簡単な説明】
【０１９８】
【図１】本発明の第１実施形態のエンジンの点火時期制御装置の概略構成図。
【図２】燃焼室内におけるタンブルの生成から消滅までをピストンの動きに合わせて表す図。
【図３】吸気弁閉時期に対するタンブル強度の特性図。
【図４】バルブリフト量に対するバルブリフト量補正係数の特性図。
【図５】燃焼圧が最大となるクランク角θｐｍａｘとθｍｂ０％、θｍｂ１０％、θｍｂ６０％との関係を表す特性図。
【図６】燃焼圧が最大となるクランク角θｐｍａｘとエンジン回転速度との関係を表す特性図。
【図７】燃焼質量割合と、燃焼質量割合が０％、１０％、２０％、３０％、４０％、５０％、６０％となるときのクランク角位置であるθｍｂ０％、θｍｂ１０％、θｍｂ２０％、θｍｂ３０％、θｍｂ４０％、θｍｂ５０％、θｍｂ６０％をそれぞれθｐｍａｘで直線近似した場合のＲの２乗値との関係を表す特性図。
【図８】コンベンショナルエンジンに適用している先行点火時期制御方法における基準クランク角と、本実施形態の変更後点火時期制御方法における基準クランク角の違いを説明するための燃焼質量割合の特性図。
【図９】吸入空気量、エンジン回転速度に対する、燃焼圧が最大となるクランク角θｐｍａｘの特性図。
【図１０】バルブリフト量に対する流速係数の特性図。
【図１１】第１実施形態の基本点火時期の算出方法を説明するための特性図。
【図１２】燃焼質量割合がｘ％となるときのクランク角位置であるθｍｂ ｘ％と、燃焼室５のθｍｂ ｘ％時における未燃ガス容積Ｖｍｂ ｘ％との関係を示す特性図。
【図１３】燃料噴射パルス幅と、燃焼によるトータルの圧力上昇分との関係を示す特性図。
【図１４】物理量の算出を説明するためのフローチャート。
【図１５】エンジンのクランクシャフトとコネクティングロッドの位置関係を説明するダイアグラム。
【図１６】内部不活性ガス率の算出を説明するためのフローチャート。
【図１７】燃焼期間の算出を説明するためのフローチャート。
【図１８】着火遅れ時間の算出を説明するためのフローチャート。
【図１９】基本点火時期の算出を説明するためのフローチャート。
【図２０】第２実施形態のクランク角度差の特性図。
【符号の説明】
【０１９９】
１  エンジン
５  燃焼室
１１  点火装置
１５  吸気弁
２１  燃料インジェクタ
２７  ＶＴＣ機構（可変動弁装置）
２８  ＶＥＬ機構（可変動弁装置）
３１  エンジンコントローラ
３３、３４  クランク角センサ
【図１】
0004747763.tif 000002
【図２】
0004747763.tif 000003
【図３】
0004747763.tif 000004
【図４】
0004747763.tif 000005
【図５】
0004747763.tif 000006
【図６】
0004747763.tif 000007
【図７】
0004747763.tif 000008
【図８】
0004747763.tif 000009
【図９】
0004747763.tif 000010
【図１０】
0004747763.tif 000011
【図１１】
0004747763.tif 000012
【図１２】
0004747763.tif 000013
【図１３】
0004747763.tif 000014
【図１４】
0004747763.tif 000015
【図１５】
0004747763.tif 000016
【図１６】
0004747763.tif 000017
【図１７】
0004747763.tif 000018
【図１８】
0004747763.tif 000019
【図１９】
0004747763.tif 000020
【図２０】
0004747763.tif 000021