JP4774893B2 - エンジンの点火時期制御方法及びエンジンの点火時期制御装置 - Google Patents

Description

本発明は、エンジン（内燃機関）の点火時期制御方法及びエンジンの点火時期制御装置に関する。

基準クランク角を、混合気の燃焼圧が最大となるクランク角θｐｍａｘとし、この基準クランク角と、燃焼室に供給された燃料に対する燃料質量の比率を表す燃焼質量割合が、ある割合に達するクランク角位置（具体的には燃料質量割合が６０％に達するクランク角位置であるθｍｂ６０％）とが等しいと仮定し、そのクランク角位置から燃焼期間ＢＵＲＮだけ進角側に計測したクランク角位置を、ＭＢＴの得られる点火時期として算出するものがある（特許文献１参照）。
特開２００４−３０１０４５公報

ところで、近年、吸気弁のリフト特性を変化させ得る可変動弁装置の開発が盛んに進められており、こうした可変動弁装置を備えるエンジンに対して、上記特許文献１の点火時期制御方法を適用したとき、燃焼圧が最大となるクランク角θｐｍａｘにおいて常に一定の燃焼質量割合になることはまれであることが実験により判明している。従って、可変動弁装置を備えるエンジンに対して、上記特許文献１の点火時期制御方法を単に適用しただけでは、点火時期算出値がＭＢＴの得られる点火時期から外れ燃費が悪化することが考えられる。

一方、机上で検討した結果によれば、燃焼圧が最大となるクランク角θｐｍａｘは、エンジンの損失との関係で変動し、エンジンの損失が大きくなるほど燃焼圧が最大となるクランク角θｐｍａｘが遅れるものと考えられる。この場合、エンジンの損失にはエンジンの冷却損失Ｑｃｏｏｌとエンジンの時間損失Ｑｔｉｍｅとが含まれる。

ここで、エンジンの冷却損失Ｑｃｏｏｌとは、燃焼室内で発生する燃焼熱のうち、エンジン冷却水により持ち去られる熱量のことである。エンジンの冷却損失Ｑｃｏｏｌを推定するための代表値としてはエンジンの負荷と回転速度Ｎｅとがあり、ウォッシーニの式を用いてエンジンの冷却損失Ｑｃｏｏｌを算出させるとすれば、エンジン回転速度Ｎｅが同じであればエンジンの負荷が大きくなるほどエンジンの冷却損失Ｑｃｏｏｌが大きくなり、またエンジン負荷が同じであればエンジン回転速度Ｎｅが大きくなるほどエンジンの冷却損失Ｑｃｏｏｌが小さくなる。このため、エンジン回転速度Ｎｅが同じであればエンジンの負荷が大きくなるほど燃焼圧が最大となるクランク角θｐｍａｘが遅れ、またエンジン負荷が同じであればエンジン回転速度Ｎｅが大きくなるほど燃焼圧が最大となるクランク角θｐｍａｘが進むこととなる。エンジンが高負荷であるほど燃焼圧が最大となるクランク角θｐｍａｘが遅れるのは、燃焼圧力に寄与すべき熱量のうちより冷却水へと奪われる熱量（つまり冷却損失）がエンジンの高負荷ほど大きくなるためである。

次に、エンジンの時間損失とは、燃焼が一瞬で終わらずクランク角で４０〜６０ｄｅｇ程度の燃焼期間を要するために発生する損失のことである。言い換えると、燃焼期間が長いほど実圧縮比が低下し、熱効率が低下するが、この熱効率の低下をエンジンの時間損失で表現している。エンジンの時間損失を推定するための代表値としては燃焼速度と燃焼期間とがあり、燃焼速度が小さくなれば、熱効率が低下するためエンジンの時間損失は大きくなり、燃焼期間が長くなれば（エンジン回転速度が小さくなる）、同じく熱効率が低下するためエンジンの時間損失は大きくなる。このため、エンジン回転速度が同じであれば燃焼速度が小さくなるほど燃焼圧が最大となるクランク角θｐｍａｘが遅れ、また燃焼速度が同じであれば燃焼期間が長い（エンジン回転速度が小さい）ほど燃焼圧が最大となるクランク角θｐｍａｘが遅れることとなる。エンジン回転速度が同じであれば燃焼速度が小さくなるほど燃焼圧が最大となるクランク角θｐｍａｘが遅れるのは、燃焼速度が小さいほど燃焼圧力に寄与しなくなるためである。

しかしながら、上記特許文献１の技術では、燃焼圧が最大となるクランク角θｐｍａｘを一定値で設定するかまたは回転速度に応じて算出している。従って、いま、燃焼圧が最大となるクランク角θｐｍａｘを設定したときのエンジン負荷及び回転速度の運転条件より回転速度は変化せずにエンジン負荷だけが大きくなった場合を考えてみると、燃焼圧が最大となるクランク角θｐｍａｘは、負荷増大前より遅角側にずれるにも拘わらず、ずれないとして点火時期を算出するため、このすれ分だけ点火時期算出値がＭＢＴの得られる点火時期より進角側の値となってしまう。この逆に、燃焼圧が最大となるクランク角θｐｍａｘを設定したときのエンジン負荷及び回転速度の運転条件より回転速度は変化せずに負荷だけが小さくなった場合を考えてみると、燃焼圧が最大となるクランク角θｐｍａｘは、負荷減少前より進角側にずれるにも拘わらず、ずれないとして点火時期を算出するため、このすれ分だけ点火時期算出値がＭＢＴの得られる点火時期より遅角側の値となる。

このように、燃焼圧が最大となるクランク角θｐｍａｘを一定値で設定したり回転速度のみの感度を持たせたのでは、ＭＢＴの得られる点火時期を精度良く算出できない事態が生じ得る。

そこで本発明は、特に吸気弁のリフト特性を変化させ得る可変動弁装置を備えるエンジンにおいても、ＭＢＴの得られる点火時期を精度良く算出することが可能な点火時期制御方法及び点火時期制御装置を提供することを目的とする。

本発明は、燃焼圧が最大となるクランク角θｐｍａｘを、エンジンの冷却損失、またはエンジンの時間損失に基づいて算出し、この燃焼圧が最大となるクランク角θｐｍａｘに基づいてＭＢＴの得られる基本点火時期を算出するように構成する。

燃焼圧が最大となるクランク角θｐｍａｘは運転条件によって変動するが、その値はエンジンの冷却損失、またはエンジンの時間損失によって決まることに鑑み、本発明によれば、燃焼圧が最大となるクランク角θｐｍａｘを、エンジンの冷却損失、またはエンジンの時間損失に基づいて算出し、この燃焼圧が最大となるクランク角θｐｍａｘに基づいてＭＢＴの得られる基本点火時期を算出するようにしたので、燃焼圧が最大となるクランク角θｐｍａｘが運転条件により変動しても、その変動するθｐｍａｘを精度よく推定できることから、ＭＢＴの得られる基本点火時期の予測精度を向上させることができる。

特に、吸気弁のリフト特性を変化させ得る可変動弁装置を備えるエンジンでは、可変動弁装置を備えていないエンジンよりもθｐｍａｘが大きく変動するので、こうしたエンジンに特に有効である。ただし、可変動弁装置を備えていないエンジンに対しても適用があることはもちろんである。

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

図１は、エンジンの点火時期制御方法の実施に直接使用するエンジンの点火時期制御装置の概略構成を示している。

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

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

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

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

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

吸気弁１５、排気弁１６は、クランクシャフト７を動力源として、各々吸気側カムシャフト２５及び排気側カムシャフト２６に設けられたカムの動作により開閉駆動される。吸気側には、吸気弁１５のバルブリフト量及び作動角を連続的に可変制御する多節リンク状の機構で構成される可変バルブ機構（以下、「ＶＥＬ機構」という。）２８を備える。このＶＥＬ機構２８には吸気弁１５のバルブリフト量及び作動角を検出する作動角センサ（図示しない）が併設されている。

同じく吸気側には、クランクシャフト７と吸気側カムシャフト２５との回転位相差を連続的に可変制御して、吸気弁１５のバルブタイミングを進遅角する可変バルブタイミング機構（以下、「ＶＴＣ機構」という。）２７を備える。また、吸気側カムシャフト２５の他端には吸気側カムシャフト２５の回転位置を検出するためのカム角センサ３４が併設されている。

これらＶＥＬ機構２８及びＶＴＣ機構２７（可変動弁装置）の具体的な構成は特開２００３−３８７２号公報により公知であるので、その詳しい説明は省略する。

ＶＥＬ機構２８、ＶＴＣ機構２７の各アクチュエータに指令して、吸気弁１５のリフト特性（バルブタイミング（開閉時期）や吸気弁１５のバルブリフト量）を変えると燃焼室５に残留する不活性ガスの量が変化する。燃焼室５内の不活性ガスの量が増えるほどポンピングロスが減って燃費がよくなるので、運転条件によりどのくらいの不活性ガスが燃焼室５内に残留したらよいかを目標吸気弁閉時期や目標バルブリフト量にして予め定めており、エンジンコントローラ３１ではそのときの運転条件（エンジンの負荷と回転速度）より目標吸気弁閉時期と目標バルブリフト量とを定め、それら目標値が得られるようにＶＴＣ機構２７及びＶＥＬ機構２８の各アクチュエータを介して吸気弁１５の閉時期とバルブリフト量とを制御する。

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

さて、ＭＢＴ（最大トルクの得られる最小進角値）で混合気に点火した場合に、混合気の燃焼圧力が最大値Ｐｍａｘとなるクランク角θｐｍａｘを基準クランク角とすると、基準クランク角は燃焼方式によらずほぼ一定である。また、燃焼室内における燃焼解析によれば、燃焼室５に供給された燃料に対する燃料質量の比率を表す燃焼質量割合は点火時に０％であり、完全燃焼によって１００％に達する。そして、基準クランク角における燃焼質量割合は一定で約６０％であるとして、そのときの運転条件（エンジンの負荷と回転速度）より燃焼速度を求め、この燃焼速度に基づいて燃焼期間ＢＵＲＮを算出し、基準クランク角よりこの燃焼期間ＢＵＲＮと着火遅れ時間相当角ＩＧＮＤＥＡＤとの合計のクランク角区間だけ進角側のクランク角位置を基本点火時期ＭＢＴＣＡＬとして算出する点火時期制御方法を提案している（特開２００４−３３２６４７号公報参照）。

燃焼解析に基づくこうした点火時期制御方法（この点火時期制御方法を、以下「先行点火時期制御方法」という。）は吸気弁１５のバルブリフト量及び吸気弁１５の開閉タイミングが変化しないエンジン（コンベンショナルエンジン）を対象に本出願人が開発してきた経過があるので、現在でもコンベンショナルエンジンに適用する限り先行点火時期制御方法でなんら問題ないのであるが、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンに対しても、この先行点火時期制御方法を適用したとき、上記の基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期と合わないことが判明している。その理由は、ＶＥＬ機構２８及びＶＴＣ機構２７が非作動状態であるときに点火時期制御に用いる各種の値を適合していれば、当然ながら、ＶＥＬ機構２８及びＶＴＣ機構２７が非作動状態であるときに基本点火時期ＭＢＴＣＡＬはＭＢＴの得られる点火時期と一致する。

この状態つまり同じ運転条件においてＶＥＬ機構２８やＶＴＣ機構２７を作動させてリフト特性（吸気弁１５のバルブリフト量や吸気弁閉時期）を変化させると、燃焼室５内のガス流動であるタンブルやスワールの各強度が変化し、これに伴って燃焼室内ガスの乱流状態での燃焼速度である乱流燃焼速度が変化する。すると、ＭＢＴの得られる点火時期も変化する。しかしながら、ＶＥＬ機構２８やＶＴＣ機構２７を作動させたからといって運転条件が同じであるため、基本点火時期ＭＢＴＣＡＬは変化しない。従って、ＶＥＬ機構２８やＶＴＣ機構２７を作動させた途端に基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期からずれて燃費が悪くなる。

例えば、ＶＥＬ機構２８の作動で吸気弁１５のバルブリフト量が大きい状態から小さい状態に切換えられたとすると、切換前よりタンブルやスワールが形成されにくくなり乱流燃焼速度が遅くなって燃焼期間が長引く。従って、同じ運転条件で考えると、バルブリフト量が小さい状態であるときには、バルブリフト量が大きい状態のときよりＭＢＴの得られる点火時期が進角側に移動するのであるが、バルブリフト量が小さい状態のときにも、バルブリフト量が大きい状態のときと同じ燃焼期間ＢＵＲＮを算出したのでは、実際より短い燃焼期間を算出してしまうことになり、基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期よりも遅すぎることとなる。

ＶＴＣ機構２７の作動で吸気弁閉時期が遅い状態から早い状態へと進角される場合も同様である。すなわち、吸気弁閉時期が遅い状態から早い状態へと進角されると、進角される前よりタンブルやスワールが形成されにくいために乱流燃焼速度が遅くなって燃焼期間が長くなる。従って、同じ運転条件で考えると、吸気弁閉時期が早い状態であるときには、吸気弁閉時期が遅い状態のときよりＭＢＴの得られる点火時期が進角側に移動するのであるが、吸気弁閉時期が早いときにも吸気弁閉時期が遅い状態のときと同じ燃焼期間ＢＵＲＮを算出したのでは、実際より短い燃焼期間を算出してしまうことになり、基本点火時期ＭＢＴＣＡＬがＭＢＴの得られる点火時期より遅すぎてしまう。

そこで本発明の参考例では、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象としていても、物理モデルに従った制御構造をもち、可能な限り適合の不要な制御とするため、コンベンショナルエンジンに適用している先行点火時期制御方法に対して〈１〉燃焼速度の算出方法、〈２〉基準クランク角の設定方法、〈３〉着火遅れ時間の算出方法を変更する。

ここで、上記〈１〉だけでなく、〈２〉と〈３〉をも追加しているのは、制御方法を再検討した結果である。まず、上記〈２〉をも追加しているのは次の理由による。すなわち、コンベンショナルエンジンに適用している先行点火時期制御方法に対してこれまで正しいとされてきた、基準クランク角と燃焼質量割合が６０％のときのクランク角位置との関係を、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンについて実験してみると、基準クランク角と燃焼質量割合が６０％のときのクランク角位置との間に大きなずれ（クランク角差）があり、そのずれがそのまま点火時期算出値（ＭＢＴＣＡＬ）の推定誤差になってしまうため、基準クランク角位置の設定方法を改める必要が生じたためである。

上記〈３〉をも追加しているのは次の理由による。〈３〉の着火遅れ時間とは、点火タイミングを起点として燃焼質量割合が０％である間の時間のことであるが、この着火遅れ時間の推定誤差が大きい。その理由は、燃焼解析装置によって燃焼質量割合が０％であるクランク角位置を計測したときの計測誤差が大きいために、この計測した値を正にして着火遅れ時間ＤＥＡＤＴＩＭＥの適合を行うと、燃焼期間ＢＵＲＮの推定誤差にこの計測誤差も含まれ誤差が大きくなってしまうためである。従って、燃焼解析装置による計測誤差が入ってこないように、着火遅れ時間の算出方法を改めるようにしたものである。

以下、項を分けて説明する。
〈１〉燃焼速度の算出方法
コンベンショナルエンジンに適用している先行点火時期制御方法では、燃焼速度Ｓｂを燃焼室内ガスの層流状態での燃焼速度である層流燃焼速度ＳＬと、燃焼室内ガスの乱流状態での燃焼速度である乱流燃焼速度ＳＴの和、つまり次の（補１）式により燃焼速度Ｓｂを求めている。

Ｓｂ＝ＳＬ＋ＳＴ …（補１）
（補１）式の層流燃焼速度ＳＬ、乱流燃焼速度ＳＴはそれぞれ次の（補２）式、（補３）式により与えている。

ＳＬ＝ＳＬ０×（Ｔ／２９８）＾ｎ×（ｐ／１０１．３２５）＾ｄ
×（１−ｂ×ＭＲＥＳＲ＾ｋ） …（補２）
ＳＴ＝Ｓｔ・Ｎｅ …（補３）
ただし、ＳＬ０ ：標準状態での層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ ：未燃ガス温度［Ｋ］、
ｐ ：燃焼室内圧力［ｋＰａ］、
ＭＲＥＳＲ ：内部不活性ガス率［％］、
Ｓｔ ：乱流燃焼速度係数、
ｎ、ｄ、ｂ、ｋ：係数、
（補３）式の乱流燃焼速度ＳＴは、燃焼室５内におけるガスの乱れ強さによって決まる。また、（補３）式のように、乱流燃焼速度Ｓｔをエンジン回転速度Ｎｅの関数としたのは、燃焼室５内におけるガスの乱れ強さは吸気弁１５のバルブリフト量、開閉タイミングが同じ場合、ピストンスピード（エンジン回転速度）に比例するため、乱流燃焼速度Ｓｔの推定にはエンジン回転速度Ｎｅの感度を持たせれば十分であると考えたものである。

しかしながら、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンでは、コンベンショナルエンジンと相違して、吸気弁１５のバルブリフト量や開閉タイミングが運転条件によって大きく変わり、それらバルブリフト量や開閉タイミングの変化による燃焼室５内におけるガスの乱れ強さの変化を無視できないので、本発明の参考例では、上記（補３）式に代え、次のようにして乱流燃焼速度を算出する。すなわち、ＶＥＬ機構２８の作動に伴う吸気弁１５のバルブリフト量の変化と、ＶＴＣ機構２７の作動に伴う吸気弁閉時期の変化とに対応するため、特にタンブルの影響を燃焼速度の推定式に加える。

まず、燃焼速度Ｓｂ［ｍ／ｓｅｃ］は次の（１）式により表される。（１）式右辺第２項は乱流燃焼速度ＳＴを与える式で、このように乱流燃焼速度ＳＴを、ｕ／ＳＬを変数としてこの変数の指数関数で与えるとする点は公知である。

Ｓｂ＝ＳＬ＋ｂ・（ｕ／ＳＬ）＾ａ …（１）
ただし、ＳＬ ：層流燃焼速度、
ｕ ：乱れ強さ、
ａ，ｂ：適合係数、
（１）式の、燃焼室内ガス流動に伴う乱れのうちの１つの特性値である乱れ強さｕ［ｍ／ｓｅｃ］は次の（２）式により表わされる。これは、燃焼室内ガスの乱れ強さｕは回転速度Ｎｅと比例の関係を持つことが、また、その比例定数は、スワール強度、タンブル強度と相関があることが知られているので、これを式に表したものである。

ｕ＝ｃ・Ｉｔ・Ｎｅ …（２）
ただし、Ｉｔ：タンブル強度［無名数］
Ｎｅ：エンジン回転速度［ｒｐｍ］、
ｃ ：適合係数、
ここで、本発明の参考例では（２）式にスワール強度は入れていない。これは、今回対象としているエンジンにおいては、ピストン６の中心軸に対して周方向に旋回する流れであるスワールは、ピストン６の中心軸に直交する軸に対して周方向に旋回する流れであるタンブルに比べ、乱流燃焼速度に与える影響が少ないと判断し、今回は考慮しないためである。

ただし、スワール強度は対象外というのではない。対象とするエンジンにおいてタンブルよりスワールのほうが乱流燃焼速度に与える影響が大きいときには、スワール強度Ｉｓを、後述するタンブル強度Ｉｔの推定方法と同様の推定方法により推定すればよい。

次に、（２）式のタンブル強度Ｉｔ（燃焼室内のガス流動）の推定方法を説明する。まず、図２（ａ）〜図２（ｅ）は燃焼室５内におけるタンブルの生成から消滅までをピストン６の動きに合わせて図解したものである。順に説明すると、図２（ａ）は吸気弁開時期ＩＶＯ、つまり吸気弁１５が開いた瞬間で、燃焼室５内においてさまざまな方向にガスが流れ込んでいる。図２（ｂ）はピストン６が下降している吸気上死点後４５ｄｅｇ付近（４５ｄｅｇＡＴＤＣ）を示し、この４５ｄｅｇＡＴＤＣ付近からタンブルが形成され始める。図２（ｃ）は吸気弁閉時期ＩＶＣ付近を示し、吸気弁閉時期ＩＶＣ付近までタンブル強度が増加していく。したがって、吸気弁閉時期ＩＶＣが遅いほどタンブルが強くなる。

図２（ｄ）は吸気上死点後３１５ｄｅｇ付近まで、つまりピストン６が上昇する圧縮行程中を示し、圧縮工程中に徐徐にタンブルが小渦（乱れ）に変換されていく。図２（ｅ）は吸気上死点後３１５ｄｅｇ付近を示し、この３１５ｄｅｇＡＴＤＣ付近でタンブルは消え全て小渦に変換される。

一方、図２（ｆ）と図２（ｇ）はＶＥＬ機構２８の作動、非作動に伴う吸気弁１５のバルブリフト量の相違でタンブルの形成がどのように違うのかを示している。このうち、左側に示す図２（ｆ）は吸気弁１５のバルブリフト量が小さい場合で、バルブリフト量が小さいときには、吸気がシリンダ壁側とシリンダ中心側の両方に分かれて流れ、お互いに流れを打ち消しあうため、タンブルが形成されにくい。これに対して右側に示す図２（ｇ）は吸気弁１５のバルブリフト量が大きい場合で、バルブリフト量が大きいときには、ほとんどの吸気がシリンダ中心側へと流れるため、タンブルを形成しやすい。

さて、タンブル強度Ｉｔ［無名数］は次式により小渦の有する角運動量の総和として算出されることがＳＡＥペーパーに記載されている（ＳＡＥ９８１０４８）。

Ｉｔ＝Σｍ_iｒ_i×ｖ_i …（３）
ただし、ｍ_i：ｉ番目の小渦におけるガスの質量、
ｒ_i：シリンダの中心からｉ番目の小渦までの距離、
ｖ_i：ｉ番目の小渦でのガス流速、
しかしながら、燃焼室内の全ての小渦についてこれらの値ｍ_i、ｒ_i、ｖ_iをエンジンコントローラ３１においてオンラインで算出することは現在のところ不可能であるので、上記（３）式に代えて次の（４）式でタンブル強度Ｉｔ［無名数］を近似する。

Ｉｔ＝ａ・（ｍ＾ｂ・Ｎｅ＾ｃ） …（４）
ただし、ｍ ：燃焼室内ガス質量［ｇ］、
Ｎｅ ：エンジン回転速度［ｒｐｍ］、
ａ ：係数、
ｂ、ｃ：適合係数、
（４）式のｍ＾ｂはタンブル強度に対する燃焼室内ガス質量ｍの、またＮｅ＾ｃはタンブル強度に対する回転速度Ｎｅの各影響を反映させたものである。ここで、燃焼室内ガス質量ｍは吸気弁閉時期ＩＶＣに燃焼室５に入っているガス質量のことである。

さらに、タンブルは上記図２（ａ）〜図２（ｅ）で示したように、吸気弁閉時期ＩＶＣが４５ｄｅｇＡＴＤＣ付近より遅いほど強くなり、また、図２（ｆ）、図２（ｇ）で示したように吸気弁１５のバルブリフト量が大きいほうが強くなると考え、（４）式の係数ａを新たに次の（５）式で表す。

ａ＝ｋＶＥＬ・（ＩＶＣ−θ０ｔ） …（５）
ただし、ｋＶＥＬ：タンブル強度のバルブリフト量補正係数［１／ｄｅｇ］、
ＩＶＣ ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θ０ｔ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
従って、（５）式を（４）式に代入した次の（６）式によりタンブル強度Ｉｔを近似する。

Ｉｔ＝ｋＶＥＬ・（ＩＶＣ−θ０ｔ）・ｍ＾ｂ・Ｎｅ＾ｃ …（６）
ただし、ｋＶＥＬ：タンブル強度のバルブリフト量補正係数［１／ｄｅｇ］、
ＩＶＣ ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θ０ｔ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｂ、ｃ ：適合係数、
（６）式の吸気弁閉時期ＩＶＣの単位は吸気上死点を起点として遅角側に計測する値［ｄｅｇＡＴＤＣ］であるので、吸気弁閉時期ＩＶＣが吸気上死点より遅れるほどＩＶＣの値が大きくなり、従って（６）式よりタンブル強度Ｉｔが大きくなる。また、（６）式のバルブリフト量補正係数ｋＶＥＬはバルブリフト量の関数であり、バルブリフト量が大きいときには、バルブリフト量が小さいときより大きくなり、従って（６）式よりタンブル強度Ｉｔが大きくなる。

（６）式によりタンブル強度Ｉｔを近似した新しい考え方が妥当かどうかを確認するため、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンについて、ＶＥＬ機構２８は非作動状態つまり吸気弁１５のバルブリフト量を一定として、ＶＴＣ機構２７のみを作動し吸気弁閉時期ＩＶＣを変化させる実験を行ったところ、図３に示すようにタンブル強度と吸気弁閉時期ＩＶＣとの関係を示す結果が得られた。図３より吸気弁１５のバルブリフト量一定の条件で、吸気弁閉時期ＩＶＣが大きくなる（つまりＩＶＣが吸気上死点より遅れる）ほどタンブル強度Ｉｔが強くなることが確かめられた。

上記（６）式のタンブル形成開始角θ０ｔは一定値（例えば４５ｄｅｇＡＴＤＣ）である。タンブル形成開始角θ０ｔはエンジン仕様に依存し、ＶＥＬ機構２８、ＶＴＣ機構２７の有無には依存しない。また、図２より、吸気弁閉時期ＩＶＣはタンブル形成開始角θ０ｔより必ず遅角側の値である。従って、上記（６）式の（ＩＶＣ−θ０ｔ）の値は必ず正の値になると考えている。

上記（６）式のタンブル強度のバルブリフト量補正係数ｋＶＥＬはバルブリフト量Ｌｉｆｔ［ｍ］（あるいはバルブ作動角）に応じて定めている。ここでは簡単のため、例えば、ＶＥＬ機構２８の作動、非作動により吸気弁１５のバルブリフト量が大小の２段に切換可能であり、通常の運転条件（運転条件１とする）ではＶＥＬ機構２８を非作動状態とし、運転条件１より所定の運転条件（運転条件２とする）に移行したときにはＶＥＬ機構２８を作動させ、図４下段に示したように実線のバルブリフト特性から破線のバルブリフト特性へとバルブリフト量を小さくしているとして説明すると、ＶＥＬ機構２８の非作動状態でバルブリフト量が第１バルブリフト量Ｌｉｆｔ１［ｍ］となり、ＶＥＬ機構２８の作動状態でバルブリフト量が第２バルブリフト量Ｌｉｆｔ２［ｍ］となる（Ｌｉｆｔ２＜Ｌｉｆｔ１）。従って、図４上段に示したように、ＶＥＬ機構２８の非作動時（あるいは運転条件１のとき）には第１バルブリフト量Ｌｉｆｔ１に対応する第１バルブリフト量補正係数ｋＶＥＬ１を、これに対してＶＥＬ機構２８の作動時（あるいは運転条件２のとき）には第２バルブリフト量Ｌｉｆｔ２に対応する第２バルブリフト量補正係数ｋＶＥＬ２をバルブリフト量補正係数ｋＶＥＬとして設定する。このように、ＶＥＬ機構２８の作動で吸気弁１５のバルブリフト量が第１バルブリフト量Ｌｉｆｔ１より第２バルブリフト量Ｌｉｆｔ２へと小さくなったときに、バルブリフト量補正係数ｋＶＥＬを第１バルブリフト量補正係数ｋＶＥＬ１より第２バルブリフト量補正係数ｋＶＥＬ２へと小さくするのは、バルブリフト量が小さいときのほうがタンブルが弱くなる（従って、タンブル強度Ｉｔが小さくなる）からである。

このように、（６）式によれば、タンブル強度Ｉｔを、新たに、吸気弁１５のバルブリフト量、吸気弁閉時期ＩＶＣ、燃焼室内ガス質量ｍ、エンジン回転速度Ｎｅの関数としている。
〈２〉基準クランク角の設定方法
コンベンショナルエンジンに適用している先行点火時期制御方法では、次の（補４）式のように、燃焼圧が最大となるクランク角θｐｍａｘを、燃焼質量割合が６０％となるときのクランク角位置に設定すると共に、この燃焼圧が最大となるクランク角θｐｍａｘを基準クランク角としている。なお、燃焼質量割合がｘ％となるときのクランク角位置を、以下、「θｍｂ ｘ％」で表す。例えば、燃焼質量割合が６０％となるときのクランク角位置は「θｍｂ６０％」である。また、θｍｂ ｘ％の起点は圧縮上死点とする。

θｍｂ６０％＝θｐｍａｘ …（補４）
また、燃焼圧が最大となるクランク角θｐｍａｘを次の（補５）式のようにエンジン回転速度Ｎｅの関数で与えている。

θｐｍａｘ＝ｈ・（Ｎｅ）＾ｉ …（補５）
ただし、ｈ、ｉ：適合係数、
一方、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象として実験してみたところ、図５に示したように、回転速度Ｎｅが一定の条件のもとで、燃焼圧が最大となるクランク角θｐｍａｘとθｍｂ ｘ％との関係を表す実験結果が得られた。ただし、図５にはθｍｂ ｘ％として燃焼質量割合が０％となるときのクランク角位置であるθｍｂ０％、燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％、燃焼質量割合が６０％となるときのクランク角位置であるθｍｂ６０％の３つの場合だけを示しており、図５において左はθｐｍａｘとθｍｂ０％の、中央はθｐｍａｘとθｍｂ１０％の、右はθｐｍａｘとθｍｂ６０％の各関係を整理したものである。実際には燃焼質量割合が０％、１０％、６０％以外の値（２０％、３０％、４０％、５０％）となるときのクランク角位置での実験結果も得ている（図７参照）。

燃焼圧が最大となるクランク角θｐｍａｘと、燃焼質量割合が６０％となるときのクランク角位置であるθｍｂ６０％とでは、図５右側に示したように最大誤差が２．６ｄｅｇもあり、θｐｍａｘとθｍｂ０％との間の誤差を示す図５左側と同様であることがわかる。

次に、図６は、吸入空気量（エンジン負荷）が一定の条件のもとでの、燃焼圧が最大となるクランク角θｐｍａｘとエンジン回転速度Ｎｅとの関係を示したものである。図示のように、傾向としては回転速度Ｎｅが大きくなるほどθｐｍａｘが進角側にずれているが、回転速度Ｎｅに対するθｐｍａｘのバラツキが大きく、燃焼圧が最大となるクランク角θｐｍａｘは上記（補５）式のようにエンジン回転速度Ｎｅだけでは表しきれないことがわかる。

図７は燃焼質量割合を横軸に、燃焼質量割合がｘ％となるときのクランク角位置であるθｍｂ ｘ％をθｐｍａｘで直線近似した場合のＲの２乗値を縦軸に採り、図５に示したθｍｂ０％時、θｍｂ１０％時、θｍｂ６０％時の３つの実験結果及び図示しない図５と同様のθｍｂ２０％時、θｍｂ３０％時、θｍｂ４０％時、θｍｂ５０％時の残りの実験結果とを整理したものである。縦軸のＲの２乗値は、この値が１．０に近づくほどθｐｍａｘとθｍｂ ｘ％（ｘ＝０、１０、２０、３０、４０、５０、６０）との間に相関があることを表すので、図７によれば、燃焼圧が最大となるクランク角θｐｍａｘはθｍｂ１０％付近で一番相関があることがわかる。

なお、図７にはＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンでもエンジン仕様の異なる２種類のエンジンの特性を示している。２つのエンジンの違いは主に燃焼室形状にあり、第２エンジンは第１エンジンよりも燃焼室が平べったい形状のものである。第１エンジンでは燃焼圧が最大となるクランク角θｐｍａｘは燃焼質量割合が１０％となるときのクランク角位置から５０％となるときのクランク角位置との間で広く相関があるのに対して、第２エンジンになると燃焼圧が最大となるクランク角θｐｍａｘは燃焼質量割合が１０％となるときのクランク角位置付近でだけ相関がある。従って、２つのエンジンに共通して燃焼圧が最大となるクランク角θｐｍａｘが高い相関を有するのは燃焼質量割合が１０％となるときのクランク角位置付近にあるときだけである。このように、エンジン機種が相違しても、燃焼圧が最大となるクランク角θｐｍａｘはθｍｂ１０％となるときのクランク角位置付近で一番相関がある、という結果が初めて得られた。

こうした実験結果を受けて本発明の参考例では、燃焼圧が最大となるクランク角θｐｍａｘに基づいて、燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％を算出し、このθｍｂ１０％を基準クランク角として設定する、つまり次の（７）式によりθｍｂ１０％［ｄｅｇＡＴＤＣ］を算出する。

θｍｂ１０％＝ａ・θｐｍａｘ−ｂ …（７）
ただし、ａ：適合係数［無名数］、
ｂ：適合クランク角［ｄｅｇ］、
（７）式は燃焼圧が最大となるクランク角θｐｍａｘをｂ（正の値）のクランク角だけ進角側にシフトしてθｍｂ１０％を求めるとする式である。この結果、図８に示したように本発明の参考例での基準クランク角（＝θｍｂ１０％）はコンベンショナルエンジンに適用している先行点火時期制御方法における基準クランク角（＝θｍｂ６０％）よりも進角側にくることとなる。

そして、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンでは、燃焼圧が最大となるクランク角θｐｍａｘは、吸気弁閉時期ＩＶＣが遅いほど、また吸気弁１５のバルブリフト量が大きいほど圧縮上死点より遅れると考え、また、燃焼室内ガス質量（負荷）によって影響されると判断し、さらにＦＴＡ（Ｆａｕｌｔ Ｔｒｅｅ Ａｎａｌｙｓｉｓ）より燃焼室内のガス流速（タンブル、スワール）の影響をもつと判断し、次の（８）式により燃焼圧が最大となるクランク角θｐｍａｘ［ｄｅｇＡＴＤＣ］を算出することとした。従って、（８）式と上記（補５）式とを比較すれば、ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｄを新たに導入したものである。

θｐｍａｘ＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｄ・Ｎｅ＾ｅ
…（８）
ただし、ｋ１ＶＥＬ：流速係数［無名数］、
ＩＶＣ ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θｔ０ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｍ ：燃焼室内ガス質量［ｇ］、
Ｎｅ ：エンジン回転速度［ｒｐｍ］、
ｄ、ｅ ：適合定数［無名数］、
（８）式に示す燃焼圧が最大となるクランク角θｐｍａｘに対するこうした新しい考え方が妥当かどうかを確認するため、吸気弁１５のバルブリフト量及び吸気弁閉時期ＩＶＣを一定として、燃焼圧が最大となるクランク角θｐｍａｘと吸入空気量（燃焼室内ガス質量ｍの代表値）、エンジン回転速度Ｎｅ（ガス流速の代表値）の関係を実験したところ、図９に示す結果が得られた。図９によれば、吸気弁１５のバルブリフト量及び吸気弁閉時期ＩＶＣが一定の条件で、吸入空気量が同じであればエンジン回転速度Ｎｅが大きいときのほうがエンジン回転速度Ｎｅが小さいときより燃焼圧が最大となるクランク角θｐｍａｘが大きくなる（圧縮上死点よりの遅角量が大きくなる）こと、またエンジン回転速度Ｎｅが同じであれば吸入空気量が多くなるほど燃焼圧が最大となるクランク角θｐｍａｘが大きくなる（圧縮上死点よりの遅角量が大きくなる）ことが確かめられた。

なお、図９の縦軸の単位［ｄｅｇＡＴＤＣ］は圧縮上死点を起点とするクランク角であり、吸気上死点を起点とするものでない。

上記（８）式のタンブル形成開始角θ０ｔについては上記（６）式のところで前述したところと同じである。すなわち、タンブル形成開始角θ０ｔは一定値（例えば４５ｄｅｇＡＴＤＣ）である。タンブル形成開始角θ０ｔはエンジンの仕様に依存し、ＶＥＬ機構２８、ＶＴＣ機構２７の有無には依存しない。図２より、吸気弁閉時期ＩＶＣはタンブル形成開始角θ０ｔより必ず遅角側の値である。従って、（８）式の（ＩＶＣ−θ０ｔ）の値は必ず正の値になる。

上記（８）式の流速係数ｋ１ＶＥＬは吸気弁１５のバルブリフト量Ｌｉｆｔ（あるいはバルブ作動角）に応じて定めている。ここでも簡単のため、前述したようにＶＥＬ機構２８の作動、非作動により吸気弁１５のバルブリフト量が大小の２段に切換可能であり、運転条件１ではＶＥＬ機構２８を非作動状態とし、運転条件１より運転条件２に移行したときにはＶＥＬ機構２８を作動させ、図１０下段に示したように実線のバルブリフト特性から破線のバルブリフト特性へとバルブリフト量を小さくしているとして説明すると、図１０上段に示したように、ＶＥＬ機構２８の非作動時（あるいは運転条件１のとき）には第１バルブリフト量Ｌｉｆｔ１に対応する第１流速係数ｋ１ＶＥＬ１を、これに対してＶＥＬ機構２８の作動時（あるいは運転条件２のとき）には第２バルブリフト量Ｌｉｆｔ２に対応する第２流速係数ｋ１ＶＥＬ２を流速係数ｋ１ＶＥＬとして設定する。このように、バルブリフト量が小さくなったときに流速係数ｋ１ＶＥＬを第１流速係数ｋ１ＶＥＬ１より第２流速係数ｋ１ＶＥＬ２へと小さくするのは、バルブリフト量が小さいときのほうがタンブルが弱くなる（従って、θｐｍａｘがより遅角側に移動する）からである。

また、上記（８）式のｍ＾ｄはθｐｍａｘに対する筒内ガス質量ｍの、Ｎｅ＾ｅはθｐｍａｘに対する回転速度Ｎｅの各影響を反映させている。

このように、（８）式によれば、燃焼圧が最大となるクランク角θｐｍａｘを、エンジン回転速度Ｎｅに加えて、燃焼室内ガス質量ｍ、吸気弁１５のバルブリフト量及び吸気弁閉時期ＩＶＣの関数でもあるとして新たに構成している。
〈３〉着火遅れ時間の算出方法
コンベンショナルエンジンに適用している先行点火時期制御方法では、点火タイミングを起点として燃焼質量割合が０％となるときのクランク角位置θｍｂ０％である間の区間を着火遅れ時間ＤＥＡＤＴＩＭＥとして算出している。

しかしながら、燃焼質量割合が０％となるときのクランク角位置であるθｍｂ０％は冷却損失の算出誤差やガス漏れの影響を受けやすく、正確に計測することが難しい。そのため、このθｍｂ０％の値を正として着火遅れ時間ＤＥＡＤＴＩＭＥの適合を行うと燃焼速度の推定誤差に着火遅れ時間の計測誤差が含まれてしまう。また、このような理由により、一般的にもθｍｂ０％を正にして燃焼シミュレーションの適合を行うこともまれである。

これに対して本発明の参考例では、点火タイミングを起点として燃焼質量割合が２％となるときのクランク角位置であるθｍｂ２％までの時間を着火遅れ時間τとみなす。このため、本発明の参考例ではθｍｂ２％よりθｍｂ１０％までのクランク角区間が燃焼期間となるが、これはコンベンショナルエンジンに適用している先行点火時期制御方法では主燃焼期間ＢＵＲＮ２に相当しているので、本発明の参考例では、燃焼期間の推定を主燃焼期間のみとすることとなる。

また、この着火遅れ時間τは純粋な意味での着火遅れ時間ではなく、点火タイミングから燃料が２％燃えるまでの期間である。そのため、コンベンショナルエンジンに適用している先行点火時期制御方法において提案しているアレニウスの式（自己着火の式）を用いることができないため、本発明の参考例では、着火遅れ時間τ［ｓｅｃ］を新たに次の（９）式により算出することとする。この（９）式は点火直後の燃焼速度を算出する式を用いて算出するものである。

τ＝（Ｄkernel−Ｄ０）
／２［（Ｔａｄ／Ｔ）・ＳＬ＋｛（２／３）・ｋ｝＾（１／２）］
…（９）
ただし、Ｄkernel：θｍｂ２％時の火炎直径［ｍ］、
Ｄ０ ：点火タイミング直後の火炎直径［ｍ］、
Ｔａｄ ：火炎温度［Ｋ］、
Ｔ ：未燃ガス温度［Ｋ］、
ＳＬ ：層流燃焼速度［ｍ／ｓｅｃ］、
ｋ ：燃焼室内ガスの運動エネルギー、
（９）式の火炎直径Ｄkernelは燃焼質量割合が２％となったときの火炎の直径であるので、次の（１０）式により新たに算出する。

Ｄkernel＝｛（Ｖｃｙｌ・［１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）］
・６）／π｝＾（１／３） …（１０）
ただし、Ｖｃｙｌ：燃焼室内容積［ｍ³］、
ｘｂ ：燃焼質量割合（＝０．０２）、
κ ：ポリトロープ指数（＝１．３４）、
（１０）式は次のようにして導いたものである。火炎の形状を球と仮定すると、火炎直径Ｄkernelを与える式は次の（１１ａ）式である。

Ｄkernel＝｛（Ｖｂ・６）／π｝＾（１／３） …（１１ａ）
ただし、Ｖｂ：火炎体積、
熱力学の公式より体積燃焼割合の式は次の（１１ｂ）式で表される（公知）。

体積燃焼割合＝１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）
…（１１ｂ）
この（１１ｂ）式と、
Ｖｂ＝Ｖｃｙｌ・体積燃焼割合 …（１１ｃ）
の関係式とを利用すれば次の（１１ｄ）式が得られる。

Ｖｂ＝Ｖｃｙｌ・［１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）］
…（１１ｄ）
この（１１ｄ）式を（１１ａ）式に代入すれば上記（１０）式が得られる。

上記（９）式の点火タイミング直後の火炎の直径Ｄ０は１ｍｍとする。火炎温度Ｔａｄと未燃ガス温度Ｔの比であるＴａｄ／Ｔは適合定数とする。層流燃焼速度ＳＬの算出方法はコンベンショナルエンジンに適用している先行点火時期制御方法と同じとする。

上記（９）式の燃焼室内ガスの運動エネルギーｋは次式により算出する。

ｋ＝（１／２）ｍ（ｆ・Ｎｅ）＾２ …（１２）
ただし、ｍ ：燃焼室内ガス質量［ｇ］、
Ｎｅ：エンジン回転速度［ｒｐｍ］、
ｆ ：適合定数、
これで、項分け説明を終了する。

次に、上記〈１〉〜〈３〉の変更を行った後の新しい点火時期制御方法（この新しい点火時期制御方法を以下「変更後点火時期制御方法」という。）を、以下にまとめる。ただし、適合係数については改めて取り直している。

ＭＢＴの得られる点火時期（この点火時期を「基本点火時期」という。）ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を次式により算出する。

ＭＢＴＣＡＬ＝−｜（τ＋ＢＴ）・Ｎｅ・６−θｍｂ１０％｜…（１３）
ただし、ＢＴ ：燃焼時間［ｓｅｃ］、
τ ：着火遅れ時間［ｓｅｃ］、
θｍｂ１０％：基準クランク角［ｄｅｇＡＴＤＣ］、
これを図解したのが図１１である。変更後点火時期制御方法では、点火タイミングからθｍｂ２％までの時間である着火遅れ時間τ［ｓｅｃ］と、θｍｂ２％からθｍｂ１０％までの時間である燃焼時間ＢＴ［ｓｅｃ］とを加算し、この加算値［ｓｅｃ］にＮｅ［ｒｐｍ］・６を乗算することによってクランク角区間［ｄｅｇ］に変換する。この換算されたクランク角区間（（τ＋ＢＴ）・Ｎｅ・６）は点火タイミングからθｍｂ１０％（＝θｐｍａｘ）までのクランク角区間である。従って、基準クランク角（＝θｐｍａｘ）であるθｍｂ１０％［ｄｅｇＡＴＤＣ］よりこの換算されたクランク角区間だけ進角側の値を基本点火時期ＭＢＴＣＡＬとして算出する。

なお、（１３）式において絶対値をとりマイナスの符号を付けているのは、θｍｂ１０％の単位は圧縮上死点より遅角側に計測するクランク角［ｄｅｇＡＴＤＣ］であるのに対して、ＭＢＴＣＡＬの単位は圧縮上死点より進角側に計測するクランク角［ｄｅｇＢＴＤＣ］であるため、θｍｂ１０％の単位をＭＢＴＣＡＬの単位へと変換してやる必要があるためである。

ここで、基準クランク角であるθｍｂ１０％［ｄｅｇＡＴＤＣ］は、次の（１４）式により算出する。

θｍｂ１０％＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｃ１・Ｎｅ＾ｃ２＋ｃ３
…（１４）
ただし、ｋ１ＶＥＬ：流速係数［無名数］、
ＩＶＣ ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θｔ０ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｍ ：燃焼室内ガス質量［ｇ］、
Ｎｅ ：エンジン回転速度［ｒｐｍ］、
ｃ１〜ｃ３：適合係数、
前述したように、吸気弁閉時期ＩＶＣとタンブル形成開始角θｔ０の単位は吸気上死点を起点としており、一方、θｍｂ１０％は圧縮上死点を起点としている。従って、吸気上死点を起点とする単位より圧縮上死点を起点とする単位への変換は適合係数ｃ３により行うこととなる。例えば４気筒エンジンでは吸気上死点と圧縮上死点の間に１８０ｄｅｇのずれがあるので、適合係数ｃ３に１８０ｄｅｇを入れてやればよい。

（１３）式の燃焼時間ＢＴ［ｓｅｃ］は次式により算出する。

ＢＴ＝０．０８・Ｖｃｙｌ／｛（（３・ｘｂ＋１）＾（１／κ））・Ａｂ・Ｓｂ｝
＝０．０７０７・Ｖｃｙｌ／（Ａｂ・Ｓｂ） …（１５）
ただし、Ｖｃｙｌ：燃焼室内容積［ｍ³］、
ｘｂ ：燃焼時間算出時の平均燃焼質量割合（＝０．０６）、
κ ：ポリトロープ指数、
Ａｂ ：火炎表面積［ｍ²］、
Ｓｂ ：燃焼速度［ｍ／ｓｅｃ］、
（１５）式はコンベンショナルエンジンに適用している先行点火時期制御方法における式と基本的に同じ式である。すなわち、（１５）式においてＶｃｙｌ／｛（３・ｘｂ＋１）＾（１／κ）｝は既燃ガスの質量割合であり、燃焼時間ＢＴはこの既燃ガスの質量割合に比例し、燃焼速度Ｓｂに反比例するとする式である。（１５）式の火炎表面積Ａｂはエンジンの仕様から算出することができる。

なお、（１５）式の燃焼質量割合ｘｂとしては、θｍｂ２％からθｍｂ１０％までのクランク角区間（燃焼期間）における値を用いる必要がある。この場合に、燃焼質量割合はθｍｂ２％からθｍｂ１０％までのクランク角区間で２％から１０％へと変化するので、ここでは、θｍｂ２％とθｍｂ１０％とを平均した値であるθｍｂ６％時の燃焼質量割合の値つまり６％を用いる。

上記（１５）式の燃焼速度Ｓｂ［ｍ／ｓｅｃ］は次式により算出する。

Ｓｂ＝ＳＬ＋［｛ｋＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｃ４・Ｎｅ＾ｃ５｝
／ＳＬ］＾ｃ６ …（１６）
ただし、ＳＬ ：層流燃焼速度［ｍ／ｓｅｃ］、
ｋＶＥＬ ：タンブル強度のバルブリフト量補正係数［１／ｄｅｇ］、
ＩＶＣ ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θ０ｔ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｃ４、ｃ５：適合係数、
ｃ６ ：適合係数、
（１６）式右辺第２項が乱流燃焼速度［ｍ／ｓｅｃ］を与える式であり、乱流燃焼速度を（１６）式右辺第２項で近似している点が新しい。つまり、本発明の参考例では、燃焼室内ガスの乱れ強さｕ［ｍ／ｓｅｃ］を次の（１６−１）式により算出し、この燃焼室内ガスの乱れ強さｕに基づいて次の（１６−２）式により乱流燃焼速度ＳＴ［ｍ／ｓｅｃ］を算出している。

ｕ＝ｋＶＥＬ・（ＩＶＣ−θｔ０）・ｍ＾ｃ４・Ｎｅ＾ｃ５ …（１６−１）
ＳＴ＝［ｕ／ＳＬ］＾ｃ６ …（１６−２）
これに対して、（１６）式の層流燃焼速度ＳＬ［ｍ／ｓｅｃ］としては、コンベンショナルエンジンに適用している先行点火時期制御方法と同じでよく、従って次の（１７）式により算出する。

ＳＬ＝ＳＬ０・（Ｔ／２９８）＾ａ・（ｐ／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ＾ｃ） …（１７）
ただし、ＳＬ０ ：標準状態での層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ ：未燃ガス温度［Ｋ］、
ｐ ：燃焼室内圧力［ｋＰａ］、
ＭＲＥＳＲ：内部不活性ガス率［％］、
ａ、ｂ、ｃ：係数、
一方、上記（１３）式の着火遅れ時間τ［ｓｅｃ］は次の（１８）式により新たに算出する。

τ＝［｛｛Ｖｃｙｌ・［１−（１−ｘｂ）／（３・ｘｂ＋１）＾（１／κ）］・６｝ ／π｝＾（１／３）−０．００１］
／｛ｃ７／ＳＬ＋ｃ８・Ｎｅ・ｍ＾（１／２）｝
＝｛０．４９０・Ｖｃｙｌ＾（１／３）｝
／｛ｃ７・ＳＬ＋ｃ８・Ｎｅ・ｍ＾（１／２）｝ …（１８）
ただし、ｘｂ ：着火遅れ時間算出時の燃焼質量割合（＝０．０２）、
κ ：ポリトロープ指数（＝１．３４）、
ＳＬ ：層流燃焼速度［ｍ／ｓｅｃ］、
ｍ ：燃焼室内ガス質量［ｇ］、
Ｎｅ ：エンジン回転速度［ｒｐｍ］、
Ｖｃｙｌ ：燃焼室内容積［ｍ³］、
ｃ７、ｃ８：適合定数、
なお、（１８）式の燃焼質量割合ｘｂとしては、θｍｂ０％からθｍｂ２％までのクランク角区間（着火遅れ期間）における値を用いる必要がある。この場合に、燃焼質量割合はθｍｂ０％からθｍｂ２％までのクランク角区間で０％より２％まで変化するのであるが、ここでは、θｍｂ０％とθｍｂ２％とを平均した値であるθｍｂ１％時の燃焼質量割合の値を用いるのではなく、着火遅れ期間の終期であるθｍｂ２％時の燃焼質量割合の値つまり２％を用いることとする。

上記（１８）式の層流燃焼速度ＳＬ［ｍ／ｓｅｃ］としては、次の（１９）式により算出する。

ＳＬ＝ＳＬ０・（Ｔ／２９８）＾ａ・（ｐ／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ＾ｃ） …（１９）
ただし、ＳＬ０ ：標準状態での層流燃焼速度［ｍ／ｓｅｃ］、
Ｔ ：未燃ガス温度［Ｋ］、
ｐ ：燃焼室内圧力［ｋＰａ］、
ＭＲＥＳＲ：内部不活性ガス率［％］、
ａ、ｂ、ｃ：係数、
ここで、着火遅れ時間τを算出するのに用いる（１９）式の層流燃焼速度ＳＬと、燃焼時間ＢＴを算出するのに用いる上記（１７）式の層流燃焼速度ＳＬとでは、式そのものは変わらず、また標準状態での層流燃焼速度ＳＬ０、内部不活性ガス率ＭＲＥＳＲの値も変わらないが、未燃ガス温度Ｔ、燃焼室内圧力ｐとして代入する値が（１９）式と（１７）式とで、後述するように相違することとなる。

上記（１７）式及び（１９）式の標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂ、ｃ［無名数］は次の（２０）式〜（２３）式により算出する。

ＳＬ０＝（０．２６３２−０．８４７２／（φ−１．１３）＾２）
…（２０）
ａ＝２．１８−０．８０・（φ−１） …（２１）
ｂ＝−０．１６＋０．２２・（φ−１） …（２２）
ｃ＝１ …（２３）
ただし、φ：当量比［無名数］、
ここで、（２０）式〜（２３）式についてはＳＡＥペーパーにより公知である（ＳＡＥ１９９９１０１７５参照）。

これで変更後点火時期制御方法をまとめたものの説明を終了する。

次に、上記（１１ｂ）式に示した体積燃焼割合の求め方と、上記（１７）式、（１９）式に示した未燃ガス温度の算出方法とについて補足説明を行う。

まず、上記（１１ｂ）式に示した体積燃焼割合の求め方を簡単に説明する。

既燃ガスの質量燃焼割合をｘ、体積燃焼割合をｙとして、燃焼室内の燃焼が定容燃焼同様に、
ｙ＝ｆ（ｘ）
＝１＋（ｘ−１）／（１＋ｘ（ｋ１−１）＾（１／κ））
…（２４）
ただし、ｋ１≒４〜５、
の式により書けると仮定すると、（２４）式のｋ１に４を代入して計算することにより、上記（１１ｂ）式が容易に得られる。

上記（１７）式、（１９）式に示した未燃ガス温度Ｔの算出方法をまとめて説明すると、断熱変化と仮定し、熱力学の次の（２５）式を用いて燃焼室５のθｍｂ ｘ％時の未燃ガス温度Ｔｍｂ ｘ％を算出する。

Ｔｍｂ ｘ％＝Ｔｉｖｃ・（Ｐｍｂ ｘ％／Ｐｉｖｃ）＾｛（κ−１）／κ｝ …（２５）
ただし、Ｔｍｂ ｘ％：θｍｂ ｘ％時の未燃ガス温度［Ｋ］、
Ｔｉｖｃ ：ＩＶＣ時の燃焼室内温度［Ｋ］、
Ｐｍｂ ｘ％：θｍｂ ｘ％時の燃焼室内圧力［ｋＰａ］、
Ｐｉｖｃ ：ＩＶＣ時の燃焼室内圧力［ｋＰａ］、
κ ：比熱比（固定値で１．３〜１．４）、
ここで、（２５）式の燃焼室５のθｍｂ ｘ％時の未燃ガス圧力Ｐｍｂ ｘ％の算出方法を説明すると、この燃焼室５のθｍｂ ｘ％時の未燃ガス圧力Ｐｍｂ ｘ％は次の（２６）式により与えられる。

Ｐｍｂ ｘ％＝Ｐｉｖｃ＋Δｐ …（２６）
ただし、Δｐ：燃焼室内未燃ガス圧力の変化量、
ここで、燃焼室内未燃ガス圧力の変化量Δｐは、容積変化による圧力変化分Δｐｖと、燃焼による圧力変化分Δｐｃとに分けて算出する。つまり、次の（２７）式により燃焼室内未燃ガス圧力の変化量Δｐを算出する。

Δｐ＝Δｐｖ＋Δｐｃ …（２７）
このうち、（２７）式の容積変化による圧力変化分Δｐｖは熱力学の式である次の（２８）式により算出することができる。

Δｐｖ＝Ｐｉｖｃ・｛（Ｖｉｖｃ／Ｖｍｂ ｘ％）＾κ−１｝ …（２８）
ただし、Ｖｉｖｃ ：ＩＶＣ時の燃焼室内容積、
Ｖｍｂ ｘ％：θｍｂ ｘ％時の燃焼室内容積、
ここで、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃはピストン６の位置から求めることができる。燃焼室５のθｍｂ ｘ％時における未燃ガス容積Ｖｍｂ ｘ％については、図１２に示したように、θｍｂ ｘ％とＶｍｂ ｘ％の関係を表す特性を予め作成しておき、θｍｂ ｘ％から当該特性を用いて、燃焼室５のθｍｂ ｘ％時における未燃ガス容積Ｖｍｂ ｘ％を求めればよい。図１２は要するにクランク角に対する燃焼室容積の特性である。

一方、上記（２７）式の燃焼による圧力変化分Δｐｃは、次の（２９）式により簡易に算出すればよい。

Δｐｃ＝（ｘ／１００）×Δｐｔｏｔａｌ …（２９）
ただし、Δｐｔｏｔａｌ：燃焼によるトータルの圧力上昇分、
ここで、燃焼によるトータルの圧力上昇分Δｐｔｏｔａｌについては、図１３に示したように、燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］（または燃料噴射量）とΔｐｔｏｔａｌの関係を表す特性を予め作成しておき、燃料噴射パルス幅Ｔｉ（または燃料噴射量）から当該特性を用いて燃焼によるトータルの圧力上昇分Δｐｔｏｔａｌを求める。

これで補足説明を終える。

次に、エンジンコントローラ３１で実行される上記変更後点火時期制御方法における基本点火時期ＭＢＴＣＡＬの算出方法を図１４、図１６、図１７、図１８、図１９のフローチャートを参照しながら詳述する。

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

まずステップ１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、バルブリフト量Ｌｉｆｔ［ｍ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、内部不活性ガス率ＭＲＥＳＲ［％］、燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］、クランク角センサ（３３、３４）により検出されるエンジン回転速度Ｎｅ［ｒｐｍ］、総ガス質量ＭＧＡＳ［ｇ］を読み込む。

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

吸気弁１５のバルブリフト量ＬｉｆｔはＶＥＬ機構２８に与える指令値から既知である。つまり、図４下段、図１０下段に示したように、ＶＥＬ機構２８の非作動時には第１バルブリフト量Ｌｉｆｔ１が、これに対してＶＥＬ機構２８の作動時には第２バルブリフト量Ｌｉｆｔ２がバルブリフト量Ｌｉｆｔである。

燃料噴射パルス幅Ｔｉは図示しない燃料噴射パルス幅の算出ルーチンにおいて算出されている。例えば、シーケンシャル噴射時に燃料インジェクタ２１に与える燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］の算出式は次のようなものである。

Ｔｉ＝Ｔｐ×Ｔｆｂｙａ×（α＋αｍ−１）×２＋Ｔｓ …（３０）
ただし、Ｔｐ ：基本噴射パルス幅［ｍｓｅｃ］、
Ｔｆｂｙａ：［無名数］、
α ：空燃比フィードバック補正係数［無名数］、
αｍ ：空燃比学習値［無名数］、
Ｔｓ ：無効パルス幅［ｍｓｅｃ］、
（３０）式の目標当量比Ｔｆｂｙａは無名数であり、理論空燃比を１４．７とすると、次の（３１）式により表される値である。

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

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

内部不活性ガス率ＭＲＥＳＦＲは燃焼室５内に残留する不活性ガス量を燃焼室５内の総ガス量で除した値で、その算出については総ガス質量ＭＧＡＳと共に図１６により後述する。

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

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

Ｖｉｖｃ＝ｆ１（θ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（＝ＩＶＣ）を代入すれば、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃを算出することができる。したがって、実用上は燃焼室５の吸気弁閉時期における容積Ｖｉｖｃは吸気弁閉時期ＩＶＣをパラメータとするテーブルで設定したものを用いる。

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

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

ステップ５では、燃料噴射パルス幅Ｔｉ（または燃料噴射量）から図１３を内容とするテーブルを検索することにより、トータルの燃焼による圧力上昇分Δｐｔｏｔａｌ［Ｐａ］を算出する。

ステップ６では、バルブリフト量Ｌｉｆｔから図４上段を内容とするテーブルを検索することにより、タンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］を算出する。すなわち、ＶＥＬ機構２８の非作動時には第１バルブリフト量補正係数ｋＶＥＬ１を、これに対してＶＥＬ機構２８の作動時には第２バルブリフト量補正係数ｋＶＥＬ２を、バルブリフト量補正係数ｋＶＥＬとして算出する。

ステップ７では、同じくバルブリフト量Ｌｉｆｔから図１０上段を内容とするテーブルを検索することにより、流速係数ｋ１ＶＥＬ［無名数］を算出する。すなわち、ＶＥＬ機構２８の非作動時には第１流速係数ｋ１ＶＥＬ１を、これに対してＶＥＬ機構２８の作動時には第２流速係数ｋ１ＶＥＬ２を、流速係数ｋ１ＶＥＬとして算出する。

ステップ８では、このうちの流速係数ｋ１ＶＥＬ、吸気弁閉時期ＩＶＣ、総ガス質量ＭＧＡＳ、エンジン回転速度Ｎｅを用いて、次の（４７）式により燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％［ｄｅｇＡＴＤＣ］を算出する。

θｍｂ１０％＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ１・Ｎｅ＾ｃ２
＋ｃ３ …（４７）
ただし、ｋ１ＶＥＬ：流速係数［無名数］、
ＩＶＣ ：吸気弁閉時期［ｄｅｇＡＴＤＣ］、
θｔ０ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ＭＧＡＳ ：総ガス質量［ｇ］、
Ｎｅ ：エンジン回転速度［ｒｐｍ］、
ｃ１〜ｃ３：適合係数、
（４７）式は、上記（１４）式において、燃焼室内ガス質量ｍとして総ガス質量ＭＧＡＳを用いたもので、基本的に上記（１４）式と変わらない。

ここで、燃焼圧が最大となるクランク角θｐｍａｘは（４７）式右辺第１項、つまり次の（４７−１）式により与えていることになる。

θｐｍａｘ＝ｋ１ＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ１・Ｎｅ＾ｃ２
…（４７−１）
ステップ９では燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％から単純にクランク角で４ｄｅｇを差し引いた値を燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％［ｄｅｇＡＴＤＣ］として求める。これはθｍｂ２％よりθｍｂ１０％までのクランク角区間で燃焼質量割合が直線的に変化するとみなしてθｍｂ６％を求めるようにしたものある。

ステップ１０では、このθｍｂ６％から図１２を内容とするテーブルを検索することにより、燃焼室５のθｍｂ６％時における容積Ｖｍｂ６％を算出する。

ステップ１１、１２はステップ９、１０と同様である。 ステップ１１では燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％から単純にクランク角で８ｄｅｇを差し引いた値を燃焼質量割合が２％となるときのクランク角位置であるθｍｂ２％［ｄｅｇＡＴＤＣ］として求める。ステップ１２ではこのθｍｂ２％から図１２を内容とするテーブルを検索することにより、燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％を算出する。

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

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

ステップ２３では、燃焼室５内の内部不活性ガス量ＭＲＥＳ［ｇ］を算出する。この内部不活性ガス量ＭＲＥＳの算出についても公知の方法を用いればよい（特開２００５−１７１８５６号公報参照）。

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

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

ステップ２５では、この総ガス質量ＭＧＡＳ、内部不活性ガス量ＭＲＥＳを用いて、次の（４９）式により内部不活性ガス率ＭＲＥＳＲ（燃焼室５内の総ガス質量に対する内部不活性ガス量の割合）［％］を算出する。

ＭＲＥＳＲ＝ＭＲＥＳ／ＭＧＡＳ …（４９）
図１７は燃焼時間ＢＴ［ｓｅｃ］を算出するためのもの、図１８は着火遅れ時間τ［ｓｅｃ］を算出するためのもので、それぞれ一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１７、図１８は図１４に続けて実行する。図１７、図１８はどちらを先に実行してもかまわない。

ここでは図１７から先に説明すると、ステップ３１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、図１４のステップ２で算出されている燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ［ｍ³］、図１４のステップ３で算出されている燃焼室５の吸気弁閉時期における温度Ｔｉｖｃ［Ｋ］、図１４のステップ４で算出されている燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ［Ｐａ］、図１４のステップ５で算出されているトータルの燃焼による圧力上昇分Δｐｃｔｏｔａｌ［ｋＰａ］、図１４のステップ６で算出されているタンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］、図１４のステップ１０で算出されている燃焼室５のθｍｂ６％時における容積Ｖｍｂ６％［ｄｅｇＡＴＤＣ］、図１６のステップ２４で算出されている総ガス質量ＭＧＡＳ［ｇ］、図１６のステップ２５で算出されている内部不活性ガス率ＭＲＥＳＲ［％］、燃料噴射パルス幅Ｔｉ［ｍｓｅｃ］、目標当量比Ｔｆｂｙａ［無名数］、エンジン回転速度Ｎｅ［ｒｐｍ］を読み込む。

ステップ３２では、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ、燃焼室５のθｍｂ６％時における容積Ｖｍｂ６％から、次の（５０）式により、容積変化による圧力変化分Δｐｖ［ｋＰ］を算出する。

Δｐｖ＝Ｐｉｖｃ・｛（Ｖｉｖｃ／Ｖｍｂ６％）＾κ−１｝ …（５０）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（５０）式は上記（２８）式と同じものである。ここでの容積変化による圧力上昇分Δｐｖはθｍｂ２％よりθｍｂ１０％までのクランク角区間（燃焼期間）における圧力上昇分であるので、θｍｂ２％とθｍｂ１０％の平均値であるθｍｂ６％を用いている。

ステップ３３では、燃焼によるトータルの圧力上昇分Δｐｃｔｏｔａｌから、次の（５１）式により、燃焼による圧力上昇分Δｐｃ［ｋＰ］を算出する。

Δｐｃ＝０．０６×Δｐｔｏｔａｌ …（５１）
（５１）式は上記（２９）式において、燃焼質量割合ｘとして６％を用いたものである。

ステップ３４では、これら燃焼による圧力上昇分Δｐｃ、容積変化による圧力上昇分Δｐｖを燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃに加算して、つまり次の（５２）式により燃焼室５のθｍｂ６％時における未燃ガス圧力Ｐｍｂ６％［ｋＰ］を算出する。

Ｐｍｂ６％＝Ｐｉｖｃ＋Δｐｖ＋Δｐｃ …（５２）
ステップ３５では、この燃焼室５のθｍｂ６％時における未燃ガス圧力Ｐｍｂ６％、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における温度Ｔｉｖｃを用いて、次の（５３）式により燃焼室５のθｍｂ６％時における未燃ガス温度Ｔｍｂ６％［Ｋ］を算出する。

Ｔｍｂ６％＝Ｔｉｖｃ・（Ｐｍｂ６％／Ｐｉｖｃ）＾｛（κ−１）／κ｝ …（５３）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（５３）式は上記（２５）式において燃焼質量割合ｘが６％のときの値である。

ステップ３６では、目標当量比Ｔｆｂｙａから次の（５４）式〜（５６）式により標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂ［無名数］を求め、ステップ３７で、これら標準状態での層流燃焼速度ＳＬ０、係数ａ、ｂと、燃焼室５のθｍｂ６％時における未燃ガス温度Ｔｍｂ６％、燃焼室５のθｍｂ６％時における未燃ガス圧力Ｐｍｂ６％、内部不活性ガス率ＭＲＥＳＲを用いて、次の（５７）式によりθｍｂ６％時の層流燃焼速度ＳＬ［ｍ／ｓ］を算出する。

ＳＬ０＝（０．２６３２−０．８４７２／（Ｔｆｂｙａ−１．１３）＾２）
…（５４）
ａ＝２．１８−０．８０・（Ｔｆｂｙａ−１） …（５５）
ｂ＝−０．１６＋０．２２・（Ｔｆｂｙａ−１） …（５６）
ＳＬ＝ＳＬ０・（Ｔｍｂ６％／２９８）＾ａ・（Ｐｍｂ６％／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ） …（５７）
（５４）式〜（５６）式は、上記（２０）式〜（２２）式において、当量比φとして目標当量比Ｔｆｂｙａを用いたもので、基本的に（２０）式〜（２２）式と変わらない。同様に、（５７）式は、上記（１７）式において、未燃ガス温度Ｔ、燃焼室内圧力ｐとしてＴｍｂ６％、Ｐｍｂ６％を用いたもので、基本的に（１７）式と変わらない。

ステップ３８では、このθｍｂ６％時の層流燃焼速度ＳＬ、タンブル強度のバルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣ、総ガス質量ＭＧＡＳ、エンジン回転速度Ｎｅを用いて、次の（５８）式により燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％時の燃焼速度Ｓｂ［ｍ／ｓ］を算出する。

Ｓｂ＝ＳＬ＋［｛ｋＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ４・Ｎｅ＾ｃ５｝
／ＳＬ］＾ｃ６ …（５８）
ただし、θｔ０ ：タンブル形成開始角［ｄｅｇＡＴＤＣ］、
ｃ４、ｃ５、ｃ６：適合係数、
（５８）式は、上記（１６）式において、燃焼室内ガス質量ｍとして総ガス質量ＭＧＡＳを用いたもので、基本的に（１６）式と変わらない。

ここで、（５８）式は詳細には次の４つの式を一つにしたものである。

Ｉｔ＝ｋＶＥＬ・（ＩＶＣ−θｔ０）・ＭＧＡＳ＾ｃ４・Ｎｅ＾（ｃ５−１）
…（５８−１）
ｕ＝Ｉｔ・Ｎｅ …（５８−２）
ＳＴ＝［ｕ／ＳＬ］＾ｃ６ …（５８−３）
Ｓｂ＝ＳＬ＋ＳＴ …（５８−４）
ただし、Ｉｔ：タンブル強度［無名数］、
ｕ：燃焼室内ガスの乱れ強さ［ｍ／ｓｅｃ］、
ＳＴ：乱流燃焼速度［ｍ／ｓｅｃ］、
ＶＴＣ機構２７の作動で吸気弁閉時期ＩＶＣが遅角側に遅らされるほど、つまりタンブル開始角を起点としてＩＶＣまでのクランク角区間（ＩＶＣ−θｔ０）が大きくなるほど実際のタンブル強度が大きくなるが、（５８−１）式よれば、吸気弁閉時期ＩＶＣが遅れるとき（ＩＶＣの値が大きくなるとき）タンブル強度Ｉｔも大きくなるのであり、算出値としてのＩｔは実際のタンブル強度と良く一致することとなる。

また、ＶＥＬ機構２８の作動で吸気弁のバルブリフト量が大きくなるほど実際のタンブル強度が大きくなるが、図４上段によればバルブリフト量Ｌｉｆｔが大きくなるときバルブリフト量補正係数ｋＶＥＬが大きくなり、（５８−１）式よれば、バルブリフト量補正係数ｋＶＥＬが大きくなるときタンブル強度Ｉｔも大きくなるのであり、算出値としてのＩｔは実際のタンブル強度と良く一致することとなる。

さらに、（５８−１）式よれば、総ガス質量ＭＧＡＳが大きくなるほどタンブル強度Ｉｔが大きくなる。

（５８−２）式は、タンブルが圧縮工程中に小渦に変換されることにより、燃焼室内ガスの乱れが強くなるので、タンブル強度Ｉｔが強いほど燃焼室内ガスの乱れが強くなると考え、タンブル強度Ｉｔ（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さｕを推定するものである。

ステップ３９では、この燃焼質量割合が６％となるときのクランク角位置であるθｍｂ６％時の燃焼速度Ｓｂ、燃焼室５のθｍｂ６％時におけ容積Ｖｍｂ６％を用いて、次の（５９）式により燃焼時間ＢＴ［ｓ］を算出する。

ＢＴ＝０．０７０７・Ｖｍｂ６％／（Ａｂ・Ｓｂ） …（５９）
ただし、Ａｂ：火炎表面積［ｍ²］、
（５９）式は、上記（１５）式において、燃焼室内容積ＶｃｙｌとしてＶｍｂ６％を用いたもので、基本的に（１５）式と変わらない。

図１８は着火遅れ時間τ［ｓｅｃ］を算出するためのもので、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。

ステップ４１〜４７は図１６のステップ３１〜３７と同様である。すなわち、ステップ４１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、図１４のステップ２で算出されている燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ［ｍ³］、図１４のステップ３で算出されている燃焼室５の吸気弁閉時期における温度Ｔｉｖｃ［Ｋ］、図１４のステップ４で算出されている燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ［Ｐａ］、図１４のステップ５で算出されているトータルの燃焼による圧力上昇分Δｐｃｔｏｔａｌ［Ｐａ］、図１４のステップ６で算出されているタンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］、図１４のステップ１２で算出されている燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％［ｄｅｇＡＴＤＣ］、図１６のステップ２４で算出されている総ガス質量ＭＧＡＳ［ｇ］、図１６のステップ２５で算出されている内部不活性ガス率ＭＲＥＳＲ［％］、燃料噴射パルス幅Ｔｉ［ｍｓ］、目標当量比Ｔｆｂｙａ［無名数］、エンジン回転速度Ｎｅ［ｒｐｍ］を読み込む。

ステップ４２では、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ、燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％から、次の（６０）式により、容積変化による圧力変化分Δｐｖ［ｋＰ］を算出する。

Δｐｖ＝Ｐｉｖｃ・｛（Ｖｉｖｃ／Ｖｍｂ２％）＾κ−１｝ …（６０）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（５９）式は上記（２８）式と同じものである。ここでの容積変化による圧力上昇分Δｐｖはθｍｂ０％よりθｍｂ２％までのクランク角区間（着火遅れ期間）における圧力上昇分であるが、θｍｂ０％とθｍｂ２％の平均値であるθｍｂ１％は用いず、着火遅れ期間の終期の値であるθｍｂ２％を用いている。

ステップ４３では、燃焼によるトータルの圧力上昇分Δｐｃｔｏｔａｌから、次の（６１）式により、燃焼による圧力上昇分Δｐｃ［ｋＰ］を算出する。

Δｐｃ＝０．０２×Δｐｔｏｔａｌ …（６１）
（６０）式は上記の（２９）式において、燃焼質量割合ｘとして２％を用いたものである。

ステップ４４では、これら燃焼による圧力上昇分Δｐｃ、容積変化による圧力上昇分Δｐｖを燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃに加算して、つまり次の（６２）式により燃焼室５のθｍｂ２％時における未燃ガス圧力Ｐｍｂ２％［ｋＰ］を算出する。

Ｐｍｂ２％＝Ｐｉｖｃ＋Δｐｖ＋Δｐｃ …（６２）
ステップ４５では、この燃焼室５のθｍｂ２％時における未燃ガス圧力Ｐｍｂ２％、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における温度Ｔｉｖｃを用いて、次の（６３）式により燃焼室５のθｍｂ２％時における未燃ガス温度Ｔｍｂ２％［Ｋ］を算出する。

Ｔｍｂ２％＝Ｔｉｖｃ・（Ｐｍｂ２％／Ｐｉｖｃ）＾｛（κ−１）／κ｝ …（６３）
ただし、κ：比熱比（固定値で１．３〜１．４）、
（６３）式は上記（２５）式において燃焼質量割合ｘが２％のときの値である。

ステップ４６では、目標当量比Ｔｆｂｙａから次の（６４）式〜（６６）式により標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂを求め、ステップ４７で、これら標準状態での層流燃焼速度ＳＬ０、係数ａ、ｂと、燃焼室５のθｍｂ２％時における未燃ガス温度Ｔｍｂ２％、燃焼室５のθｍｂ２％時における未燃ガス圧力Ｐｍｂ２％、内部不活性ガス率ＭＲＥＳＲを用いて、次の（６７）式によりθｍｂ２％時の層流燃焼速度ＳＬ［ｍ／ｓ］を算出する。

ＳＬ０＝（０．２６３２−０．８４７２／（Ｔｆｂｙａ−１．１３）＾２）
…（６４）
ａ＝２．１８−０．８０・（Ｔｆｂｙａ−１） …（６５）
ｂ＝−０．１６＋０．２２・（Ｔｆｂｙａ−１） …（６６）
ＳＬ＝ＳＬ０・（Ｔｍｂ２％／２９８）＾ａ・（Ｐｍｂ２％／１０１．３２５）＾ｂ
・（１−２．１×ＭＲＥＳＲ） …（６７）
（６４）式〜（６６）式は、上記（２０）式〜（２２）式において、当量比φとして目標当量比Ｔｆｂｙａを用いたもので、基本的に（２０）式〜（２２）式と変わらない。同様に、（６７）式は、上記（１９）式において、未燃ガス温度Ｔ、燃焼室内圧力ｐとしてＴｍｂ２％、Ｐｍｂ２％を用いたもので、基本的に（１９）式と変わらない。

ステップ４８では、このθｍｂ２％時の層流燃焼速度ＳＬ、燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％、エンジン回転速度Ｎｅ、総ガス質量ＭＧＡＳを用いて、次の（６８）式により着火遅れ時間τ［ｓｅｃ］を算出する。

τ＝｛０．４９０・Ｖｍｂ２％＾（１／３）｝
／｛ｃ７・ＳＬ＋ｃ８・Ｎｅ・ＭＧＡＳ＾（１／２）｝
…（６８）
ただし、ｃ７、ｃ８：適合係数、
（６８）式は、上記（１８）式において、燃焼室内容積Ｖｃｙｌ、燃焼室内ガス質量ｍとして、それぞれ燃焼室５のθｍｂ２％時における容積Ｖｍｂ２％、総ガス質量ＭＧＡＳを用いたもので、基本的に（１８）式と変わらない。

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

ステップ５１では、図１４のステップ８で算出されている燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％［ｄｅｇＡＴＤＣ］、図１７のステップ３９で算出されている燃焼時間ＢＴ［ｓｅｃ］、図１８のステップ４８で算出されている着火遅れ時間τ［ｓｅｃ］、エンジン回転速度Ｎｅ［ｒｐｍ］を読み込む。

ステップ５２では、着火遅れ時間τ、燃焼期間ＢＴ、燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％、エンジン回転速度Ｎｅを用いて、次の（６９）式により基本点火時期ＭＢＴＣＡＬ［ｄｅｇＡＴＤＣ］を算出し、ステップ５３で次の（７０）式によりその絶対値を採り、マイナスの符号を付けたものを改めて基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］として算出する。

ＭＢＴＣＡＬ＝（τ＋ＢＴ）・Ｎｅ・６−θｍｂ１０％ …（６９）
ＭＢＴＣＡＬ＝−｜ＭＢＴＣＡＬ｜ …（７０）
（６９）式、（７０）式は、上記（１３）式と同じ式である。

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

ここで、本発明の参考例の作用効果を説明する。

タンブル強度Ｉｔ（燃焼室内のガス流動）は、燃焼室内のさまざまな場所に分布する小渦の角運動量の総和で表されるので、例えばシリンダ（燃焼室５）の中心からｉ番目の小渦までの距離をｒ_i、ｉ番目の小渦のガス質量をｍ_i、ｉ番目の小渦のガス流速をｖ_iとすると、上記（３）式によりタンブル強度Ｉｔ（燃焼室内のガス流動）が表されるのであるが、燃焼室内の全ての小渦についてこれらの値ｍ_i、ｒ_i、ｖ_iをエンジンコントローラ３１においてオンラインで算出することは現在のところ不可能である。

そこで、本発明者は、タンブルが圧縮工程中に小渦に変換されることにより、燃焼室内ガスの乱れｕが強くなるため、タンブル強度Ｉｔが強いほど燃焼室内ガスの乱れが強くなると考え、タンブル強度Ｉｔ（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さｕを推定するものとした。すなわち、本発明の参考例によれば、タンブル強度Ｉｔ（燃焼室内のガス流動）に基づいて燃焼室内ガスの乱れ強さｕを算出し（上記（５８−２）式を参照）、この燃焼室内ガスの乱れ強さｕに基づいて乱流燃焼速度ＳＴを算出し（上記（５８−３）式を参照）、この乱流燃焼速度ＳＴに基づいてＭＢＴの得られる基本点火時期を算出する（上記（５８−４）式及び図１７のステップ３８、図１９のステップ５１〜５３参照）ので、タンブル強度Ｉｔ（燃焼室内のガス流動）の影響を受けて乱流燃焼速度ＳＴが変動しても、その変動する乱流燃焼速度ＳＴを精度よく算出できることから、基本点火時期ＭＢＴＣＡＬ（ＭＢＴの得られる基本点火時期）の予測精度を向上させることができる。

上記のように、燃焼室内の全ての小渦についてｍ_i、ｒ_i、ｖ_iをエンジンコントローラ３１においてオンラインで算出することは現在のところ不可能であるので、タンブル強度Ｉｔを求めるには何らかの近似が必要である。そこで、本発明者は燃焼室内のガス流動を新たに吸気弁１５のリフト特性とエンジン回転速度Ｎｅとに基づいて推定するものとした。ここで、吸気弁１５のリフト特性とは、吸気弁１５のバルブリフト量（あるいは吸気弁１５のバルブ作動角）または吸気弁１５の開閉時期のことである。

すなわち、本発明の参考例によれば、タンブル強度Ｉｔ（燃焼室内のガス流動）を吸気弁１５のバルブリフト量Ｌｉｆｔ及び吸気弁閉時期ＩＶＣ（リフト特性）とエンジン回転速度Ｎｅとに基づいて推定するので（上記（５８−１）式を参照）、簡易な構成でタンブル強度Ｉｔ（燃焼室内のガス流動）を精度よく算出できることになった。

ＶＥＬ機構２８を備えるときには、ＶＥＬ機構２８の作動状態で吸気弁１５のリフト量が小さいときと、ＶＥＬ機構２８の非作動状態で吸気弁１５のリフト量が大きいときでは燃焼室内のガス流動の向きが変わるため（図２（ｆ）と図２（ｇ）を参照）、燃焼室内ガスの乱れ強さｕも変わってしまうのであるが、本発明の参考例によれば、乱流燃焼速度ＳＴを少なくとも吸気弁１５のバルブリフト量Ｌｉｆｔ（リフト特性）に基づいて算出するので（図１４のステップ６、上記（５８−１）式〜（５８−３）式を参照）、ＶＥＬ機構２８の作動、非作動に伴う吸気弁１５のバルブリフト量Ｌｉｆｔの大小に関係なく乱流燃焼速度ＳＴの推定精度を向上できる。

ＶＴＣ機構２７の働きにより、吸気弁閉時期ＩＶＣが遅いほど燃焼室内ガスの乱れを形成する期間が長くなり、燃焼室内ガスの乱れ強さｕが強くなることに対応し、本発明の参考例によれば、乱流燃焼速度ＳＴを少なくとも吸気弁閉時期ＩＶＣに基づいて算出するので（上記（５８−１）式〜（５８−３）式を参照）、吸気弁閉時期ＩＶＣの進角、遅角に関係なく乱流燃焼速度ＳＴの推定精度を向上できる。

燃焼室５の総ガス質量ＭＧＡＳ（燃焼室内のガス量）が多いほどガスの持つ運動エネルギーが多くなり燃焼室内ガスの乱れ強さｕも強くなり燃焼が早くなることに対応し、本発の参考例によれば、乱流燃焼速度ＳＴを燃焼室５の総ガス質量ＭＧＡＳ（燃焼室内のガス量）に基づいても算出するので（上記（５８−１）式〜（５８−３）式を参照）、燃焼室５の総ガス質量ＭＧＡＳ（燃焼室内のガス量）の多少に関係なく乱流燃焼速度ＳＴの推定精度を向上できる。

燃焼室内ガスの乱れが少ない状態では層流燃焼速度ＳＬが乱流燃焼速度ＳＴに与える影響が大きいことに対応し、本発明の参考例によれば、層流燃焼速度ＳＴを算出する層流燃焼速度算出処理手順を含み（図１７のステップ３７参照）、乱流燃焼速度ＳＴをこの層流燃焼速度ＳＬに基づいても算出するので（上記（５８−３）式を参照）、燃焼室内ガスの乱れが少ない状態においても乱流燃焼速度ＳＴの推定精度を向上できる。

さて、上記〈２〉では基準クランク角の設定方法を説明したが、その考え方の基本となったアイデアを次に説明する。

燃焼圧が最大となるクランク角θｐｍａｘは、机上で検討した結果によれば、エンジン損失の影響を受けると考えられる。エンジン損失にはエンジンの冷却損失Ｑｃｏｏｌとエンジンの時間損失Ｑｔｉｍｅとが含まれる。このうち燃焼圧が最大となるクランク角θｐｍａｘはエンジンの冷却損失Ｑｃｏｏｌの影響が殆どで、エンジンの時間損失Ｑｔｉｍｅの影響をわずかに受ける。従って、エンジンの冷却損失Ｑｃｏｏｌ、エンジンの時間損失Ｑｔｉｍｅをパラメータとする補正係数ｆ１（Ｑｃｏｏｌ）、ｇ１（Ｑｔｉｍｅ）を導入すれば、次の（７１）式により燃焼圧が最大となるクランク角θｐｍａｘを算出することができる。

θｐｍａｘ＝ＴＨ０・ｆ１（Ｑｃｏｏｌ）・ｇ１（Ｑｔｉｍｅ）
…（７１）
（７１）式のＴＨ０は燃焼圧が最大となるクランク角の基準値で、例えば圧縮上死点［ＴＤＣ］にとっておき、補正係数ｆ１（Ｑｃｏｏｌ）、ｇ１（Ｑｔｉｍｅ）により圧縮上死点より遅角側の値を算出させるようにする。

ここで、エンジンの冷却損失Ｑｃｏｏｌが大きいほど燃焼圧が最大となるクランク角θｐｍａｘは遅角側にずれると考えられるので、補正係数ｆ１（Ｑｃｏｏｌ）は、エンジンの冷却損失Ｑｃｏｏｌが大きくなるほど大きくなるように、また、エンジンの時間損失Ｑｔｉｍｅが大きいほど燃焼圧が最大となるクランク角θｐｍａｘが遅角側にずれると考えられるので、補正係数ｇ１（Ｑｔｉｍｅ）は、エンジンの時間損失Ｑｔｉｍｅが大きくなるほど大きくなるようにそれぞれ設定する。

エンジンの冷却損失Ｑｃｏｏｌを推定するための代表値はエンジン負荷であるので、簡単にはこのエンジン負荷（例えば燃焼室内ガス質量ｍ）をそのままエンジンの冷却損失として用いることができる。また、エンジンの冷却損失Ｑｃｏｏｌを推定するための代表値にはエンジン負荷以外にも燃焼期間があり、燃焼期間が長くなれば（エンジン回転速度Ｎｅが小さくなる）、エンジンの冷却損失Ｑｃｏｏｌは大きくなるため、燃焼期間と逆数の関係にあるエンジン回転速度Ｎｅに応じた回転速度補正率ｆ２（Ｎｅ）を導入し、この回転速度補正率ｆ２（Ｎｅ）で冷却損失補係数ｆ１（Ｑｃｏｏｌ）を補正する。

一方、エンジンの時間損失Ｑｔｉｍｅを推定するための代表値は燃焼速度の逆数であるので、簡単には燃焼速度の逆数をそのままエンジンの時間損失Ｑｔｉｍｅとして用いることができる。また、エンジンの時間損失Ｑｔｉｍｅを推定するための代表値には燃焼速度の逆数以外にも燃焼期間があり、燃焼期間が長くなれば（エンジン回転速度Ｎｅが小さくなる）、時間損失Ｑｔｉｍｅは大きくなるため、燃焼期間と逆数の関係にあるエンジン回転速度Ｎｅに応じた回転速度補正率ｇ２（Ｎｅ）を導入し、この回転速度補正率ｇ２（Ｎｅ）で時間損失補正係数ｇ１（Ｑｔｉｍｅ）を補正する。

従って、このときには上記（７１）式に代えて次の（７２）式により燃焼圧が最大となるクランク角θｐｍａｘを算出することとなる。

θｐｍａｘ＝ＴＨ０・ｆ１（Ｑｃｏｏｌ）・ｆ２（Ｎｅ）
・ｇ１（Ｑｔｉｍｅ）・ｇ２（Ｎｅ） …（７２）
図２０のフローチャートは以上の基本的なアイデアを盛り込んで燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％を算出するためのもので（本発明の第１実施形態）、一定時間毎（例えば１０ｍｓｅｃ毎）に実行する。図１４〜図１９で前述した具体例（本発明の参考例）との関係では、図１９において、図１４のステップ８で算出している本発明の参考例のθｍｂ１０％に代えて、図２０で算出した第１実施形態のθｍｂ１０％を用いることとなる。

フローの説明に入ると、ステップ６１では、図１４のステップ１と同様に、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、バルブリフト量Ｌｉｆｔ［ｍ］、回転速度Ｎｅ［ｒｐｍ］、燃焼室内ガス質量ｍ［ｇ］を読み込む。

ここで、燃焼室内ガス質量ｍとしては、エアフローメータ３２により検出される吸入空気量ＭＡ［ｇ／ｓｅｃ］を用いる。

ステップ６２では、図１４のステップ６と同様に、バルブリフト量Ｌｉｆｔから図４上段を内容とするテーブルを検索することにより、タンブル強度のバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］を算出する。すなわち、ＶＥＬ機構２８の非作動時には第１バルブリフト量補正係数ｋＶＥＬ１を、これに対してＶＥＬ機構２８の作動時には第２バルブリフト量補正係数ｋＶＥＬ２を、バルブリフト量補正係数ｋＶＥＬとして算出する。

ステップ６３では、燃焼室内ガス質量ｍをそのままエンジンの冷却損失Ｑｃｏｏｌ［Ｊ］とし、ステップ６４において、このエンジンの冷却損失Ｑｃｏｏｌから図２１を内容とするテーブルを検索することにより、基本冷却損失補正係数ｆ１（Ｑｃｏｏｌ）［無名数］を算出する。

図２１のように、基本冷却損失補正係数ｆ１（Ｑｃｏｏｌ）はエンジンの冷却損失Ｑｃｏｏｌが大きくなるほど大きくなる値である。すなわち、基本冷却損失補正係数ｆ１（Ｑｃｏｏｌ）は、エンジンの冷却損失Ｑｃｏｌが大きくなるほど燃焼圧が最大となるクランク角θｐｍａｘが圧縮上死点より遅角側に向かうように定めた値である。

ステップ６５では、エンジン回転速度Ｎｅから図２２を内容とするテーブルを検索することにより、冷却損失補正係数の回転速度補正率ｆ２（Ｎｅ）［無名数］を算出し、ステップ６６において、この回転速度補正率ｆ２（Ｎｅ）を上記の基本冷却損失補正係数ｆ１（Ｑｃｏｏｌ）に乗算した値を、冷却損失補正係数ｆ（Ｑｃｏｏｌ）として、つまり次の（７３）式により冷却損失補正係数ｆ（Ｑｃｏｏｌ）［無名数］を算出する。

ｆ（Ｑｃｏｏｌ）＝ｆ１（Ｑｃｏｏｌ）・ｆ２（Ｎｅ） …（７３）
図２２のように、回転速度補正率ｆ２（Ｎｅ）は回転速度Ｎｅが大きくなるほど小さくなる値である。すなわち、回転速度補正率ｆ２（Ｎｅ）は、燃焼期間が長くなるほど（エンジン回転速度が小さくなるほど）燃焼圧が最大となるクランク角θｐｍａｘがより遅角側に向かうように定めた値である。

ステップ６７では、バルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣ、燃焼室内ガス質量ｍ、エンジン回転速度Ｎｅを用いて次の（７４）式により燃焼速度［ｍ／ｓｅｃ］を算出する。

燃焼速度＝ｋＶＥＬ・ＩＶＣ・ｍ＾ｃ４・Ｎｅ＾ｃ５ …（７４）
ただし、ｃ４、ｃ５：適合定数、
（７４）式を上記（１６）式と比較すれば分かるように、（７４）式は正確には乱れ強さｕを求める式であるが、ここでは、乱れ強さｕを簡易的にそのまま燃焼速度として設定するようにしたものである。

ステップ６８では、この燃焼速度の逆数をそのままエンジンの時間損失Ｑｔｉｍｅ［Ｊ］とし、ステップ６９において、このエンジンの時間損失Ｑｔｉｍｅから図２３を内容とするテーブルを検索することにより、基本時間損失補正係数ｇ１（Ｑｔｉｍｅ）［無名数］を算出する。

図２３のように、基本時間損失補正係数ｇ１（Ｑｔｉｍｅ）はエンジンの時間損失Ｑｔｉｍｅが大きくなるほど大きくなる値である。すなわち、基本時間損失補正係数ｇ１（Ｑｔｉｍｅ）は、エンジンの時間損失Ｑｔｉｍｅが大きくなるほど燃焼圧が最大となるクランク角θｐｍａｘが遅角側に向かうように定めた値である。

ステップ７０では、エンジン回転速度Ｎｅから図２４を内容とするテーブルを検索することにより、時間損失補正係数の回転速度補正率ｇ２（Ｎｅ）［無名数］を算出し、ステップ７１において、この回転速度補正率ｇ２（Ｎｅ）を上記の基本時間損失補正係数ｇ１（Ｑｔｉｍｅ）に乗算した値を、時間損失補正係数ｇ（Ｑｔｉｍｅ）として、つまり次の（７５）式により時間損失補正係数ｇ（Ｑｔｉｍｅ）［無名数］を算出する。

ｇ（Ｑｔｉｍｅ）＝ｇ１（Ｑｔｉｍｅ）・ｇ２（Ｎｅ） …（７５）
図２４のように、回転速度補正率ｇ２（Ｎｅ）は回転速度Ｎｅが大きくなるほど小さくなる値である。すなわち、回転速度補正率ｆ２（Ｎｅ）は、燃焼期間が長くなるほど（エンジン回転速度が小さくなるほど）燃焼圧が最大となるクランク角θｐｍａｘがより遅角側に向かうように定めた値である。

ステップ７２では、この時間損失補正係数ｇ（Ｑｔｉｍｅ）と、上記の冷却損失補正係数ｆ（Ｑｃｏｏｌ）とで、圧縮上死点である燃焼圧が最大となるクランク角の基準値ＴＨ０［ＴＤＣ］を補正して、つまり次の（７６）式により燃焼圧が最大となるクランク角θｐｍａｘ［ｄｅｇＡＴＤＣ］を算出する。

θｐｍａｘ＝ＴＨ０・ｆ（Ｑｃｏｏｌ）・ｇ（Ｑｔｉｍｅ） …（７６）
ここで、基準とする運転状態でのエンジンの冷却損失、エンジンの時間損失、エンジン回転速度、基本冷却損失補正係数、基本時間損失補係数をそれぞれＱｃｏｏｌ０、Ｑｔｉｍｅ０、Ｎｅ０、ｆ１（Ｑｃｏｏｌ０）、ｇ１（Ｑｔｉｍｅ０）とすると、適合により基準とする運転状態での基本冷却損失補正係数ｆ１（Ｑｃｏｏｌ０）、基準とする運転状態での基本時間損失補係数ｇ１（Ｑｔｉｍｅ１）を定めることができる。この基準とする運転状態でのθｐｍａｘをθｐｍａｘ０とすれば、（７６）式より
θｐｍａｘ０＝ＴＨ０・ｆ（Ｑｃｏｏｌ０）・ｇ（Ｑｔｉｍｅ０）
…（７７）
であり、θｐｍａｘ０は圧縮上死点後の所定のクランク角位置に定まる。

いま、エンジンの時間損失は一定であるとして、基準とする運転状態でのエンジンの冷却損失Ｔｃｏｏｌ０よりエンジンの冷却損失Ｑｃｏｏｌが大きくなると、燃焼圧が最大となるクランク角θｐｍａｘはθｐｍａｘ０より遅角側に移動し、この逆に基準とする運転状態でのエンジンの冷却損失Ｔｃｏｏｌ０よりエンジンの冷却損失Ｑｃｏｏｌが小さくなると燃焼圧が最大となるクランク角θｐｍａｘはθｐｍａｘ０より進角側に移動する。上記の図２１はこうしたエンジンの冷却損失Ｑｃｏｏｌの、Ｑｃｏｏｌ０からの相違をθｐｍａｘに反映させるものである。

同様にして、エンジンの冷却損失Ｑｃｏｏｌは一定であるとして、基準とする運転状態でのエンジンの時間損失Ｑｔｉｍｅ０よりエンジンの時間損失Ｑｔｉｍｅが大きくなると、燃焼圧が最大となるクランク角θｐｍａｘはθｐｍａｘ０より遅角側に移動し、この逆に基準とする運転状態でのエンジンの時間損失Ｑｔｉｍｅ０よりエンジンの時間損失Ｑｔｉｍｅが小さくなると、燃焼圧が最大となるクランク角θｐｍａｘはθｐｍａｘ０より進角側に移動する。上記の図２３はこうしたエンジンの時間損失Ｑｔｉｍｅの、Ｑｔｉｍｅ０からの相違をθｐｍａｘに反映させるものである。

一方、回転速度補正率ｆ２（Ｎｅ）、ｇ２（Ｎｅ）は燃焼期間の相違を燃焼圧が最大となるクランク角θｐｍａｘに反映させるためのものである。具体的には基準とする運転状態での回転速度Ｎｅ０のときに回転速度補正率ｆ２（Ｎｅ）、ｇ２（Ｎｅ）はいずれも１．０である。基準とする運転状態での回転速度Ｎｅ０より回転速度Ｎｅが大きくなると燃焼圧が最大となるクランク角θｐｍａｘはθｐｍａｘ０より進角側に移動し、この逆に基準とする運転状態での回転速度Ｎｅ０より回転速度Ｎｅが小さくなると燃焼圧が最大となるクランク角θｐｍａｘはθｐｍａｘ０より遅角側に移動する。そこで、図２２、図２４によりこうしたエンジン回転速度Ｎｅの、Ｎｅ０からの相違（エンジン回転速度と逆数の関係にある燃焼期間の、Ｎｅ０での燃焼期間からの相違）を冷却損失補正係数ｆ（Ｑｔｉｍｅ）、時間損失補正係数ｇ（Ｑｔｉｍｅ）にそれぞれ反映させる。

ステップ７３では、このようにして算出した燃焼圧が最大となるクランク角θｐｍａｘを用いて次の（７８）式により燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％［ｄｅｇＡＴＤＣ］を算出する。

θｍｂ１０％＝ａ・θｐｍａｘ−ｂ …（７８）
ただし、ａ、ｂ：適合値、
そして、本発明の参考例の図１９において本発明の参考例の図１４のステップ８で算出しているθｍｂ１０％代えて、（７８）式で算出したθｍｂ１０％を用いて基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出する。すなわち、図１４のステップ８を除く残りの部分、図１６、図１７、図１８、図１９は第１実施形態でも流用する。

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

燃焼圧が最大となるクランク角θｐｍａｘは運転条件によって変動するが、その値はエンジンの損失によって決まることに鑑み、第１実施形態（請求項１、１１に記載の発明）によれば、燃焼圧が最大となるクランク角θｐｍａｘを、エンジンの損失に基づいて算出し（図２０のステップ７２参照）、この燃焼圧が最大となるクランク角θｐｍａｘに基づいてＭＢＴの得られる基本点火時期を算出するようにしたので（図２０のステップ７３、図１９のステップ５１〜５３参照）、燃焼圧が最大となるクランク角θｐｍａｘが運転条件により変動しても、その変動するθｐｍａｘを精度よく推定できることから、ＭＢＴの得られる基本点火時期の予測精度を向上させることができる。

特に、ＶＥＬ機構２８やＶＴＣ機構２７（可変動弁装置）を備えるエンジンでは、ＶＥＬ機構２８及びＶＴＣ機構２７を備えていないエンジンよりも燃焼圧が最大となるクランク角θｐｍａｘが大きく変動するので、こうしたエンジンに特に有効である。ただし、ＶＥＬ機構２８及びＶＴＣ機構２７（可変動弁装置）を備えていないエンジンに対しても適用があることはもちろんである。

第１実施形態（請求項２、１２に記載の発明）によれば、燃焼圧が最大となるクランク角θｐｍａｘに基づいて、燃焼室に供給された燃料に対する燃料質量の比率を表す燃焼質量割合が１０％（所定の基準割合）に達するクランク角位置であるθｍｂ１０％を算出し（図２０のステップ７３参照）、そのθｍｂ１０％（基準割合に達するクランク角位置）に基づいて基本点火時期ＭＢＴＣＡＬ（ＭＢＴの得られる基本点火時期）を算出している（図１９のステップ５１、５２、５３参照）。つまり、所定の基準割合に達するクランク角位置とは、燃焼圧が最大となるクランク角θｐｍａｘと、燃焼質量割合がある値に達するクランク角位置とで一番相関のあるクランク角位置（具体的には燃焼質量割合が１０％に達するときのクランク角位置であるθｍｂ１０％）であり、この一番相関のあるクランク角位置に基づいて基本点火時期ＭＢＴＣＡＬ（ＭＢＴの得られる基本点火時期）を算出することで、基本点火時期ＭＢＴＣＡＬの算出精度を向上できる。

燃焼圧が最大となるクランク角θｐｍａｘが変動するのは、エンジンの冷却損失Ｑｃｏｏｌによる影響がほとんどであることに鑑み、第１実施形態（請求項１、１１に記載の発明）によれば、エンジンの冷却損失Ｑｃｏｏｌに基づいて燃焼圧が最大となるクランク角θｐｍａｘを算出するので（図２０のステップ６４、６６、７２参照）、燃焼圧が最大となるクランク角θｐｍａｘの算出を容易に行うことができる。

燃焼圧が最大となるクランク角θｐｍａｘの変動には、エンジンの冷却損失Ｑｃｏｏｌの他に、エンジンの時間損失Ｑｔｉｍｅによる影響もあることに鑑み、第１実施形態（請求項１、１１に記載の発明）によれば、エンジンの冷却損失Ｑｃｏｏｌに加えてエンジンの時間損失Ｑｔｉｍｅに基づいても燃焼圧が最大となるクランク角θｐｍａｘを算出するので（図２０のステップ６９、７１、７２参照）、燃焼圧が最大となるクランク角θｐｍａｘの推定精度をさらに向上できる。

次に、図２５のフローチャートは第２実施形態の燃焼質量割合が１０％に達するときのクランク角位置であるθｍｂ１０％を算出するためのもので、第１施形態の図２０と置き換わるものである。図２０と同一部分には同一のステップ番号を付けている。

第１実施形態では、燃焼室内ガス質量ｍをそのままエンジンの冷却損失Ｑｃｏｏｌとし、燃焼速度の逆数をそのままエンジンの時間損失Ｑｔｉｍｅとして簡素に構成したが（図２０のステップ６３、６８参照）、第２実施形態では、エンジンの冷却損失Ｑｃｏｏｌをウォッシーニの式により精密に、かつ燃焼速度についても精密に算出するようにしたものである。

第１実施形態との違いを主に説明すると、ステップ８２では、エンジンの冷却損失Ｑｃｏｏｌ［Ｊ］を算出する。このエンジンの冷却損失Ｑｃｏｏｌの算出については図２６のフローにより説明する。

図２６（図２５のステップ８２のサブルーチン）において、ステップ９１では、図１４のステップ１と同様に、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、図２５のステップ６２で算出されているバルブリフト量補正係数ｋＶＥＬ［１／ｄｅｇ］、温度センサ４３により検出されるコレクタ内温度ＴＣＯＬ［Ｋ］、圧力センサ４４により検出されるコレクタ内圧力ＰＣＯＬ［Ｐａ］、温度センサ４５により検出される排気温度ＴＥＸＨ［Ｋ］、内部不活性ガス率ＭＲＥＳＲ［％］、クランク角センサにより検出されるエンジン回転速度Ｎｅ［ｒｐｍ］、エアフローメータ３２により検出される吸入空気量ＭＡ［ｇ／ｓｅｃ］を読み込む。

ステップ９２、９３、９４では、図１４のステップ２、３、４と同様に、燃焼室５の吸気弁閉時期ＩＶＣにおける容積Ｖｉｖｃ［ｍ³］、燃焼室５の吸気弁閉時期ＩＶＣにおける温度Ｔｉｖｃ［Ｋ］、燃焼室５の吸気弁閉時期ＩＶＣにおける圧力Ｐｉｖｃ［ｋＰａ］を算出する。

ステップ９５〜１０１はウォッシーニの式に必要となる値（ガス流速ＶＳ、係数ｃ１、燃焼室内圧力Ｐ０、モータリング時圧力ＰＭ、Ｗ、未燃ガス温度ＴＵＧＡＳ０、ＨＣ）を算出する部分である。

まず、ステップ９５では、バルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣ、回転速度Ｎｅを用いて、次の（８１）式により燃焼室内ガス流速ＶＳ［ｍ／ｓｅｃ］を算出する。

ＶＳ＝ｋＶＥＬ・ＩＶＣ・Ｎｅ …（８１）
ここでは、バルブリフト量補正係数ｋＶＥＬを本発明の参考例から流用しているが、ガス流速ＶＳの単位は［ｍ／ｓｅｃ］であるので、バルブリフト量補正係数ｋＶＥＬの単位は［ｍ／ｄｅｇ］でなければならない。この点を明確にするには、バルブリフト量Ｌｉｆｔから図４上段を内容とするテーブルと同様のテーブルを検索することにより、第２バルブリフト量補正係数ｋ２ＶＥＬ［ｍ／ｄｅｇ］を算出し、（８１）式のバルブリフト量補正係数ｋＶＥＬに代えて、この第２バルブリフト量補正係数ｋ２ＶＥＬを用いればよい。

本発明の参考例では、バルブリフト量補正係数ｋＶＥＬはタンブルのバルブリフト量補正係数であったが、第２実施形態では、燃焼室内ガス流動を必ずしもタンブルに限定しているわけでないので、燃焼室内ガス流動がスワールの場合にも適用がある。

また、（８１）式によれば、バルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣが入っているので、第２実施形態においてもＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象としているかのようにみえるが、そうではなく、第２実施形態はＶＥＬ機構２８及びＶＴＣ機構２７を備えないエンジンをも対象としている。ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象とするときには、バルブリフト量補正係数ｋＶＥＬが切換えられる２値となり、かつ吸気弁閉時期ＶＴＣが可変値となるが、ＶＥＬ機構２８及びＶＴＣ機構２７を備えないエンジンを対象とするときにはバルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣとも固定値になるだけである。

ステップ９６では、この燃焼室内ガス流速ＶＳとエンジン回転速度Ｎｅを用いて次の（８２）式により係数ｃ１を算出する。

ｃ１＝２．２８＋０．３０８・ＶＳ／２／ＳＴ／Ｎｅ …（８２）
ただし、ＳＴ：ピストンの全ストローク［ｍ］、
ステップ９７、９８では、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ、吸入空気量ＭＡを用いて次の（８３）式、（８４）式により燃焼室内圧力Ｐ０、燃焼室５のモータリング時における圧力ＰＭを算出する。

Ｐ０＝Ｐｉｖｃ・（Ｖ／Ｖｉｖｃ）＾κ＋ＭＡ・ＭＢＡＶＥ・ＤＰＣ
…（８３）
ＰＭ＝Ｐｉｖｃ・（Ｖｉｖｃ／Ｖ）＾κ …（８４）
ただし、Ｖ ：燃焼室内容積［ｍ³］、
κ ：比熱比（固定値で１．３〜１．４）、
ＭＢＡＶＥ：基準燃焼質量割合［％］、
ＤＰＣ ：吸入空気１ｇ当たりの燃焼による圧力上昇分
［ｋＰａ／ｇ］、
（８３）式右辺第１項は容積変化による圧力変化分（つまり燃焼室５のモータリング時における圧力）、第２項は燃焼質量割合が１０％に達するクランク角位置（θｍｂ１０％）までの燃焼による圧力上昇分である。従って、（８３）式の基準燃焼質量割合ＭＢＡＶＥとは１０％である。

上記（８３）式、（８４）式の燃焼室内容積Ｖは適合値とする。

ステップ９９では、係数ｃ１、回転速度Ｎｅ、燃焼室５の吸気弁閉時期における温度Ｔｉｖｃ、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室５の吸気弁閉時期における容積Ｖｉｖｃ、燃焼室内圧力Ｐ０、燃焼室５のモータリング時における圧力ＰＭを用いて次の（８５）式によりＷを算出する。

Ｗ＝ｃ１・２・ＳＴ・Ｎｅ
＋ｃ２・Ｖｃ・Ｔｉｖｃ／Ｐｉｖｃ／Ｖｉｖｃ・（Ｐ０−ＰＭ）
…（８５）
ただし、ＳＴ ：ピストンの全ストローク［ｍ］、
Ｖｃ ：隙間容積［ｍ³］、
ｃ１、ｃ２：係数、
ステップ１００では、燃焼室５の吸気弁閉時期における温度Ｔｉｖｃ、燃焼室５の吸気弁閉時期における圧力Ｐｉｖｃ、燃焼室内圧力Ｐ０を用いて次の（８６）式により未燃ガス温度ＴＵＧＡＳ０［Ｋ］を算出する。

ＴＵＧＡＳ０＝Ｔｉｖｃ・（Ｐ０／Ｐｉｖｃ）＾（（κ−１）／κ）
…（８６）
ただし、κ：比熱比（固定値で１．３〜１．４）、
ステップ１０１では、上記のＷ、燃焼室内圧力Ｐ０、未燃ガス温度ＴＵＧＡＳ０を用いて次の（８７）式によりＨＣを算出する。

ＨＣ＝３．２６・Ｄ＾−０．２・Ｐ０＾０．８・ＴＵＧＡＳ０＾−０．５５
・Ｗ＾０．８
…（８７）
ただし、Ｄ：シリンダボア径［ｍ］、
ステップ１０２では、このＨＣ、未燃ガス温度ＴＵＧＡＳ０、回転速度Ｎｅを用いて冷却損失の式（ウォッシーニの式）である次の（８８）式によりエンジンの冷却損失Ｑｃｏｏｌ［Ｊ］を算出する。

Ｑｃｏｏｌ＝ＡＣＹＬ・ＨＣ・（ＴＵＧＡＳ０−Ｔｗ）・ｄθ・６／Ｎｅ
…（８８）
ただし、ＡＣＹＬ：燃焼室内壁面積［ｍ²］、
Ｔｗ ：燃焼室壁温［Ｋ］、
ｄθ ：平均燃焼期間［ｄｅｇ］、
燃焼室内壁面積ＡＣＹＬ、燃焼室壁温Ｔｗ、平均燃焼期間ｄθは適合値である。

このようにしてエンジンの冷却損失Ｑｃｏｏｌを算出したら、図２５に戻り、ステップ６４〜６６の操作を第１実施形態の図２０のステップ６４〜６６と同じに実行する。

図２５のステップ８３では、エンジンの時間損失Ｑｔｉｍｅ［Ｊ］を算出する。このエンジンの時間損失Ｑｔｉｍｅの算出については図２７のフローにより説明する。

図２７（図２５のステップ８３のサブルーチン）において、ステップ１１１では、吸気弁閉時期ＩＶＣ［ｄｅｇＡＴＤＣ］、内部不活性ガス率ＭＲＥＳＲ［％］、目標当量比Ｔｆｂｙａ、クランク角センサにより検出されるエンジン回転速度Ｎｅ［ｒｐｍ］、エアフローメータ３２により検出される吸入空気量ＭＡ［ｇ／ｓｅｃ］、図２５のステップ６２で算出しているバルブリフト量補正係数ｋＶＥＬ、図２６のステップ９７で算出している燃焼室内圧力Ｐ０［ｋＰａ］、図２６のステップ１００で算出している未燃ガス温度ＴＵＧＡＳ０［Ｋ］を読み込む。

ステップ１１２では、目標当量比Ｔｆｂｙａから次の（８９）式〜（９１）式により標準状態での層流燃焼速度ＳＬ０［ｍ／ｓｅｃ］、係数ａ、ｂ［無名数］を求め、ステップ１１３でこれら標準状態での層流燃焼速度ＳＬ０、係数ａ、ｂと、未燃ガス温度ＴＵＧＡＳ０、燃焼室内圧力Ｐ０、内部不活性ガス率ＭＲＥＳＲを用いて次の（９２）式により層流燃焼速度ＳＬ［ｍ／ｓｅｃ］を算出する。

ＳＬ０＝０．３０５−０．５４９（Ｔｆｂｙａ−１．２１）＾２
…（８９）
ａ＝２．４−０．２７１・Ｔｆｂｙａ＾３．５１ …（９０）
ｂ＝−０．３５７＋０．１４・Ｔｆｂｙａ＾２．７７ …（９１）
ＳＬ＝ＳＬ０・（ＴＵＧＡＳ０／２９８）＾ａ・（Ｐ０／１０１）＾ｂ
・（１−２．０６・ＭＲＥＳＲ＾０．７７）
…（９２）
（８９）式〜（９２）式は、上記（５４）式〜（５７）式と同様の式である。

ステップ１１４では、この層流燃焼速度ＳＬ、バルブリフト量補正係数ｋＶＥＬ、吸気弁閉時期ＩＶＣ、吸入空気量ＭＡ、エンジン回転速度Ｎｅを用いて次の（９３）式により燃焼速度Ｓｂ［ｍ／ｓｅｃ］を算出する。

Ｓｂ＝ＳＬ＋（ｋＶＥＬ・ＩＶＣ・ＭＡ＾ｅ・Ｎｅ＾ｆ／ＳＬ）＾ｇ
…（９３）
ただし、ｆ、ｇ：適合値、
（９３）式は上記（５８）式と同様の式である。

このように、第２実施形態でも、乱流燃焼速度の算出式（（９３）式右辺第２項）に新たにバルブリフト量補正係数ｋＶＥＬと吸気弁閉時期ＩＶＣと吸入空気量ＭＡとを新たに導入している。ただし、前述したように、ＶＥＬ機構２８及びＶＴＣ機構２７を備えるエンジンを対象とするだけでなく、ＶＥＬ機構２８及びＶＴＣ機構２７を備えないエンジンをも対象としている。

ステップ１１５ではこの燃焼速度Ｓｂ、回転速度Ｎｅを用いて、次の（９４）式によりエンジンの時間損失Ｑｔｉｍｅ［Ｊ］を算出する。

Ｑｔｉｍｅ＝ｄ／（Ｓｂ・Ｎｅ） …（９４）
ただし、ｄ：適合値、
（９４）式は、エンジンの時間損失Ｑｔｉｍｅが燃焼速度Ｓｂ及び回転速度Ｎｅに反比例するとする式である。

このようにしてエンジンの時間損失Ｑｔｉｍｅを算出したら、図２５に戻り、ステップ６９〜７３の操作を第１実施形態の図２０のステップ６９〜７３と同じに実行して燃焼質量割合が１０％に達するときのクランク角位置であるθｍｂ１０％を算出する。

そして、本発明の参考例の図１９において本発明の参考例の図１４のステップ８で算出しているθｍｂ１０％代えて、図２５のステップ７３で算出したθｍｂ１０％を用いて基本点火時期ＭＢＴＣＡＬ［ｄｅｇＢＴＤＣ］を算出する。すなわち、図１４のステップ８を除く残りの部分、図１６、図１７、図１８、図１９は第２実施形態でも流用する。

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

エンジンの冷却損失Ｑｃｏｏｌは、未燃ガス温度ＴＵＧＡＳ０（燃焼室内の温度）、燃焼室内圧力Ｐ０（燃焼室内の圧力）及び燃焼室内ガス流速ＶＳ（燃焼室内のガス流速）と燃焼期間（回転速度Ｎｅより換算）による影響が大きい。第２実施形態（請求項３、１３に記載の発明）によれば、それらに基づいて、ウォッシーニの式によりエンジンの冷却損失Ｑｃｏｏｌを算出するので（図２６のステップ９５〜１０２参照）、エンジンの冷却損失Ｑｃｏｏｌを正確に算出できる。

燃焼室内ガス流速ＶＳの算出を正確に行うのは非常に困難であるが、燃焼室内ガス流速ＶＳはほぼエンジン回転速度Ｎｅに比例することに鑑み、第２実施形態（請求項４、１４に記載の発明）によれば、少なくともエンジン回転速度Ｎｅに基づいて燃焼室内ガス流速ＶＳを算出するので（図２６のステップ９５参照）、燃焼室内ガス流速ＶＳを簡素に推定できる。

燃焼室内ガス流速ＶＳは、ＶＥＬ機構２８が作動して、ＶＥＬ機構２８の非作動時よりバルブリフト量が小さくその吸気弁１５の開口面積が小さいと、その開口面積で吸気弁１５を通過して燃焼室５に流れ込む吸気の流速が制限を受けてしまうため、ＶＥＬ機構２８の非作動状態にあって吸気弁１５の開口面積が大きい場合とは違った燃焼室内ガス流速ＶＳとなる。つまり、ＶＥＬ機構２８（可変動弁装置）の作動、非作動により運転条件に応じてバルブリフト量（吸気弁１５のリフト特性）が違ったものとなり、その運転条件に応じたバルブリフト量（吸気弁１５のリフト特性）の相違により燃焼室内ガス流速ＶＳが、ＶＥＬ機構２８（可変動弁装置）を備えないエンジンの場合よりも大きく変化する。第２実施形態（請求項５、１５に記載の発明）によれば、ＶＥＬ機構２８（可変動弁装置）を備える場合に、燃焼室内ガス流速ＶＳをエンジン回転速度Ｎｅに加えて、バルブリフト量（吸気弁のリフト特性）に基づいても算出するので（図２６のステップ９５参照）、運転条件に応じたバルブリフト量（吸気弁のリフト特性）の相違で吸気弁１５の開口面積が変化する場合においても、燃焼室内ガス流速ＶＳを精度良く算出できる。

エンジンの時間損失Ｑｔｉｍｅそのものをエンジンコントローラ３１で計算することは困難であるが、第２実施形態（請求項６、１６に記載の発明）によれば、燃焼速度が小さいほどエンジンの時間損失Ｑｔｉｍｅが大きくなるという関係を利用して、エンジンの時間損失Ｑｔｉｍｅを、燃焼室内ガスの燃焼速度Ｓｂに基づいて算出するので（図２７のステップ１１５参照）、エンジンの時間損失Ｑｔｉｍｅの算出が容易になる。

燃焼が進んだ期間では、燃焼室内ガスの燃焼速度Ｓｂはエンジン回転速度Ｎｅでほぼ決まることに鑑み、第２実施形態（請求項７、１７に記載の発明）によれば、燃焼室内ガスの燃焼速度Ｓｂを、少なくともエンジン回転速度Ｎｅに基づいて算出するので（図２７のステップ１１４参照）、比較的推定精度を落とすことなく、容易にエンジンの時間損失Ｑｔｉｍｅを算出できる。

ＶＥＬ機構２８が作動してバルブリフト量Ｌｉｆｔが大きい状態から小さい状態へと変化しあるいはこの逆にＶＥＬ機構２８の非作動でバルブリフト量Ｌｉｆｔが小さい状態から大きい状態へと変化したり、ＶＴＣ機構２７により吸気弁１５の開閉タイミングが遅角側や進角側へと変化すると（可変動弁装置が働いて吸気弁１５のリフト特性が変化すると）、燃焼室内ガス流動が変化するため、燃焼室内ガスの燃焼速度Ｓｂが、ＶＥＬ機構２８、ＶＴＣ機構２７を備えないエンジンの場合よりも大きく変化する。この場合に、第２実施形態（請求項８、１８に記載の発明）によれば、エンジン回転速度Ｎｅに加えて、バルブリフト量（吸気弁１５のリフト特性）に基づいても燃焼室内ガスの燃焼速度Ｓｂを算出するので（図２７のステップ１１４のｋＶＥＬ、ＩＶＣ参照）、吸気弁１５のバルブリフト量や吸気弁１５の開閉タイミング（吸気弁のリフト特性）が相違する場合であっても、燃焼速度Ｓｂを精度良く算出できる。

燃焼初期では燃焼室内ガスの燃焼速度Ｓｂは、燃焼室内ガスの温度、圧力、不活性ガス率の影響が大きいことに鑑み、第２実施形態（請求項９、１９に記載の発明）によれば、エンジン回転速度Ｎｅ、吸気弁１５のバルブリフト量や開閉タイミング（吸気弁のリフト特性）に加えて、未燃ガス温度ＴＵＧＡＳ０（燃焼室内ガスの温度）、燃焼室内圧力Ｐ０（燃焼室内の圧力）、不活性ガス率ＭＲＥＳＲに基づいても層流燃焼速度ＳＬ（燃焼室内ガスの燃焼速度Ｓｂ）を算出するので（図２７のステップ１１３参照）、特に燃焼初期における層流燃焼速度ＳＬ（燃焼室内ガスの燃焼速度Ｓｂ）の推定精度が向上する。

第２実施形態に対して、さらに、スワールコントロールバルブ、タンブルコントロールバルブのように、燃焼室５への吸気流動を変化させ得る可変吸気流動装置を備えさせることが考えられる。例えばスワールコントロールバルブを作動させた状態では、スワールコントロールバルブを作動させていない状態より、燃焼室内ガスの乱流燃焼速度ＳＴが大きくなるので、スワールコントロールバルブを作動させるときには、バルブの作動に伴って乱流燃焼速度ＳＴが大きくなる分だけ、乱流燃焼速度ＳＴが大きくなる側に補正してやればよい。すなわち、燃焼室内ガスの燃焼速度Ｓｂをスワールコントロールバルブの作動状態（作動状態、非作動状態）（可変吸気流動装置の状態）に応じても算出させる（請求項１０、２０に記載の発明）。これにより、可変吸気流動装置を備えている場合においても、燃焼室内ガスの燃焼速度Ｓｂを精度よく算出できる。

実施形態では、ＶＥＬ機構２８及びＶＴＣ機構２７を備える場合で説明したが、本発明は、ＶＥＬ機構２８、ＶＴＣ機構２７のいずれかを備える場合あるいはいずれも備えない場合にも適用がある。

請求項１に記載の発明において、燃焼圧最大クランク角算出処理手順は図２０のステップ７２により、基本点火時期算出処理手順は図１９のステップ５１、５２、５３によりそれぞれ果たされている。

請求項１３に記載の発明において、燃焼圧最大クランク角算出手段の機能は図２０のステップ７２により、基本点火時期算出手段の機能は図１９のステップ５１、５２、５３によりそれぞれ果たされている。

本発明の参考例のエンジンの点火時期制御装置の概略構成図。 燃焼室内におけるタンブルの生成から消滅までをピストンの動きに合わせて表す図。 吸気弁閉時期に対するタンブル強度の特性図。 バルブリフト量に対するバルブリフト量補正係数の特性図。 燃焼圧が最大となるクランク角θｐｍａｘとθｍｂ０％、θｍｂ１０％、θｍｂ６０％との関係を表す特性図。 燃焼圧が最大となるクランク角θｐｍａｘとエンジン回転速度との関係を表す特性図。 燃焼質量割合と、燃焼質量割合が０％、１０％、２０％、３０％、４０％、５０％、６０％となるときのクランク角位置であるθｍｂ０％、θｍｂ１０％、θｍｂ２０％、θｍｂ３０％、θｍｂ４０％、θｍｂ５０％、θｍｂ６０％をそれぞれθｐｍａｘで直線近似した場合のＲの２乗値との関係を表す特性図。 コンベンショナルエンジンに適用している先行点火時期制御方法における基準クランク角と、本発明の参考例の変更後点火時期制御方法における基準クランク角の違いを説明するための燃焼質量割合の特性図。 吸入空気量、エンジン回転速度に対する、燃焼圧が最大となるクランク角θｐｍａｘの特性図。 バルブリフト量に対する流速係数の特性図。 本発明の参考例の基本点火時期の算出方法を説明するための特性図。 燃焼質量割合がｘ％となるときのクランク角位置であるθｍｂ ｘ％と、燃焼室５のθｍｂ ｘ％時における未燃ガス容積Ｖｍｂ ｘ％との関係を示す特性図。 燃料噴射パルス幅と、燃焼によるトータルの圧力上昇分との関係を示す特性図。 物理量の算出を説明するためのフローチャート。 エンジンのクランクシャフトとコネクティングロッドの位置関係を説明するダイアグラム。 内部不活性ガス率の算出を説明するためのフローチャート。 燃焼期間の算出を説明するためのフローチャート。 着火遅れ時間の算出を説明するためのフローチャート。 基本点火時期の算出を説明するためのフローチャート。 本発明の第１実施形態の燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％の算出を説明するためのフローチャート。 基本冷却損失補正係数の特性図。 冷却損失補正係数の回転速度補正率の特性図。 基本時間損失補正係数の特性図。 時間損失補正係数の回転速度補正率の特性図。 第２実施形態の燃焼質量割合が１０％となるときのクランク角位置であるθｍｂ１０％の算出を説明するためのフローチャート。 第２実施形態の冷却損失の算出を説明するためのフローチャート。 第２実施形態の時間損失の算出を説明するためのフローチャート。

符号の説明

１ エンジン
５ 燃焼室
１１ 点火装置
１５ 吸気弁
２１ 燃料インジェクタ
２７ ＶＴＣ機構（可変動弁装置）
２８ ＶＥＬ機構（可変動弁装置）
３１ エンジンコントローラ
３３、３４ クランク角センサ

Claims (20)

1. 燃焼圧が最大となるクランク角を、エンジンの冷却損失、またはエンジンの時間損失に基づいて算出する燃焼圧最大クランク角算出処理手順と、
   この燃焼圧が最大となるクランク角に基づいてＭＢＴの得られる基本点火時期を算出する基本点火時期算出処理手順と
   を含むことを特徴とするエンジンの点火時期制御方法。
2. 前記燃焼圧が最大となるクランク角に基づいて、燃焼室に供給された燃料に対する燃料質量の比率を表す燃焼質量割合が所定の基準割合に達するクランク角位置を算出し、その基準割合に達するクランク角位置に基づいて前記ＭＢＴの得られる基本点火時期を算出することを特徴とする請求項１に記載のエンジンの点火時期制御方法。
3. 前記エンジンの冷却損失を、燃焼室内温度と、燃焼室内圧力と、燃焼室内ガス流速と、エンジン回転速度とに基づいて算出することを特徴とする請求項１に記載のエンジンの点火時期制御方法。
4. 前記燃焼室内ガス流速を少なくともエンジン回転速度に基づいて算出することを特徴とする請求項３に記載のエンジンの点火時期制御方法。
5. 吸気弁のリフト特性を変化させ得る可変動弁装置を備え、
   前記燃焼室内ガス流速を少なくともエンジン回転速度と、前記吸気弁のリフト特性とに基づいて算出することを特徴とする請求項３に記載のエンジンの点火時期制御方法。
6. 前記エンジンの時間損失を、燃焼室内ガスの燃焼速度に基づいて算出することを特徴とする請求項１に記載のエンジンの点火時期制御方法。
7. 前記燃焼室内ガスの燃焼速度を、少なくともエンジン回転速度に基づいて算出することを特徴とする請求項６に記載のエンジンの点火時期制御方法。
8. 吸気弁のリフト特性を変化させ得る可変動弁装置を備え、
   前記燃焼室内ガスの燃焼速度を、少なくともエンジン回転速度と前記吸気弁のリフト特性とに基づいて算出することを特徴とする請求項６に記載のエンジンの点火時期制御方法。
9. 吸気弁のリフト特性を変化させ得る可変動弁装置を備え、
   前記燃焼室内ガスの燃焼速度を、燃焼室内ガスの温度、燃焼室ガスの圧力、不活性ガス率、エンジン回転速度及び前記吸気弁のリフト特性に基づいて算出することを特徴とする請求項６に記載のエンジンの点火時期制御方法。
10. 燃焼室への吸気流動を変化させ得る可変吸気流動装置を備え、
    前記燃焼室内ガスの燃焼速度をこの可変吸気流動装置の状態に応じても算出することを特徴とする請求項８または９に記載のエンジンの点火時期制御方法。
11. 燃焼圧が最大となるクランク角を、エンジンの冷却損失、またはエンジンの時間損失に基づいて算出する燃焼圧最大クランク角算出手段と、
    この燃焼圧が最大となるクランク角に基づいてＭＢＴの得られる基本点火時期を算出する基本点火時期算出手段と
    を含むことを特徴とするエンジンの点火時期制御装置。
12. 前記燃焼圧が最大となるクランク角に基づいて、燃焼室に供給された燃料に対する燃料質量の比率を表す燃焼質量割合が所定の基準割合に達するクランク角位置を算出し、その基準割合に達するクランク角位置に基づいて前記ＭＢＴの得られる基本点火時期を算出することを特徴とする請求項１１に記載のエンジンの点火時期制御装置。
13. 前記エンジンの冷却損失を、燃焼室内温度と、燃焼室内圧力と、燃焼室内ガス流速と、エンジン回転速度とに基づいて算出することを特徴とする請求項１１に記載のエンジンの点火時期制御装置。
14. 前記燃焼室内ガス流速を少なくともエンジン回転速度に基づいて算出することを特徴とする請求項１３に記載のエンジンの点火時期制御装置。
15. 吸気弁のリフト特性を変化させ得る可変動弁装置を備え、
    前記燃焼室内ガス流速を少なくともエンジン回転速度と、前記吸気弁のリフト特性とに基づいて算出することを特徴とする請求項１３に記載のエンジンの点火時期制御装置。
16. 前記エンジンの時間損失を、燃焼室内ガスの燃焼速度に基づいて算出することを特徴とする請求項１１に記載のエンジンの点火時期制御装置。
17. 前記燃焼室内ガスの燃焼速度を、少なくともエンジン回転速度に基づいて算出することを特徴とする請求項１６に記載のエンジンの点火時期制御装置。
18. 吸気弁のリフト特性を変化させ得る可変動弁装置を備え、
    前記燃焼室内ガスの燃焼速度を、少なくともエンジン回転速度と前記吸気弁のリフト特性とに基づいて算出することを特徴とする請求項１６に記載のエンジンの点火時期制御装置。
19. 吸気弁のリフト特性を変化させ得る可変動弁装置を備え、
    前記燃焼室内ガスの燃焼速度を、燃焼室内ガスの温度、燃焼室ガスの圧力、不活性ガス率、エンジン回転速度及び前記吸気弁のリフト特性に基づいて算出することを特徴とする請求項１６に記載のエンジンの点火時期制御装置。
20. 燃焼室への吸気流動を変化させ得る可変吸気流動装置を備え、
    前記燃焼室内ガスの燃焼速度をこの可変吸気流動装置の状態に応じても算出することを特徴とする請求項１８または１９に記載のエンジンの点火時期制御装置。

Images