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