JP4692179B2 - 内燃機関の混合気濃度推定装置 - Google Patents

Description

本発明は、内燃機関の燃焼室内に噴射された燃料と燃焼室内に吸入されている筒内ガスとが混ざり合って形成される混合気の濃度（或いは、同濃度に相当する値）を推定する内燃機関の混合気濃度推定装置に関する。

火花点火式内燃機関、ディーゼル機関等の内燃機関から排出されるＮＯｘ等のエミッションの排出量は、着火後の火炎温度（燃焼温度）に強い相関を有している。また、この火炎温度は、燃焼室内の混合気の温度に大きく依存する。従って、ＮＯｘ等のエミッションの排出量を低減するためには、混合気温度を所定温度に制御することが有効である。

一般に、混合気温度は検出され得ないから、混合気温度を所定温度に制御するためには混合気温度を推定する必要がある。一方、混合気温度は、混合気の濃度（混合気の燃料濃度。以下、単に「混合気濃度」と云う。）に強い相関を有している。よって、混合気温度を推定するためには混合気濃度を推定することが効果的である。

また、特に、圧縮による自己着火により混合気が燃焼を開始するディーゼル機関においては、機関の運転状態に応じて着火時期を適切に制御する必要がある。この着火時期は混合気濃度に大きく依存する。従って、着火時期を適切に制御するためにも混合気濃度を推定する必要がある。

係る混合気濃度を推定する手法としては種々のものが知られている。例えば、下記非特許文献１には、運動量保存則を利用して混合気濃度（空気過剰率）を噴射後の経過時間の関数として求める式が紹介されていて、係る式に従って混合気濃度を求める手法が紹介されている。
機械学会論文集 25-156(1959年)，820頁 「ディーゼル機関の噴霧到達距離に関する研究」 和栗雄太郎，藤井勝，網谷竜夫，恒屋礼次郎

ところで、混合気が燃焼室内を移動していく過程では、実際には、ブレークアップ、蒸発等のような様々な現象が混合気に対して起こり得る。係る現象が発生すると、燃料と筒内ガスとが混ざり合う程度が変化し、以降、混合気濃度が大きく変化することが考えられる。

しかしながら、上記文献に紹介された運動量保存則を利用した手法では、係る現象が考慮されずに混合気濃度が算出されていく。従って、混合気濃度を精度良く推定することができないという問題がある。

本発明の目的は、混合気が燃焼室内を移動していく過程において混合気に対して起こり得る様々の現象を考慮して混合気濃度を精度良く推定できる混合気濃度推定装置を提供することにある。

本発明に係る混合気濃度推定装置は、内燃機関の燃焼室内に噴射された後において前記燃焼室内に吸入されている筒内ガス中にて時々刻々と拡散していく燃料と、前記燃料と混ざり合って混合気を形成する前記筒内ガスの一部分であってその量が前記燃料の拡散の進行に伴って時々刻々と増大していく混合気形成筒内ガスと、からなる混合気の濃度に相当する値である混合気濃度相当値の基本値を、前記燃料の噴射後において発生する前記混合気のブレークアップを考慮することなく且つ少なくとも運動量保存則を利用して、前記燃料の噴射後において逐次算出する混合気濃度算出手段を備えている。

混合気濃度算出手段は、例えば、上述した非特許文献１に記載された手法により運動量保存則を利用して混合気濃度を算出する。また、混合気濃度相当値とは、混合気の燃料濃度そのもの、或いは、空気過剰率等に代表される混合気濃度と等価な値を意味する。

本発明に係る混合気濃度推定装置の特徴は、前記混合気の状態に応じて前記混合気濃度算出手段により算出される混合気濃度相当値を補正する混合気濃度補正手段を更に備えたことにある。

これによれば、混合気が燃焼室内を移動していく過程において混合気に対して起こり得る様々の現象を考慮して混合気濃度算出手段により算出される混合気濃度相当値を補正することができる。従って、様々の現象を考慮して混合気濃度を精度良く推定できる。

なお、混合気濃度算出手段が、前記燃料と前記筒内ガスとが混ざり合う速度（混合気希釈速度）を算出し、算出された混合気希釈速度を積分していくことで混合気濃度相当値を算出するように構成されている場合、混合気濃度補正手段は、混合気濃度算出手段により算出された前記混合気希釈速度を補正して、補正された混合気希釈速度を積分していくことで混合気濃度相当値を補正してもよい。

より具体的には、上記本発明に係る混合気濃度推定装置においては、前記混合気のブレークアップが発生したか否かを前記燃料の噴射後において逐次判定するブレークアップ判定手段を更に備え、前記混合気濃度補正手段は、前記ブレークアップが発生したと判定された後において、前記混合気濃度算出手段により算出される混合気濃度相当値の基本値を前記ブレークアップに基づく分だけ前記混合気の濃度の減少方向に補正して最終的な混合気濃度相当値を逐次算出する。

混合気が燃焼室内を移動していく過程において混合気にブレークアップが発生すると、混合気の分散の程度がより一層進み易くなり、以降、ブレークアップが発生していない場合に比して混合気濃度が小さくなる（空気過剰率が大きくなる）傾向があると考えられる。しかしながら、例えば、上述した非特許文献１に記載された手法により運動量保存則を利用して混合気濃度を算出する場合、係るブレークアップが考慮されていない。

これに対し、上記のように、ブレークアップが発生したと判定された後において混合気濃度算出手段により算出される混合気濃度相当値を補正するように構成すれば、例えば、ブレークアップが発生したと判定された後において、上述した非特許文献１に記載された手法により算出される混合気濃度をより小さくなる方向（空気過剰率をより大きくなる方向）に補正していくことができる。従って、ブレークアップを考慮して混合気濃度を精度良く推定できる。

前記ブレークアップ判定手段は、例えば、前記燃焼室内を移動していく前記混合気の位置を算出する混合気位置算出手段と、前記燃焼室内において前記混合気のブレークアップの発生位置を算出するブレークアップ位置算出手段とを備え、前記算出された混合気の位置が前記算出されたブレークアップの発生位置に達したとき、前記ブレークアップが発生したと判定するように構成される。

ここにおいて、混合気位置算出手段は、例えば、上述した非特許文献１に記載された、燃焼室内を移動していく混合気の先頭部の位置を噴射後の経過時間の関数として求める式を利用して混合気位置を算出できる。また、ブレークアップ位置算出手段も、周知の一般式を利用してブレークアップ発生位置を算出できる。

また、上記本発明に係る混合気濃度推定装置においては、前記混合気内の（液滴）燃料が蒸発したか否かを前記燃料の噴射後において逐次判定する蒸発判定手段を更に備え、前記混合気濃度補正手段は、前記燃料が蒸発したと判定された後において、前記混合気濃度算出手段により算出される混合気濃度相当値の基本値を前記燃料の蒸発に基づく分だけ前記混合気の濃度の増大方向に補正して最終的な混合気濃度相当値を逐次算出するように構成されることが好適である。

混合気が燃焼室内を移動していく過程において混合気中の燃料（液滴燃料）が蒸発すると、混合気の粘度が大きくなって燃料（燃料蒸気）と筒内ガスとが混ざり合い難くなる。換言すれば、混合気中の燃料（液滴燃料）が蒸発した以降、蒸発が発生していない場合に比して混合気濃度が大きくなる（空気過剰率が小さくなる）傾向があると考えられる。しかしながら、例えば、上述した非特許文献１に記載された手法により運動量保存則を利用して混合気濃度を算出する場合、係る蒸発が考慮されていない。

これに対し、上記のように、液滴燃料が蒸発したと判定された後において混合気濃度算出手段により算出される混合気濃度相当値を補正するように構成すれば、例えば、液滴燃料が蒸発したと判定された後において、上述した非特許文献１に記載された手法により算出された混合気濃度をより大きくなる方向（空気過剰率をより小さくなる方向）に補正していくことができる。従って、液滴燃料の蒸発を考慮して混合気濃度を精度良く推定できる。

前記蒸発判定手段は、例えば、前記混合気の温度を算出する混合気温度算出手段を備え、前記算出された混合気の温度が前記燃料の蒸発温度に達したとき、前記燃料が蒸発したと判定するように構成される。ここにおいて、混合気温度算出手段は、例えば、上記推定される混合気濃度を利用して（或いは、更に、混合気内での種々の反応により発生する熱を考慮して）混合気温度を算出できる。

前記混合気濃度補正手段は、前記混合気の温度が高いほど、前記燃料の蒸発に基づく前記混合気の濃度の増大方向への補正の量を大きくするように構成されるとより好適である。混合気温度が高くなるほど、気体としての混合気の粘度が高くなって、燃料（燃料蒸気）と筒内ガスとがより混ざり合い難くなる。

上記構成によれば、混合気温度が高いほど、上述した非特許文献１に記載された手法により算出される混合気濃度を大きくなる方向（空気過剰率を小さくなる方向）に補正する程度を大きくすることができる。従って、液滴燃料の蒸発を考慮して混合気濃度をより一層精度良く推定できる。

また、上記本発明に係る混合気濃度推定装置においては、前記混合気濃度補正手段は、前記燃料が蒸発したと判定された後において、前記混合気濃度算出手段により算出される混合気濃度相当値の基本値を前記混合気内の反応により発生する熱による混合気の膨張に基づく分だけ前記混合気の濃度の減少方向に補正して最終的な混合気濃度相当値を逐次算出するように構成されることが好適である。

燃焼室内を移動していく混合気は、混合気内にて発生する燃焼に代表される種々の反応により発生する熱により熱膨張し得る。このように混合気が熱膨張すると、混合気に乱れが発生し、燃料（燃料蒸気）と筒内ガスとが混ざり合い易くなる。換言すれば、係る熱膨張が発生しない場合に比して混合気濃度が小さくなる（空気過剰率が大きくなる）傾向があると考えられる。しかしながら、例えば、上述した非特許文献１に記載された手法により運動量保存則を利用して混合気濃度を算出する場合、係る熱膨張が考慮されていない。

これに対し、上記のように、熱膨張を考慮して混合気濃度算出手段により算出される混合気濃度相当値を補正するように構成すれば、例えば、混合気の熱膨張の程度（例えば、混合気の熱膨張速度等）に応じて上述した非特許文献１に記載された手法により算出される混合気濃度をより小さくなる方向（空気過剰率をより大きくなる方向）に補正していくことができる。従って、混合気の熱膨張を考慮して混合気濃度を精度良く推定できる。

液滴燃料が蒸発した後では、混合気は気体のみから構成される。従って、これによれば、気体の状態方程式を利用して混合気の膨張の程度（膨張速度等）を算出することができる。

前記混合気濃度補正手段は、前記混合気の膨張の程度が大きいほど、前記混合気の膨張に基づく前記混合気の濃度の減少方向への補正の量を大きくするように構成されることがより好適である。混合気の膨張の程度（例えば、膨張速度等）が大きくなるほど、混合気の乱れの程度が大きくなって、燃料（燃料蒸気）と筒内ガスとが混ざり合い易くなる。

上記構成によれば、混合気の膨張の程度が大きいほど、上述した非特許文献１に記載された手法により算出される混合気濃度を小さくなる方向（空気過剰率を大きくなる方向）に補正する程度を大きくすることができる。従って、混合気の膨張を考慮して混合気濃度をより一層精度良く推定できる。

以下、本発明による内燃機関（ディーゼル機関）の混合気濃度推定装置の実施形態について図面を参照しつつ説明する。

図１は、本発明の実施形態に係る内燃機関の混合気濃度推定装置を含んだ混合気温度推定装置を４気筒内燃機関（ディーゼル機関）１０に適用したシステム全体の概略構成を示している。このシステムは、燃料供給系統を含むエンジン本体２０、エンジン本体２０の各気筒の燃焼室（筒内）にガスを導入するための吸気系統３０、エンジン本体２０からの排ガスを放出するための排気系統４０、排気還流を行うためのＥＧＲ装置５０、及び電気制御装置６０を含んでいる。

エンジン本体２０の各気筒の上部には燃料噴射弁（噴射弁、インジェクタ）２１が配設されている。各燃料噴射弁２１は、図示しない燃料タンクと接続された燃料噴射用ポンプ２２に燃料配管２３を介して接続されている。燃料噴射用ポンプ２２は、電気制御装置６０と電気的に接続されていて、電気制御装置６０からの駆動信号（後述する指令最終燃料噴射圧力Pcrfinに応じた指令信号）により燃料の実際の噴射圧力（吐出圧力）が指令最終燃料噴射圧力Pcrfinになるように同燃料を昇圧するようになっている。

これにより、燃料噴射弁２１には、燃料噴射用ポンプ２２から指令最終燃料噴射圧力Pcrfinまで昇圧された燃料が供給されるようになっている。また、燃料噴射弁２１は、電気制御装置６０と電気的に接続されていて、電気制御装置６０からの駆動信号（指令燃料噴射量（質量）Qfinに応じた指令信号）により噴射期間TAUだけ開弁し、これにより各気筒の燃焼室内に指令最終燃料噴射圧力Pcrfinにまで昇圧された燃料を指令燃料噴射量Qfinだけ直接噴射するようになっている。

吸気系統３０は、エンジン本体２０の各気筒の燃焼室にそれぞれ接続された吸気マニホールド３１、吸気マニホールド３１の上流側集合部に接続され吸気マニホールド３１とともに吸気通路を構成する吸気管３２、吸気管３２内に回動可能に保持されたスロットル弁３３、電気制御装置６０からの駆動信号に応答してスロットル弁３３を回転駆動するスロットル弁アクチュエータ３３ａ、スロットル弁３３の上流において吸気管３２に順に介装されたインタクーラー３４と過給機３５のコンプレッサ３５ａ、及び吸気管３２の先端部に配設されたエアクリーナ３６とを含んでいる。

排気系統４０は、エンジン本体２０の各気筒にそれぞれ接続された排気マニホールド４１、排気マニホールド４１の下流側集合部に接続された排気管４２、排気管４２に配設された過給機３５のタービン３５ｂ、及び排気管４２に介装されたディーゼルパティキュレートフィルタ（以下、「ＤＰＮＲ」と称呼する。）４３を含んでいる。排気マニホールド４１及び排気管４２は排気通路を構成している。

ＥＧＲ装置５０は、排気ガスを還流させる通路（ＥＧＲ通路）を構成する排気還流管５１と、排気還流管５１に介装されたＥＧＲ制御弁５２と、ＥＧＲクーラー５３とを備えている。排気還流管５１はタービン３５ｂの上流側排気通路（排気マニホールド４１）とスロットル弁３３の下流側吸気通路（吸気マニホールド３１）を連通している。ＥＧＲ制御弁５２は電気制御装置６０からの駆動信号に応答し、再循環される排気ガス量（排気還流量、ＥＧＲガス流量）を変更し得るようになっている。

電気制御装置６０は、互いにバスで接続されたＣＰＵ６１、ＣＰＵ６１が実行するプログラム、テーブル（ルックアップテーブル、マップ）、及び定数等を予め記憶したＲＯＭ６２、ＣＰＵ６１が必要に応じてデータを一時的に格納するＲＡＭ６３、電源が投入された状態でデータを格納するとともに同格納したデータを電源が遮断されている間も保持するバックアップＲＡＭ６４、並びにＡＤコンバータを含むインターフェース６５等からなるマイクロコンピュータである。

インターフェース６５は、吸気管３２に配置された熱線式エアフローメータ７１、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に設けられた吸気温センサ７２、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気管圧力センサ７３、クランクポジションセンサ７４、アクセル開度センサ７５、燃料噴射用ポンプ２２の吐出口の近傍の燃料配管２３に配設された燃料温度センサ７６、及び、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気酸素濃度センサ７７と接続されていて、これらのセンサからの信号をＣＰＵ６１に供給するようになっている。

また、インターフェース６５は、燃料噴射弁２１、燃料噴射用ポンプ２２、スロットル弁アクチュエータ３３ａ、及びＥＧＲ制御弁５２と接続されていて、ＣＰＵ６１の指示に応じてこれらに駆動信号を送出するようになっている。

熱線式エアフローメータ７１は、吸気通路内を通過する吸入空気の質量流量（単位時間当りの吸入空気量、単位時間あたりの新気量）を計測し、質量流量Ga（空気流量Ga）を表す信号を発生するようになっている。吸気温センサ７２は、エンジン１０のシリンダ（即ち、燃焼室、筒内）に吸入されるガスの温度（即ち、吸気温度）を検出し、吸気温度Tbを表す信号を発生するようになっている。吸気管圧力センサ７３は、エンジン１０のシリンダに吸入されるガスの圧力（即ち、吸気管圧力）を検出し、吸気管圧力Pbを表す信号を発生するようになっている。

クランクポジションセンサ７４は、各気筒の絶対クランク角度を検出し、実クランク角度CAactを表すとともにエンジン１０の回転速度であるエンジン回転速度NEをも表す信号を発生するようになっている。アクセル開度センサ７５は、アクセルペダルＡＰの操作量を検出し、アクセル操作量Accpを表す信号を発生するようになっている。燃料温度センサ７６は、燃料配管２３を通過する燃料の温度を検出し、燃料温度Tcrを表す信号を発生するようになっている。吸気酸素濃度センサ７７は、吸気中の酸素濃度を検出し、吸気酸素濃度RO2inを表す信号を発生するようになっている。

（混合気の空気過剰率の推定方法の概要）
次に、上記のように構成された混合気濃度推定装置を含んだ混合気温度推定装置（以下、「本装置」と云う。）による、混合気濃度相当値としての混合気の空気過剰率の推定方法について説明する。

図２は、或る一つの気筒のシリンダ内(筒内、燃焼室内)に吸気マニホールド３１からガスが吸入され、燃焼室内に吸入されたガスが排気マニホールド４１へ排出される様子を模式的に示した図である。図２に示したように、燃焼室内に吸入されるガス（従って、筒内ガス）には、吸気管３２の先端部からスロットル弁３３を介して吸入された新気と、排気還流管５１からＥＧＲ制御弁５２を介して吸入されたＥＧＲガスが含まれる。

吸入される新気量（新気質量）と吸入されるＥＧＲガス量（ガス質量）の和に対するＥＧＲガス量の割合（即ち、ＥＧＲ率）は、運転状態に応じて電気制御装置６０（ＣＰＵ６１）により適宜制御されるスロットル弁３３の開度、及びＥＧＲ制御弁５２の開度に応じて変化する。

かかる新気、及びＥＧＲガスは、吸気行程において開弁している吸気弁Vinを介してピストンの下降に伴って燃焼室内に吸入されて筒内ガスとなる。筒内ガスは、ピストンが圧縮下死点に達する時点近傍で吸気弁Vinが閉弁することにより燃焼室内に密閉され、その後の圧縮行程においてピストンの上昇に伴って圧縮される。

そして、ピストンが圧縮上死点近傍に達すると（具体的には、後述する最終燃料噴射開始時期（クランク角度）CAinjfinが到来すると）、本装置は、指令燃料噴射量Qfinに応じた噴射期間TAUだけ燃料噴射弁２１を開弁することで燃料を燃焼室内に直接噴射する。この結果、燃料噴射弁２１の噴孔から噴射された液体の（液滴）燃料は、時間の経過に伴って筒内ガスを取り込みながら混合気となって燃焼室内において円錐状に拡散していく。

ここで、指令燃料噴射量Qfinの燃料は、実際には、最終燃料噴射開始時期CAinjfinから噴射期間TAUだけ連続して噴射されるが、以下においては便宜上、指令燃料噴射量Qfinの燃料が最終燃料噴射開始時期CAinjfinにて一時（瞬時）に噴射されるものとして説明を続ける。

図３（ａ）は、燃料噴射弁２１の噴孔から、指令燃料噴射量（質量）Qfinの燃料が一時に噴射された時点（即ち、噴射後経過時間ｔ＝０）での質量Qfinの液滴燃料の様子を模式的に示した図である。図３（ｂ）は、その後の或る時点（任意の噴射後経過時間ｔ）での図３（ａ）に示した質量Qfinの燃料の様子を模式的に示した図である。

図３（ｂ）に示すように、質量Qfinの燃料は、最終燃料噴射開始時期CAinjfin（即ち、噴射後経過時間ｔ＝０）において噴射された後、噴霧角θをもって円錐状に拡散しながら筒内ガスを順次取り込んでいく。そして、質量Qfinの燃料は、任意の噴射後経過時間ｔにおいて、噴射後経過時間ｔの関数である質量Gの筒内ガス（以下、「混合気形成筒内ガス」と云うこともある。）と混ざり合って質量（Qfin＋G）の混合気となっているものと仮定する。本装置は、係る混合気の任意の噴射後経過時間ｔにおける混合気濃度相当値としての空気過剰率λを噴射後経過時間ｔ＝０から微小時間Δｔ毎に推定するものである。

噴射後経過時間ｔにおける空気過剰率λは、下記(1)式に示すように定義される。下記(1)式において、stoichは、単位質量の燃料の燃焼に必要な筒内ガスの質量（以下、「筒内ガス理論空燃比stoich」と呼ぶ。）である。筒内ガス理論空燃比stoichの値は、吸気中の酸素濃度に応じて変化すると考えられるから、上記吸気酸素濃度RO2inを引数とする所定の関数に従って取得され得る。また、上述したように、Gは混合気形成筒内ガスの質量であり、Qfinは指令燃料噴射量（質量）である（図３を参照）。

このように定義される空気過剰率λは、例えば、上述した非特許文献１にて紹介された、運動量保存則を利用して得られる式である下記(2)式、及び下記(3)式に基づいて噴射後経過時間ｔの関数として求めることができる。

上記(3)式において、ｔは上記噴射後経過時間であり、(dλ/dt)baseは噴射後経過時間ｔの関数である混合気希釈速度（基本混合気希釈速度）である。また、ｃは収縮係数、ｄは燃料噴射弁２１の噴孔径、ρfは液滴燃料の密度、Ｌは論理希釈ガス量であって、これらの各値は全て定数である。

上記(3)式において、ΔPは有効噴射圧力であって、上記最終燃料噴射圧力Pcrfinから噴射開始時点（即ち、噴射後経過時間ｔ＝０）での筒内ガス圧力Pg0を減じた値である。筒内ガス圧力Pg0は、圧縮行程（及び膨張行程）における筒内ガスの状態が吸気弁Vinの閉弁時（即ち、筒内ガスが密閉された時点。以下、「ＩＶＣ」と呼ぶ。）以降、断熱変化するとの仮定のもと、下記(4)式に従って求めることができる。

上記(4)式において、PgivcはＩＶＣにおける筒内ガス圧力である。上述したように、ＩＶＣは圧縮下死点近傍であるから、ＩＶＣにおいて筒内ガス圧力は吸気管圧力Pbと略等しいと考えられる。従って、ＩＶＣにおいて吸気管圧力センサ７３により検出される吸気管圧力PbがPgivcとして使用され得る。Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積であり、Vg(CAinj)は噴射後経過時間ｔ＝０におけるクランク角度CAに対応する筒内容積である。筒内容積Vgは機関１０の設計諸元に基づいてクランク角度CAの関数Vg(CA)として取得することができるから、Vg(Caivc)，Vg(CAinj)も取得することができる。κは筒内ガスの比熱比（本例では、一定）である。

また、上記(3)式において、θは図３（ｂ）に示した噴霧角である。噴霧角θは、噴射開始時点（即ち、噴射後経過時間ｔ＝０）における筒内ガスの密度ρg0、及び上記有効噴射圧力ΔPに応じて変化すると考えられるから、筒内ガスの密度ρg0、及び有効噴射圧力ΔPと噴霧角θとの関係を予め規定したテーブルMapθに基づいて取得することができる。筒内ガスの密度ρg0は、筒内ガスの全質量Mgを、噴射後経過時間ｔ＝０における上記筒内容積Vg(CAinj)で除することで取得することができる。筒内ガスの全質量Mgは、ＩＶＣにおける気体の状態方程式に基づく下記(5)式に従って取得され得る。下記(5)式において、TgivcはＩＶＣにおける筒内ガス温度である。ＩＶＣは圧縮下死点近傍であるから、ＩＶＣにおいて筒内ガス温度は吸気温度Tbと略等しいと考えられる。従って、ＩＶＣにおいて吸気温センサ７２により検出される吸気温度TbがTgivcとして使用され得る。Rは筒内ガスのガス定数（本例では、一定）である。

また、上記(3)式において、ρgは噴射後経過時間ｔにおける筒内ガス密度であって、前記筒内ガスの全質量Mgを、噴射後経過時間ｔにおける上記筒内容積Vg(CA)で除することで、噴射後経過時間ｔの関数として取得することができる。

以上のことから、有効噴射圧力ΔPと噴霧角θとを上述のようにして求めれば、噴射後経過時間ｔの値と同噴射後経過時間ｔの関数である筒内ガス密度ρgの値とを使用して、上記(3)式に従って基本混合気希釈速度(dλ/dt)baseが噴射後経過時間ｔの関数として求められる。そして、噴射後経過時間ｔ＝０から微小時間Δｔ（例えば、０．１msec）毎に求めた基本混合気希釈速度(dλ/dt)baseの値を上記(2)式に従って時間で積分（積算）していくことで噴射後経過時間ｔにおける空気過剰率λ（後述する補正がなされていない空気過剰率）を噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。

なお、上記(3)式から取得される基本混合気希釈速度(dλ/dt)baseの値は常に正の値となることから上記(2)式から取得される空気過剰率λの値は噴射後経過時間ｔの増大に従って増加していく。そうすると、上記(1)式から理解できるように、混合気形成筒内ガスの質量Gが噴射後経過時間ｔの増大に従って増加していく。このことは、噴射後の燃料蒸気が円錐状に拡散していくことに伴って燃料蒸気と混ざり合う（燃料蒸気が取り込む）筒内ガス（従って、混合気形成筒内ガス）の量が増大していくことに対応している。

（空気過剰率（混合気希釈速度）の補正）
混合気が燃焼室内を移動していく過程では、実際には、燃料と筒内ガスとが混ざり合う程度を変化させるような様々な現象が混合気に対して起こり得る。係る現象が発生すると、混合気希釈速度（従って、空気過剰率λ）が大きく変化することが考えられる。本例では、係る現象として、以下のように、ブレークアップ、蒸発、及び熱膨張が考慮される。

<ブレークアップを考慮した空気過剰率（混合気希釈速度）の補正>
混合気が燃焼室内を移動していく過程において混合気にブレークアップが発生することが知られている。このブレークアップが発生すると、混合気内の液滴燃料の運動エネルギーの一部が内部エネルギーに変換されること等に起因して、混合気の分散の程度がより一層進み易くなる。即ち、ブレークアップが発生すると、液滴燃料が筒内ガスを取り込み易くなり、以降、ブレークアップが発生していない場合に比して混合気希釈速度（従って、空気過剰率λ）が大きくなる傾向があると考えられる。

しかしながら、上記(3)式に従って運動量保存則を利用して算出される基本混合気希釈速度(dλ/dt)baseは、係るブレークアップが全く考慮されずに算出される値である。従って、ブレークアップを考慮すると、ブレークアップが発生した時点以降、上記(3)式に従って算出される基本混合気希釈速度(dλ/dt)baseの値を大きくなる方向に補正する必要がある。

他方、混合気のブレークアップが発生する位置（燃料噴射弁２１の噴孔からの距離。以下、「ブレークアップ発生位置Xbreak」と称呼する。）は、例えば、周知の下記(6)式に従って求めることができる。

上記(6)式において、C0は燃料噴射弁２１の形状定数（ノズル形状定数）である。また、上記(3)式と同様、ｄは燃料噴射弁２１の噴孔径であり、ρfは液滴燃料の密度であり、ρgは噴射後経過時間ｔにおける筒内ガス密度である。このように、ブレークアップ発生位置Xbreakは、上記(6)式に従って筒内ガス密度ρgの関数（従って、噴射後経過時間ｔの関数）として求められる。

また、噴射開始時点以降における混合気（先頭部）の燃料噴射弁２１の噴孔からの到達距離（以下、「混合気到達距離X」と称呼する。）は、例えば、上記非特許文献１にて紹介された式である下記(7)式、(8)式に基づいて噴射後経過時間ｔの関数として求めることができる。下記(8)式において、dX/dtは噴射後経過時間ｔの関数である混合気移動速度である。なお、下記(8)式の右辺に示される各種値は、上記(3)式の右辺に示されるものと同一である。

即ち、噴射後経過時間ｔの値と噴射後経過時間ｔの関数である筒内ガス密度ρgの値とにより、上記(8)式に従って混合気移動速度dX/dtを逐次求めていき、逐次求めた混合気移動速度dX/dtの値を上記(7)式に従って時間で積分していくことで噴射後経過時間ｔにおける混合気到達距離Xを取得することができる。

そこで、本装置は、上記(7)式に従って算出される混合気到達距離Xが上記(6)式に従って算出されるブレークアップ発生位置Xbreakに達した時点で混合気にブレークアップが発生したと判定する。本装置は、希釈係数K1を導入し、ブレークアップが発生したと判定する前の段階では希釈係数K1の値を「１」に設定し、ブレークアップが発生したと判定した時点以降、希釈係数K1の値を一定値α（＞１）に設定する。

そして、本装置は、上記基本混合気希釈速度(dλ/dt)baseの値に希釈係数K1を乗じることで基本混合気希釈速度(dλ/dt)baseを補正する。これにより、ブレークアップが発生した時点以降、係るブレークアップの観点から、上記基本混合気希釈速度(dλ/dt)baseの値が大きくなる方向に補正されていく。

<蒸発を考慮した空気過剰率（混合気希釈速度）の補正>
混合気が燃焼室内を移動していく過程において混合気の温度が燃料の蒸発温度に達すると、混合気内の液滴燃料が蒸発する。液滴燃料が蒸発すると、混合気の粘度が大きくなり、燃料（燃料蒸気）が筒内ガスを取り込み難くなる。即ち、液滴燃料が蒸発した時点以降、係る蒸発が発生していない場合に比して混合気希釈速度（従って、空気過剰率λ）が小さくなる傾向があると考えられる。

しかしながら、上記(3)式に従って運動量保存則を利用して算出される基本混合気希釈速度(dλ/dt)baseは、係る液滴燃料の蒸発が全く考慮されずに算出される値である。従って、液滴燃料の蒸発を考慮すると、液滴燃料が蒸発した時点以降、上記(3)式に従って算出される基本混合気希釈速度(dλ/dt)baseの値を小さくなる方向に補正する必要がある。

更には、液滴燃料が蒸発した時点以降、混合気の温度が高くなるほど、気体としての混合気の粘度が高くなって、燃料（燃料蒸気）と筒内ガスとがより混ざり合い難くなる。換言すれば、混合気の温度が高くなるほど、上記基本混合気希釈速度(dλ/dt)baseの値を小さくなる方向に補正する程度を大きくする必要がある。他方、混合気の温度（混合気温度Tmix）は、後述するように、噴射後経過時間ｔの関数として求めることができる。

そこで、本装置は、後述する混合気温度Tmixが燃料の蒸発温度に達した時点で混合気内の液滴燃料が蒸発したと判定する。本装置は、希釈係数K2を導入し、液滴燃料が蒸発したと判定する前の段階では希釈係数K2の値を「１」に設定し、液滴燃料が蒸発したと判定した時点以降、希釈係数K2の値を、図４に示したように、「１」未満の値であって、混合気温度Tmixが高いほどより小さい値に設定する。

そして、本装置は、上記基本混合気希釈速度(dλ/dt)baseの値に希釈係数K2を乗じることで基本混合気希釈速度(dλ/dt)baseを補正する。これにより、液滴燃料が蒸発した時点以降、係る蒸発の観点から、上記基本混合気希釈速度(dλ/dt)baseの値が小さくなる方向に補正されていく。

<熱膨張を考慮した空気過剰率（混合気希釈速度）の補正>
混合気内の液滴燃料が蒸発した時点以降、気体として燃焼室内を移動していく混合気は、混合気内にて発生する燃焼に代表される種々の反応により発生する熱により熱膨張し得る。このように混合気が熱膨張すると、混合気と混合気を取り巻く筒内ガスとの相対速度が大きくなり、混合気における筒内ガスとの境界部分に乱れが発生し易くなる。このことは、燃料（燃料蒸気）と筒内ガスとが混ざり合い易くなることを意味する。即ち、液滴燃料が蒸発した時点以降、係る熱膨張が発生しない場合に比して混合気希釈速度（従って、空気過剰率λ）が大きくなる傾向があると考えられる。

しかしながら、上記(3)式に従って運動量保存則を利用して算出される基本混合気希釈速度(dλ/dt)baseは、係る熱膨張が全く考慮されずに算出される値である。従って、熱膨張を考慮すると、液滴燃料が蒸発した時点以降、上記(3)式に従って算出される基本混合気希釈速度(dλ/dt)baseの値を大きくなる方向に補正する必要がある。

更には、混合気の熱膨張の速度（膨張速度）が大きくなるほど、上述した混合気の乱れの程度が大きくなって、燃料（燃料蒸気）と筒内ガスとがより混ざり合い易くなる。換言すれば、混合気の膨張速度が大きくなるほど、上記基本混合気希釈速度(dλ/dt)baseの値を大きくなる方向に補正する程度を大きくする必要がある。他方、図５に示すように、気体としての混合気が球形状を呈していると仮定すれば、混合気の膨張速度dr/dtは、下記(9)式にて表すことができる。

上記(9)式において、r0は、混合気内の液滴燃料が蒸発した時点以降の任意の噴射後経過時間t＝tにおける球形状の混合気の半径である。この半径r0は、気体の状態方程式に基づく下記(10)式から得られる球形状の混合気の体積Vmixを使用して、下記(11)式に従って求めることができる。

上記(10)式において、Mmixは混合気の総質量（混合気質量）、Tmixは上記混合気温度、Pgは筒内ガス圧力、Rは混合気のガス定数（本例では、一定）である。混合気質量Mmix、混合気温度Tmix、及び筒内ガス圧力Pgは、それぞれ後述するように噴射後経過時間ｔ＝０から微小時間Δｔ毎に求めることができる。従って、半径r0も、微小時間Δｔ毎に求めることができる。

他方、上記(9)式において、rは、噴射後経過時間ｔ＝ｔからの微小時間Δｔの間において混合気内にて発生する化学反応による反応熱Hr（詳細は、後述する。）により熱膨張した混合気の噴射後経過時間ｔ＝(ｔ＋Δｔ)における半径（r＞r0）である。この微小時間Δｔの間において混合気内にて発生する反応熱Hrにより混合気温度TmixがΔTだけ上昇したものとすると、熱膨張後の混合気の半径rは、微小時間Δｔに亘って筒内ガス圧力Pg、及び混合気質量Mmixが一定であるとの仮定の下、気体の状態方程式を利用して、下記(12)式に従って求めることができる。

上記(12)式において、反応熱Hrによる混合気の温度上昇量ΔTは、下記(13)式に従って求めることができる。下記(13)式において、Cmixは混合気の定圧比熱であり、後述するように噴射後経過時間ｔ＝０から微小時間Δｔ毎に求めることができる。従って、半径rも、微小時間Δｔ毎に求めることができる。

以上より、混合気内の液滴燃料が蒸発した時点以降、上記(9)式にて表される混合気の膨張速度dr/dtも微小時間Δｔ毎に求めていくことができる。

そこで、本装置は、希釈係数K3を導入し、液滴燃料が蒸発したと判定する前の段階では希釈係数K3の値を「１」に設定し、液滴燃料が蒸発したと判定した時点以降、希釈係数K3の値を、図６に示したように、「１」より大きい値であって、上記(9)式に従って求められる混合気の膨張速度dr/dtに係数Kr（一定）を乗じた値（以下、「乱れ強度Ｂ」と称呼する。）が大きいほどより大きい値に設定する。

そして、本装置は、上記基本混合気希釈速度(dλ/dt)baseの値に希釈係数K3を乗じることで基本混合気希釈速度(dλ/dt)baseを補正する。これにより、液滴燃料が蒸発した時点以降、混合気の熱膨張の観点から、上記基本混合気希釈速度(dλ/dt)baseの値が大きくなる方向に補正されていく。

以上のように、本装置は、実際には、下記(14)式に従って、上記基本混合気希釈速度(dλ/dt)baseの値にブレークアップ、蒸発、及び熱膨張にそれぞれ対応する上記希釈係数K1,K2,K3を乗じた値である補正後混合気希釈速度dλ/dtを、噴射後経過時間ｔ＝０から微小時間Δｔ（例えば、０．１msec）毎に求め、それぞれの補正後混合気希釈速度dλ/dtの値を下記(15)式に従って時間で積分（積算）していくことで噴射後経過時間ｔにおける空気過剰率λ（補正後の空気過剰率）を噴射後経過時間ｔ＝０から微小時間Δｔ毎に求める。

図７は、本装置により上記(14)、(15)式に従って求められた噴射後経過時間ｔに対する補正後の空気過剰率λの変化の一例を示したグラフである。この例では、噴射後経過時間ｔ＝ｔ１にて混合気のブレークアップが発生し、その後のｔ＝ｔ２にて混合気内の液滴燃料が蒸発し、その後のｔ＝ｔ３にて着火（燃焼。熱炎）が発生した場合が示されている。なお、図７における破線は、上記(1)、(2)式に従って求められる補正前の空気過剰率λbaseの変化を示している。

図７に示すように、この場合、ｔ＝ｔ１以前では、希釈係数K1,K2,K3の全てが「１」に設定される。これにより、補正後混合気希釈速度dλ/dtが基本混合気希釈速度(dλ/dt)baseと一致するから、補正後の空気過剰率λも補正前の空気過剰率λbaseと一致する。

混合気のブレークアップが発生するｔ＝ｔ１以降、希釈係数K1がα(＞１、一定値)に設定される。混合気内の液滴燃料の蒸発が発生するｔ＝ｔ２以降、希釈係数K2が「１」未満の値（図４を参照）に設定されるとともに希釈係数K3が「１」より大きい値（図６を参照）に設定される。この結果、ｔ＝ｔ１以降、補正後混合気希釈速度dλ/dtが基本混合気希釈速度(dλ/dt)baseと一致しなくなるから、補正後の空気過剰率λも補正前の空気過剰率λbaseから偏移していく。

なお、着火が発生するｔ＝ｔ３以降の短期間に亘って、補正後の空気過剰率λが急激に増大している。これは、着火の発生により上記反応熱Hrが一時的に非常に大きい値になって上記温度上昇量ΔT（上記(13)式を参照）が一時的に非常に大きい値になることで、膨張速度dr/dt（従って、乱れ強度Ｂ）、ひいては希釈係数K3が一時的に非常に大きい値になることに基づく。

なお、図７では、混合気のブレークアップが発生した後に混合気内の液滴燃料が蒸発する例が示されているが、機関の状態によっては、混合気内の液滴燃料が蒸発した後に混合気のブレークアップが発生する場合も発生し得る。以上が、空気過剰率（混合気希釈速度）の補正の概要である。

（混合気温度Tmixの取得）
次に、上述のように噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得される補正後の空気過剰率λの値を利用して任意の噴射後経過時間ｔにおける混合気温度Tmixを取得する方法について説明する。一般に、混合気の熱エネルギー（エンタルピ）Hmixは、混合気温度Tmixを用いて下記(16)式に従って表すことができる。

上記(16)式において、Mmixは上述した混合気の総質量（混合気質量）、Cmixは上述した混合気の定圧比熱である。従って、混合気のエンタルピHmix、混合気質量Mmix、及び混合気の定圧比熱Cmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていく（更新していく）ことで、下記(17)式に従って混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。以下、先ず、混合気質量Mmixの求め方について説明する。

<混合気質量Mmix>
上述したように、質量Qfinの燃料蒸気は、任意の噴射後経過時間ｔにおいて質量Gの混合気形成筒内ガスと混ざり合って質量（Qfin＋G）の混合気となっているから、任意の噴射後経過時間ｔにおける混合気質量Mmixは(Qfin＋G)である。ここで、上記(1)式より「G＝stoich・λ・Qfin」と表すことができるから、混合気質量Mmixは、空気過剰率λを用いて下記(18)式にて表すことができる。

よって、上述したように噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得され得る空気過剰率λの値を上記(18)式に順次適用していくことで、混合気質量Mmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。

<混合気の定圧比熱Cmix>
次に、混合気の定圧比熱Cmixの求め方について説明する。一般に、混合気の定圧比熱Cmixは、同混合気内の酸素濃度[O2]mix、及び混合気温度Tmixに大きく依存すると考えられる。ここで、混合気内の酸素濃度[O2]mixは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。従って、混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができれば、下記(19)式に従って、混合気の定圧比熱Cmixを微小時間Δｔ毎に求めることができる。

上記(19)式において、funcCmixは、混合気の酸素濃度[O2]mix、及び混合気温度Tmixを引数とする混合気の定圧比熱Cmixを求めるための関数である。なお、上記(19)式を使用して混合気の定圧比熱Cmixを微小時間Δｔ毎に求めていく際における[O2]mix，Tmixの引数値としては、それぞれ現時点（即ち、噴射後経過時間ｔ）よりも微小時間Δｔ前の値が使用される。

<混合気のエンタルピHmix>
次に、混合気のエンタルピHmixの求め方について説明する。いま、噴射後経過時間（ｔ−Δt）における混合気のエンタルピHmix(ｔ−Δｔ)が既知である場合において、同噴射後経過時間(ｔ−Δｔ)から噴射後経過時間ｔまでの微小時間Δｔの間における混合気のエンタルピの増加分ΔHmixについて考える。この混合気のエンタルピの増加分ΔHmixとしては、微小時間Δｔの間に混合気に新たに取り込まれる筒内ガスの熱エネルギーΔHgと、同微小時間Δｔの間において混合気内で発生する化学反応による上記反応熱Hrとが挙げられる。

先ず、上記筒内ガスの熱エネルギーΔHgは、下記(20)式により表すことができる。ここにおいて、gは微小時間Δｔの間に混合気に新たに取り込まれる筒内ガスの質量である。この質量gは、噴射後経過時間ｔにおける混合気形成筒内ガスの質量から噴射後経過時間（ｔ−Δｔ）における混合気形成筒内ガスの質量を減じた値である。従って、上述した関係「G＝stoich・λ・Qfin」を利用して下記(21)式により求めることができる。(21)式において、λ(ｔ)，λ(ｔ−Δｔ)はそれぞれ、噴射後経過時間ｔ，(ｔ−Δｔ)における空気過剰率であり、上記(14)式、(15)式から求めることができる。

また、上記(20)式において、Tgは、噴射後経過時間ｔにおける筒内ガスの温度であり、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(22)式に従って求めることができる。下記(22)式において、上述したように、TgivcはＩＶＣにおける筒内ガス温度であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。

また、上記(20)式において、Cgは、噴射後経過時間ｔにおける筒内ガスの定圧比熱であり、混合気の定圧比熱Cmixを求める上記(19)式と同様、下記(23)式に従って求めることができる。下記(23)式において、funcCgは、吸気中の酸素濃度[O2]in、及び筒内ガス温度Tgを引数とする筒内ガスの定圧比熱Cgを求めるための関数である。

なお、上記(23)式を使用して筒内ガスの定圧比熱Cgを微小時間Δｔ毎に求めていく際における[O2]inの引数値としては、吸気酸素濃度センサ７７により検出される上記吸気酸素濃度RO2inが使用される。また、上記筒内ガス温度Tgの引数値としては、現時点（即ち、噴射後経過時間ｔ）における値が使用される。以上により、上記(20)式の右辺の項の全てを求めることができるから、同(20)式に従って上記筒内ガスの熱エネルギーΔHgを求めることができる。

次に、上記微小時間Δｔの間において混合気内で発生する化学反応による反応熱Hrは、下記(24)式で表すことができる。下記(24)式において、Hfは所定の定数であり、qrは上記微小時間Δｔの間において混合気内で発生する化学反応による燃料消費量である。

上記燃料消費量qrの対象となる化学反応としては、前駆物質Ｐの生成反応、着火反応（熱炎反応）、低温酸化反応（冷炎反応）のみならず、その他の種々の化学反応が含まれる。この燃料消費量qrは、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、混合気温度Tmixに大きく依存すると考えられるから、下記(25)式に従って表すことができる。

上記(25)式において、funcqrは、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、及び混合気温度Tmixを引数とする上記燃料消費量qrを求めるための関数である。混合気内の燃料濃度[Fuel]mixも、上記混合気内の酸素濃度[O2]mixと同様、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。上記(25)式を使用して上記燃料消費量qrを微小時間Δｔ毎に求めていく際における[O2]mix，[Fuel]mixの引数値としては、ぞれぞれ現時点（即ち、噴射後経過時間ｔ）よりも微小時間Δｔ前の値が使用される。

また、上記(25)式の混合気温度Tmixの引数値としては化学反応前の混合気温度Tpreが使用される。この化学反応前の混合気温度Tpreは、上記(21)式にて算出される質量gの筒内ガスが混合気に新たに取り込まれた後であって、且つ、噴射後経過時間(ｔ−Δｔ)からの上記微小時間Δｔの間における化学反応が発生する前の段階における混合気温度であり、下記(26)式に従って求めることができる。

上記(26)式において、Mmix(ｔ−Δｔ)，Cmix(ｔ−Δｔ)はそれぞれ、噴射後経過時間(ｔ−Δｔ)における混合気質量、及び混合気の定圧比熱であって、上記(18)式、及び上記(19)式によりそれぞれ取得され得る。更には、噴射後経過時間(ｔ−Δｔ)における混合気のエンタルピHmix(ｔ−Δｔ)は既知である。以上により、上記化学反応前の混合気温度Tpreを求めることができる。従って、上記(25)式の右辺の引数値を全て取得することができるから、(25)式，(24)式に従って上記化学反応による反応熱Hrを求めることができる。

以上のことから、噴射後経過時間（ｔ−Δt）における混合気のエンタルピHmix(ｔ−Δｔ)が既知である場合において、同噴射後経過時間(ｔ−Δｔ)から噴射後経過時間ｔまでの微小時間Δｔの間における混合気のエンタルピの増加分ΔHmix(＝ΔHg＋Hr)が求められるから、噴射後経過時間ｔにおける混合気のエンタルピHmix(ｔ)（＝Hmix(ｔ−Δｔ)＋ΔHmix）を求めることができる。

更には、噴射後経過時間ｔ＝０における混合気は、筒内ガスが取り込まれる前の状態（即ち、液滴燃料のみ）であるから（図３（ａ）を参照）、この時点での混合気のエンタルピHmix(0)は、下記(27)式にて求めることができる。ここにおいて、Cfは燃料（蒸気）の定圧比熱（ここでは、定数）である。

また、Tfは、燃料蒸気そのものの温度であり、液体の燃料が噴射直後に燃料蒸気に変化する際の単位質量当たりの潜熱Qvaporを考慮して下記(28)式に従って求めることができる。下記(28)式において、Tcrは噴射後経過時間ｔ＝０において燃料温度センサ７６により検出される液体の燃料温度である。αcrは燃料が燃料噴射用ポンプ２２の吐出口近傍から燃料噴射弁２１までの燃料配管２３を通過する際の熱損失分を考慮するための補正係数である。

従って、噴射後経過時間ｔ＝０における混合気のエンタルピHmix(0)も求めることができる。以上より、混合気のエンタルピHmixを、噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。

このようにして、混合気のエンタルピHmix、混合気質量Mmix、及び混合気の定圧比熱Cmixがそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求められるから、上記(17)式に従って混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。

<混合気内の燃料濃度[Fuel]mixの取得>
次に、混合気内の燃料濃度[Fuel]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の燃料濃度[Fuel]mixは、上記(18)式により取得される噴射後経過時間ｔにおける混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する燃料の質量」の割合である。

ここで、「噴射後経過時間ｔにおいて混合気内に存在する燃料の質量」は、噴射後経過時間ｔ＝０にて噴射された燃料量（指令燃料噴射量Qfin）から、噴射から現時点（噴射後経過時間ｔ）までの間に化学反応で消費された燃料分を減じた値である。従って、噴射後経過時間ｔにおける混合気内の燃料濃度[Fuel]mixは、下記(29)式で表すことができる。

上記(29)式において、「Σqr」は、噴射から現時点（噴射後経過時間ｔ）まで微小時間Δｔ毎に上記(25)式に従って逐次取得・更新されていくそれぞれの燃料消費量qrの和である。このように、燃料消費量qr、及び混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(29)式に従って混合気内の燃料濃度[Fuel]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。

<混合気内の酸素濃度[O2]mixの取得>
次に、混合気内の酸素濃度（質量濃度）[O2]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の酸素濃度[O2]mixは、噴射後経過時間ｔにおける上記混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する酸素の質量」の割合である。

ここで、「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」は、上述した噴射後経過時間ｔにおける混合気形成筒内ガスの質量Gから、噴射から現時点（噴射後経過時間ｔ）までの間に化学反応で消費された筒内ガス分を減じた値である。上記微小時間Δｔにおける燃料消費量qrの燃料と反応して同微小時間Δｔの間において混合気内で化学反応により消費される筒内ガスの消費量grは、下記(30)式にて表すことができる。

従って、「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」は、「G−Σgr」と表すことができる。ここで、「Σgr」は、噴射から現時点（噴射後経過時間ｔ）まで微小時間Δｔ毎に上記(30)式に従って逐次取得・更新されていくそれぞれの筒内ガス消費量grの和である。

「噴射後経過時間ｔにおいて混合気内に存在する酸素の質量」は、上記「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」に、筒内ガス内の酸素濃度（従って、上記吸気中の酸素濃度[O2]in）を乗じることで求めることができる。以上のことから、噴射後経過時間ｔにおける混合気内の酸素濃度[O2]mixは、下記(31)式で表すことができる。

このように、混合気形成筒内ガスの質量G、筒内ガス消費量gr、及び混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(31)式に従って混合気内の酸素濃度[O2]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。

<筒内ガス圧力Pgの取得>
次に、筒内ガス圧力Pgを取得する方法について説明する。筒内ガス圧力Pgは、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(32)式に従って求めることができる。下記(32)式において、上述したように、PgivcはＩＶＣにおける筒内ガス圧力であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。以上、混合気温度Tmixの推定に必要となる各種物理量の推定方法について説明した。以上が、上述した補正後の空気過剰率λの値を利用した、任意の噴射後経過時間ｔにおける混合気温度Tmixの取得方法の概要である。

（燃料噴射制御の概要）
本装置は、目標着火時期を所定時期（CAref、例えば、ATDC10°）に設定し、混合気（先頭部）の着火時期を目標着火時期CArefに一致させるため、クランク角度CAがCArefとなる時点において取得された混合気温度Tmix（以下、「制御用混合気温度」と称呼する。）が所定の目標混合気温度Tmixrefになるように、燃料噴射開始時期、及び燃料噴射圧力をフィードバック制御する。

具体的には、本装置は、上述した空気過剰率λ、混合気温度Tmix等の計算を筒内ガスの量が確定するＩＶＣの直後に開始し、燃料噴射開始時期の到来前に上記制御用混合気温度の推定を完了する。

そして、本装置は、上記制御用混合気温度が目標混合気温度Tmixrefよりも高いとき、今回の燃料噴射気筒についての燃料噴射開始時期を基本燃料噴射開始時期よりも所定量だけ遅らせ、且つ、燃料噴射圧力を基本燃料噴射圧力よりも所定量だけ低くする。これにより、今回の燃料噴射気筒についての目標着火時期CArefにおける実際の混合気温度が低くなるように制御され、この結果、今回の燃料噴射気筒の実際の着火時期が目標着火時期CArefに一致せしめられる。

一方、上記制御用混合気温度が目標混合気温度Tmixrefよりも低いとき、今回の燃料噴射気筒についての燃料噴射開始時期を基本燃料噴射開始時期よりも所定量だけ早め、且つ、燃料噴射圧力を基本燃料噴射圧力よりも所定量だけ高くする。これにより、今回の燃料噴射気筒についての目標着火時期CArefにおける実際の混合気温度が高くなるように制御され、この結果、今回の燃料噴射気筒の実際の着火時期が目標着火時期CArefに一致せしめられる。以上が、燃料噴射制御の概要である。

（実際の作動）
次に、上記のように構成された内燃機関の混合気濃度推定装置を含んだ混合気温度推定装置の実際の作動について説明する。

<空気過剰率、及び混合気温度の算出>
ＣＰＵ６１は、図８〜図１２に一連のフローチャートにより示した空気過剰率、及び混合気温度の算出を行うためのルーチンを所定時間の経過毎に、気筒毎に、繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１はステップ８００から処理を開始し、ステップ８０５に進んで吸気弁Vinが開状態から閉状態へと変化したか否か（即ち、ＩＶＣが到来したか否か）を判定し、「Ｎｏ」と判定する場合、ステップ８９５に直ちに進んで本ルーチンを一旦終了する。

いま、或る気筒においてＩＶＣが到来したものとすると、ＣＰＵ６１はステップ８０５に進んだとき「Ｙｅｓ」と判定してステップ８１０に進み、ＩＶＣ時クランク角度CAivcをクランクポジションセンサ７４から取得される現時点での実クランク角度CAactの値に設定し、ＩＶＣ時筒内ガス圧力Pgivcを吸気管圧力センサ７３から得られる現時点での吸気管圧力Pbの値に設定し、ＩＶＣ時筒内ガス温度Tgivcを吸気温センサ７２から得られる現時点での吸気温度Tbの値に設定し、吸気酸素濃度[O2]inを吸気酸素濃度センサ７７から得られる現時点での吸気酸素濃度RO2inの値に設定する。

続いて、ＣＰＵ６１はステップ８１５に進んで、上記設定されたＩＶＣ時筒内ガス圧力Pgivcと、上記設定されたＩＶＣ時筒内ガス温度Tgivcと、上記(5)式とに基づいて筒内ガスの全質量Mgを求める。

次いで、ＣＰＵ６１はステップ８２０に進み、アクセル開度センサ７５により得られる現時点でのアクセル開度Accp、クランクポジションセンサ７４から取得される現時点でのエンジン回転速度NE、及び図１３に示したテーブル（マップ）MapQfinから指令燃料噴射量Qfin（実際には、燃料噴射期間TAU）を求める。テーブルMapQfinは、アクセル開度Accp及びエンジン回転速度NEと指令燃料噴射量Qfinとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。

次に、ＣＰＵ６１はステップ８２５に進み、指令燃料噴射量Qfin、エンジン回転速度NE、及び図１４に示したテーブルMapCAinjbaseから基本燃料噴射開始時期（クランク角度）CAinjbaseを決定する。テーブルMapCAinjbaseは、指令燃料噴射量Qfin及びエンジン回転速度NEと基本燃料噴射開始時期CAinjbaseとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。

続いて、ＣＰＵ６１はステップ８３０に進んで、指令燃料噴射量Qfin、エンジン回転速度NE、及び図１５に示したテーブルMapPcrbaseから基本燃料噴射圧力Pcrbaseを決定する。テーブルMapPcrbaseは、指令燃料噴射量Qfin及びエンジン回転速度NEと基本燃料噴射圧力Pcrbaseとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。

次に、ＣＰＵ６１はステップ８３５に進み、燃料温度センサ７６から得られる現時点での燃料温度Tcrと、上記(28)式とに基づいて燃料温度Tfを求める。続いて、ＣＰＵ６１はステップ８４０に進んで、上記吸気酸素濃度[O2]inと、[O2]inを引数とする筒内ガス理論空燃比stoichを求めるための関数funcstoichとに基づいて筒内ガス理論空燃比stoichを求める。

次に、ＣＰＵ６１はステップ８４５に進み、現時点でのエンジン回転速度NEと、微小時間Δｔ（例えば、０．１msec）と、NE，Δｔを引数とする微小クランク角度ΔCAを求めるための関数funcΔCAとに基づいて、同微小時間Δｔに相当するクランク角度である微小クランク角度ΔCAを求める。この微小クランク角度ΔCAは、エンジン回転速度NEが現時点（即ち、ＩＶＣ直後）での値である場合における、微小時間Δｔに相当するクランク角度である。

続いて、ＣＰＵ６１はステップ８５０に進み、先のステップ８１５にて求めた筒内ガスの全質量Mgを、先のステップ８２５にて求めた基本燃料噴射開始時期CAinjbaseから得られる燃料噴射開始時筒内容積Vg(CAinjbase)で除することで燃料噴射開始時筒内ガス密度ρg0を求める。

続いて、ＣＰＵ６１はステップ８５５に進み、先のステップ６１０にて求めたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記燃料噴射開始時筒内容積Vg(CAinjbase)と、上記(4)式に相当する式とに基づいて燃料噴射開始時筒内ガス圧力Pg0を求める。

次に、ＣＰＵ６１はステップ８６０に進んで、先のステップ８３０にて求めた基本燃料噴射圧力Pcrbaseから上記燃料噴射開始時筒内ガス圧力Pg0を減じることで有効噴射圧力ΔPを求め、続くステップ８６５にて、上記求めた有効噴射圧力ΔPと、筒内ガス密度ρg0と、テーブルMapθとに基づいて噴霧角θ（図３を参照）を求める。

次いで、ＣＰＵ６１は図９のステップ９０５に進み、空気過剰率前回値λbを初期値「０」に設定し、続くステップ９１０にて混合気形成筒内ガス質量Gの値を初期値「０」に設定するとともに、続くステップ９１５にて燃料消費量積算値sumqr、及び筒内ガス消費量積算値sumgrを共に初期値「０」に設定する。

続いて、ＣＰＵ６１はステップ９２０に進んで、混合気のエンタルピHmixを、上記(27)式に相当する式に従って、初期値（即ち、先のステップ８２０にて求めた指令燃料噴射量Qfinと、燃料の定圧比熱Cfと、先のステップ８３５にて求めた燃料温度Tfの積）に設定する。

次に、ＣＰＵ６１はステップ９２５に進み、混合気の定圧比熱Cmixを初期値である上記燃料の定圧比熱Cfに設定し、続くステップ９３０にて混合気質量Mmixを初期値である指令燃料噴射量Qfinに設定する。

次いで、ＣＰＵ６１はステップ９３５に進んで、酸素濃度[O2]mixを初期値「０」に設定するとともに、混合気内の燃料濃度[Fuel]mixを初期値「１」に設定する。続いて、ＣＰＵ６１はステップ９４０に進み、噴射後経過時間ｔを初期値「０」に設定するとともに、クランク角度CAを初期値である基本燃料噴射開始時期CAinjbaseに設定する。これにより、噴射後経過時間ｔが基本燃料噴射開始時期CAinjbaseからカウントされることになる。

続いて、ＣＰＵ６１はステップ９４５に進んで、混合気温度Tmixを初期値である、ステップ８３５にて求められている燃料温度Tfに設定し、続くステップ９５０にて混合気到達距離Xを初期値「０」に設定する。

そして、ＣＰＵ６１はステップ９５５に進み、フラグＸＢＵ、フラグＸＶＡの値を共に初期値「０」に設定する。ここで、フラグＸＢＵは、その値が「１」のとき混合気のブレークアップが発生していることを示し、その値が「０」のとき混合気のブレークアップが発生していないことを示す。フラグＸＶＡは、その値が「１」のとき混合気内の液滴燃料が蒸発していることを示し、その値が「０」のとき混合気内の液滴燃料が蒸発していないことを示す。このようにして、各種初期値が決定される。

次に、ＣＰＵ６１は図１０のルーチンに進み、（補正後の）空気過剰率λの算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ１００５に進み、噴射後経過時間ｔの値（現時点では、「０」）を微小時間Δｔだけ増大・更新させるとともに、クランク角度CAの値（現時点では、「CAinjbase」）を微小クランク角度ΔCAだけ増大・更新させる。このように、クランク角度CAの値が噴射後経過時間ｔに対応する値に維持されていく。これにより、以降、噴射後経過時間ｔ＝Δｔとなり、クランク角度CA＝CAinjbase＋ΔCAとなる。

続いて、ＣＰＵ６１はステップ１０１０に進んで、先のステップ８１５にて求めた筒内ガスの全質量Mgを、上記ステップ１００５にて更新されたクランク角度CAに対応する筒内容積Vg(CA)で除することで噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス密度ρgを求める。

次に、ＣＰＵ６１はステップ１０１５に進んで、上記ステップ８１０にて求められたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記クランク角度CAに対応する筒内容積Vg(CA)と、上記(32)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス圧力Pgを求める。

次いで、ＣＰＵ６１はステップ１０２０に進み、上記ステップ８１０にて求められたＩＶＣ時筒内ガス温度Tgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記クランク角度CAに対応する筒内容積Vg(CA)と、上記(22)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス温度Tgを求める。

次に、ＣＰＵ６１はステップ１０２５に進んで、上記ステップ８１０にて求められた上記吸気酸素濃度[O2]inと、ステップ１０２０にて求めた筒内ガス温度Tgと、[O2]in，Tgを引数とする筒内ガスの定圧比熱Cgを求めるための関数funcCgと、上記(23)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガスの定圧比熱Cgを求める。

続いて、ＣＰＵ６１はステップ１０３０に進み、先のステップ１０１０にて求めた筒内ガス密度ρgと、先のステップ８６５にて求めた噴霧角θと、先のステップ８６０にて求めた有効噴射圧力ΔPと、先のステップ１００５にて更新した噴射後経過時間ｔと、上記(3)式とに基づいて基本混合気希釈速度(dλ/dt)baseを求める。

次いで、ＣＰＵ６１はステップ１０３５を経由して図１８にフローチャートにより示した希釈係数K1を決定するためのルーチンの処理を開始する。即ち、ＣＰＵ６１はステップ１８００からステップ１８０５に進むと、フラグＸＢＵの値が「０」であるか否かを判定する。

現時点では、先のステップ９５５の処理によりフラグＸＢＵの値は「０」になっている。従って、ＣＰＵ６１はステップ１８０５にて「Ｙｅｓ」と判定してステップ１８１０に進み、先のステップ１０１０にて求めた筒内ガス密度ρgと、ρgを引数とする上記(6)式に対応する関数funcXbreakとに基づいてブレークアップ発生位置Xbreakを求める。

続いて、ＣＰＵ６１はステップ１８１５に進んで、先のステップ１０１０にて求めた筒内ガス密度ρgと、先のステップ８６５にて求めた噴霧角θと、先のステップ８６０にて求めた有効噴射圧力ΔPと、先のステップ１００５にて更新した噴射後経過時間ｔと、上記(8)式とに基づいて混合気移動速度dX/dtを求める。

次に、ＣＰＵ６１はステップ１８２０に進み、上記(7)式に従って、混合気到達距離Xを、その時点での値（現時点では、ステップ９５０の処理により「０」）に、上記求めた混合気移動速度dX/dtに微小時間Δｔを乗じた値「dX/dt・Δｔ」を加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気到達距離Xが求められる。

次いで、ＣＰＵ６１はステップ１８２５に進んで、上記更新された混合気到達距離Xがステップ１８１０にて求めたブレークアップ発生位置Xbreak未満である（即ち、ブレークアップが発生していない）か否かを判定する。現時点は噴射開始直後に対応するから混合気到達距離Xはブレークアップ発生位置Xbreak未満である。従って、ＣＰＵ６１はステップ１８２５にて「Ｙｅｓ」と判定してステップ１８３０に進み、希釈係数K1の値を「１」に設定してステップ１８９５に進む。

一方、（その後、図１８のルーチン（特に、ステップ１８２０）が繰り返し実行されて）混合気到達距離Xがブレークアップ発生位置Xbreak以上となっている場合、ＣＰＵ６１はステップ１８２５にて「Ｎｏ」と判定してステップ１８３５に進み、フラグＸＢＵの値を「０」から「１」に変更し、その後、ステップ１８４０に進んで、希釈係数K1の値を値α（＞１）に設定してステップ１８９５に進む。この場合、以降、フラグＸＢＵの値が「１」になっているから、ＣＰＵ６１はステップ１８０５に進んだとき、「Ｎｏ」と判定してステップ１８４０に直ちに進むようになる。このようにして、希釈係数K1の値は、ブレークアップが発生する前では「１」に設定され、ブレークアップが発生した後は値αに設定される。

ＣＰＵ６１はステップ１８９５を経由して図１０のステップ１０３５に戻ると、次に、ステップ１０４０を経由して図１９にフローチャートにより示した希釈係数K2を決定するためのルーチンの処理を開始する。即ち、ＣＰＵ６１はステップ１９００からステップ１９０５に進むと、フラグＸＶＡの値が「０」であるか否かを判定する。

現時点では、先のステップ９５５の処理によりフラグＸＶＡの値は「０」になっている。従って、ＣＰＵ６１はステップ１９０５にて「Ｙｅｓ」と判定してステップ１９１０に進み、先のステップ１０１５にて求めた筒内ガス圧力Pgと、Pgを引数とする燃料の蒸発温度Tvapを求めるための関数funcTvapとに基づいて燃料の蒸発温度Tvapを求める。

続いて、ＣＰＵ６１はステップ１９１５に進んで、混合気温度Tmixの値（現時点では、ステップ９４５の処理により燃料温度Tfの値と等しい）がステップ１９１０にて求めた燃料の蒸発温度Tvap未満である（即ち、液滴燃料が蒸発していない）か否かを判定する。現時点は噴射開始直後に対応するから混合気温度Tmixは蒸発温度Tvap未満である。従って、ＣＰＵ６１はステップ１９１５にて「Ｙｅｓ」と判定してステップ１９２０に進み、希釈係数K2の値を「１」に設定してステップ１９９５に進む。

一方、（その後、図１９のルーチンが繰り返し実行されて）後述するステップ１１４５にて繰り返し更新されている混合気温度Tmixがステップ１９１０にて求めた蒸発温度Tvap以上となっている場合、ＣＰＵ６１はステップ１９１５にて「Ｎｏ」と判定してステップ１９２５に進み、フラグＸＶＡの値を「０」から「１」に変更し、その後、ステップ１９３０に進んで、混合気温度Tmixと、図４に示したグラフに対応するテーブルとに基づいて希釈係数K2の値（＜１）を決定してステップ１９９５に進む。この場合、以降、フラグＸＶＡの値が「１」になっているから、ＣＰＵ６１はステップ１９０５に進んだとき、「Ｎｏ」と判定してステップ１９３０に直ちに進むようになる。このようにして、希釈係数K2の値は、液滴燃料が蒸発する前では「１」に設定され、液滴燃料が蒸発した後は「１」より小さい値に設定される。

ＣＰＵ６１はステップ１９９５を経由して図１０のステップ１０４０に戻ると、次に、ステップ１０４５を経由して図２０にフローチャートにより示した希釈係数K3を決定するためのルーチンの処理を開始する。即ち、ＣＰＵ６１はステップ２０００からステップ２００５に進むと、フラグＸＶＡの値が「０」であるか否かを判定する。

現時点では、上述のごとく、フラグＸＶＡの値は「０」になっている。従って、ＣＰＵ６１はステップ２００５にて「Ｙｅｓ」と判定してステップ２０１０に進んで、希釈係数K3の値を「１」に設定してステップ２０９５に進む。

一方、（その後、上述したステップ１９２５の処理が実行されて）フラグＸＶＡの値が「０」から「１」に変更されている場合、ＣＰＵ６１はステップ２００５にて「Ｎｏ」と判定してステップ２０１５に進み、ステップ１０１５、後述するステップ１１１５、１１４５にてそれぞれ繰り返し更新されている筒内ガス圧力Pg、混合気質量Mmix、混合気温度Tmixと、及び上記(10)式、(11)式に相当するステップ２０１５内に記載の式とに基づいて半径r0を求める。

続いて、ＣＰＵ６１はステップ２０２０に進み、後述するステップ１１１５、１１３５、１１７５にてそれぞれ繰り返し更新されている混合気質量Mmix、反応熱Hr、混合気の定圧比熱Cmixと、上記(13)式とに基づいて温度上昇量ΔTを求める。次に、ＣＰＵ６１はステップ２０２５に進み、上記混合気温度Tmixと、上記温度上昇量ΔTと、ステップ２０２５内に記載の式とに基づいて値Ａを求め、続くステップ２０３０にて、上記半径r0と、値Ａと、上記(9)式、(12)式に相当するステップ２０３０内に記載の式とに基づいて膨張速度dr/dtを求める。

次いで、ＣＰＵ６１はステップ２０３５に進み、上記求めた膨張速度dr/dtに係数Krを乗じることで乱れ強度Ｂを求め、続くステップ２０４０にて、乱れ強度Ｂの値と、図６に示したグラフに対応するテーブルとに基づいて希釈係数K3の値（＞１）を決定してステップ２０９５に進む。この場合、以降も、フラグＸＶＡの値が「１」になっているから、ＣＰＵ６１はステップ２００５に進んだとき「Ｎｏ」と判定して、ステップ２０１５〜２０４０の処理を繰り返し実行する。このようにして、希釈係数K3の値は、液滴燃料が蒸発する前では「１」に設定され、液滴燃料が蒸発した後は「１」より大きい値に設定される。

ＣＰＵ６１はステップ２０９５を経由して図１０のステップ１０４５に戻ると、次に、ステップ１０５０に進んで、先のステップ１０３０にて求めた基本混合気希釈速度(dλ/dt)baseに上記希釈係数K1,K2,K3を乗じることで補正後混合気希釈速度dλ/dtを求める。

続いて、ＣＰＵ６１はステップ１０５５に進んで、上記(15)式に従って、補正後の空気過剰率λを、その時点での空気過剰率前回値λb（現時点では、ステップ９０５の処理により「０」）に、上記求めた補正後混合気希釈速度dλ/dtに微小時間Δｔを乗じた値「dλ/dt・Δｔ」を加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気の補正後の空気過剰率λが求められる。

次に、ＣＰＵ６１は図１１のステップ１１０５に進み、先のステップ８４０にて求めた筒内ガス理論空燃比stoichと、ステップ１０５５にて求めた補正後の空気過剰率λと、空気過剰率前回値λb（現時点では、ステップ９０５の処理により「０」。次回からは、後述するステップ１１５０にて設定されている値）と、ステップ８２０にて設定された指令燃料噴射量Qfinと、上記(21)式に相当する式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気に新たに取り込まれた筒内ガス質量gを求める。

続いて、ＣＰＵ６１はステップ１１１０に進んで、混合気形成筒内ガス質量Gを、その時点での値（現時点では、ステップ９１０の処理により「０」）に上記新たに取り込まれた筒内ガス質量gを加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気形成筒内ガス質量Gが求められる。

次に、ＣＰＵ６１はステップ１１１５に進み、混合気質量Mmixを、その時点での値（現時点では、ステップ９３０の処理により指令燃料噴射量Qfinと等しい）に上記新たに取り込まれた筒内ガス質量gを加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気質量Mmixが求められる。

続いて、ＣＰＵ６１はステップ１１２０に進んで、化学反応前の混合気のエンタルピHpreを、その時点での混合気のエンタルピHmix（現時点では、ステップ９２０の処理により値「Qfin・Cf・Tf」）に、上記(20)式に従って求められる「上記新たに取り込まれた筒内ガスの熱エネルギーΔHg＝g・Cg・Tg」を加えた値に設定する。

次に、ＣＰＵ６１はステップ１１２５に進み、上記求めた化学反応前の混合気のエンタルピHpreを、上記ステップ１１１５にて求めた混合気質量Mmixにその時点での混合気の定圧比熱Cmix（現時点では、ステップ９２５の処理により燃料の定圧比熱Cf。次回からは、後述するステップ１１７５にて設定されている値）を乗じた値で除することにより上記(26)式に相当する式に従って化学反応前の混合気の温度Tpreを求める。

次いで、ＣＰＵ６１はステップ１１３０に進んで、その時点での混合気内の酸素濃度[O2]mix（現時点では、ステップ９３５の処理により「０」。次回からは、後述するステップ１１７０にて設定されている値）と、燃料濃度[Fuel]mix（現時点では、ステップ９３５の処理により「１」。次回からは、後述するステップ１１６５にて設定されている値）と、上記求めた化学反応前の混合気の温度Tpreと、上記(25)式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrを求める。

続いて、ＣＰＵ６１はステップ１１３５に進み、上記求めた消費燃料量qrと、上記(24)式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による反応熱Hrを求め、続くステップ１１４０にて混合気のエンタルピHmixを、上記求めた化学反応前の混合気のエンタルピHpreに上記求めた反応熱Hrを加えた値に設定・更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気のエンタルピHmixが求められる。

そして、ＣＰＵ６１はステップ１１４５に進んで、上記ステップ１１４０にて求めた混合気のエンタルピHmixと、上記ステップ１１１５にて求めた混合気質量Mmixと、その時点での混合気の定圧比熱Cmix（現時点では、ステップ９２５の処理により燃料の定圧比熱Cf。次回からは、後述するステップ１１７５にて設定されている値）と、上記(17)式とに基づいて混合気温度Tmixを算出する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj＋ΔCA）における混合気温度Tmixが求められる。

次に、ＣＰＵ６１はステップ１１５０に進んで、次回の計算の準備のため、空気過剰率前回値λbを上記ステップ１０５５にて求めた空気過剰率λの値に設定する。この値は、以降、上述したステップ１０５５にて使用される。このようにして、噴射後経過時間ｔにおける混合気温度Tmixが算出される。

次に、ＣＰＵ６１はステップ１１５５に進み、上記ステップ１１３０にて求めた微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrと、ステップ８４０にて求めた筒内ガス理論空燃比stoichと、上記(30)式とに基づいて、微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による筒内ガス消費量grを求める。

続いて、ＣＰＵ６１はステップ１１６０に進んで、燃料消費量積算値sumqrを、その時点での値（現時点では、ステップ９１５の処理により「０」）にステップ１１３０にて求めた上記燃料消費量qrを加えた値に設定・更新し、筒内ガス消費量積算値sumgrを、その時点での値（現時点では、ステップ９１５の処理により「０」）にステップ１１５５にて求めた上記筒内ガス消費量grを加えた値に設定・更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における燃料消費量積算値sumqr、及び筒内ガス消費量積算値sumgrが求められる。

次いで、ＣＰＵ６１はステップ１１６５に進み、ステップ８２０にて求めた指令燃料噴射量Qfinと、上記求めた燃料消費量積算値sumqrと、ステップ１１１５にて求めた混合気質量Mmixと、上記(29)式に相当する式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における混合気内の燃料濃度[Fuel]mixを求める。

次に、ＣＰＵ６１はステップ１１７０に進んで、ステップ１１１０にて求めた混合気形成筒内ガス質量Gと、上記求めた筒内ガス消費量積算値sumgrと、ステップ８１０にて設定された吸気酸素濃度[O2]inと、ステップ１１１５にて求めた混合気質量Mmixと、上記(31)式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における混合気内の酸素濃度[O2]mixを求める。

続いて、ＣＰＵ６１はステップ１１７５に進み、上記ステップ１１７０にて求めた混合気内の酸素濃度[O2]mixと、ステップ１１４５にて求めた混合気温度Tmixと、上記(19)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における混合気の定圧比熱Cmixを求める。この値は、以降、ステップ１１２５、１１４５にて使用される。

次いで、ＣＰＵ６１は図１２のステップ１２０５に進み、クランク角度CAが目標着火時期CArefに一致したか否かを判定する。現時点では、上述のごとく、クランク角度CAは、基本燃料噴射開始時期CAinjbaseに微小クランク角度ΔCAを加えた値であるから目標着火時期CArefに達していない。

従って、現時点では、ＣＰＵ６１はステップ１２０５にて「Ｎｏ」と判定して図１０のステップ１００５に戻り、噴射後経過時間ｔ（現時点では、「1・Δｔ」）を微小時間Δｔだけ増大・更新させるとともに、クランク角度CA（現時点では、「CAinjbase＋ΔCA」）を微小クランク角度ΔCAだけ増大・更新させた後、上述した図１０のステップ１０１０〜図１２のステップ１２０５の処理を再び実行する。

これにより、噴射後経過時間ｔ＝２・Δｔ（クランク角度CA＝CAinjbase＋２・ΔCA）における、混合気の空気過剰率λ（ステップ１０５５を参照。）、混合気温度Tmix（ステップ１１４５を参照。）、及び各種濃度（[Fuel]mix，[O2]mix。ステップ１１６５、１１７０を参照。）が算出される。

そして、図１２のステップ１２０５の判定において「Ｎｏ」と判定される毎に、図１０のステップ１００５〜図１２のステップ１２０５の処理が繰り返し実行されていく。これにより、図１２のステップ１２０５の判定において「Ｎｏ」と判定される限りにおいて、混合気の空気過剰率λ、混合気温度Tmix、及び各種濃度（[Fuel]mix，[O2]mix）が微小時間Δｔ毎に（即ち、CAinjbaseから微小クランク角度ΔCA毎に）更新されていく。

そして、クランク角度CAが目標着火時期CArefに達した場合、ＣＰＵ６１は図１２のステップ１２０５に進んだとき「Ｙｅｓ」と判定してステップ１２１０以降に進む。

ＣＰＵ６１はステップ１２１０に進むと、混合気温度偏差ΔTmixを、上記目標混合気温度Tmixrefからステップ１１４５にて求められている現時点での混合気温度Tmix（上記制御用混合気温度）を減じた値に設定する。

続いて、ＣＰＵ６１はステップ１２１５に進んで、噴射時期補正値Δθを、上記混合気温度偏差ΔTmixと図１６に示したテーブルMapΔθとから決定する。テーブルMapΔθは、混合気温度偏差ΔTmixと噴射時期補正値Δθとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。

その後、ＣＰＵ６１はステップ１１２０に進んで、噴射圧力補正値ΔPcrを、上記混合気温度偏差ΔTmixと図１７に示したテーブルMapΔPcrとから決定する。テーブルMapΔPcrは、混合気温度偏差ΔTmixと噴射圧力補正値ΔPcrとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。

次いで、ＣＰＵ６１はステップ１２２５に進み、基本燃料噴射開始時期CAinjbaseを噴射時期補正値Δθで補正して最終燃料噴射開始時期CAinjfinを決定する。これにより、混合気温度偏差ΔTmixに応じて噴射時期が補正されることになる。この場合、図１６から明らかなように、混合気温度偏差ΔTmixが正の大きい値になるほど噴射時期補正値Δθが正の大きな値となって最終燃料噴射開始時期CAinjfinが進角側となり、混合気温度偏差ΔTmixが負の大きい値（絶対値が大きい値）になるほど噴射時期補正値Δθは負の大きな値となって最終燃料噴射開始時期CAinjfinが遅角側に移行される。

続いて、ＣＰＵ６１はステップ１２３０に進み、基本燃料噴射圧力Pcrbaseを噴射圧力補正値ΔPcrで補正して指令最終燃料噴射圧力Pcrfinを決定する。これにより、混合気温度偏差ΔTmixに応じて噴射圧力が補正されることになる。この場合、図１７から明らかなように、混合気温度偏差ΔTmixが正の大きい値になるほど噴射圧力補正値ΔPcrが正の大きな値となって指令最終燃料噴射圧力Pcrfinが高圧側となり、同混合気温度偏差ΔTmixが負の大きい値（絶対値が大きい値）になるほど噴射圧力補正値ΔPcrは負の大きな値となって指令最終燃料噴射圧力Pcrfinが低圧側に移行される。

そして、ＣＰＵ６１はステップ１２３５に進むと、実際の燃料噴射圧力が上記設定された指令最終燃料噴射圧力Pcrfinになるように燃料噴射用ポンプ２２（の駆動回路）に対して制御指示を行い、ステップ８９５に進んで図８〜図１２の一連の本ルーチンを一旦終了する。以降、ＣＰＵ６１は、次のＩＶＣが到来するまでの間、ステップ８０５に進む毎に「Ｎｏ」と判定し続ける。

この結果、本ルーチンの実行により、ＩＶＣが到来する毎に、燃料噴射形態（噴射量、噴射圧力、噴射時期）が決定されるとともに空気過剰率λ、混合気温度Tmix等が直ちに推定され、係る推定結果に基づいて噴射時期、噴射圧力が補正されていく。

また、ＣＰＵ６１は、図２１にフローチャートにより示した燃料噴射制御を行うためのルーチンを所定時間の経過毎に、気筒毎に、繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１はステップ２１００から処理を開始し、ステップ２１０５に進んで実クランク角度CAactが先のステップ１２２５にて決定されている最終燃料噴射開始時期CAinjfinに一致したか否かを判定し、「Ｎｏ」と判定する場合、ステップ２１９５に直ちに進んで本ルーチンを一旦終了する。

いま、実クランク角度CAactが最終燃料噴射開始時期CAinjfinに一致したものとすると、ＣＰＵ６１はステップ２１１０に進んで、対応する燃料噴射弁２１に対してステップ８２０にて決定されている指令燃料噴射量Qfinの燃料の噴射指示（具体的には、燃料噴射期間TAUに亘る開弁指示）を行い、ステップ２１９５に進んで本ルーチンを一旦終了する。これにより、指令燃料噴射量Qfinの燃料が上記指令最終燃料噴射圧力Pcrfinをもって噴射される。

以上、説明したように、本発明による混合気濃度推定装置の実施形態を含んだ混合気温度推定装置によれば、運動量保存則を利用して得られる式（上記(3)式を参照）から算出される基本混合気希釈速度(dλ/dt)baseの値に、ブレークアップ、蒸発、及び熱膨張にそれぞれ対応する希釈係数K1,K2,K3を乗じた値である補正後混合気希釈速度dλ/dt(＝(dλ/dt)base・K1・K2・K3)を、噴射後経過時間ｔ＝０から微小時間Δｔ（例えば、０．１msec）毎に求める。

希釈係数K1,K2,K3は、噴射後経過時間ｔ＝０では全て「１」に設定される。混合気のブレークアップが発生すると（図７のｔ＝ｔ１を参照）、それ以降、K1のみが値α（＞１、一定値）に変更される。混合気内の液滴燃料が蒸発すると（図７のｔ＝ｔ２）、それ以降、希釈係数K2は「１」未満の値であって混合気温度Tmixが高いほどより小さい値（図４を参照）に設定されるとともに、希釈係数K3は「１」より大きい値であって膨張速度dr/dt（乱れ強度Ｂ）が大きいほどより大きい値（図６を参照）に設定される。

これにより、混合気が燃焼室内を移動していく過程において燃料と筒内ガスとが混ざり合う程度を変化させる現象となる混合気のブレークアップ、混合気内の液滴燃料の蒸発、及び気体としての混合気の熱膨張が考慮されて補正後混合気希釈速度dλ/dtが噴射後経過時間ｔ＝０から微小時間Δｔ毎に算出される。

そして、それぞれの補正後混合気希釈速度dλ/dtの値を時間で順に積分（積算）していくことで混合気の空気過剰率λが噴射後経過時間ｔ＝０から微小時間Δｔ毎に求められる。従って、混合気のブレークアップ、混合気内の液滴燃料の蒸発、及び気体としての混合気の熱膨張が考慮されて混合気の空気過剰率λが精度良く算出され得る。

更には、このように精度良く算出された空気過剰率λを利用して噴射後経過時間ｔにおける混合気温度Tmixが算出されるから、混合気温度Tmixも精度良く算出され得る。

本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記実施形態では、基本混合気希釈速度(dλ/dt)baseの値に、ブレークアップ、蒸発、及び熱膨張にそれぞれ対応する希釈係数K1,K2,K3を全て乗じることで基本混合気希釈速度(dλ/dt)baseを補正しているが、希釈係数K1,K2,K3のうちの少なくとも１つを基本混合気希釈速度(dλ/dt)baseに乗じることで基本混合気希釈速度(dλ/dt)baseを補正してもよい。

また、上記実施形態においては、混合気内の液滴燃料が蒸発した時点以降、希釈係数K2を「１」未満の値であって混合気温度Tmixに応じた値に設定しているが、希釈係数K2を「１」未満の一定値に設定してもよい。同様に、混合気内の液滴燃料が蒸発した時点以降、希釈係数K3を「１」より大きい値であって乱れ強度Ｂに応じた値に設定しているが、希釈係数K3を「１」より大きい一定値に設定してもよい。

また、上記実施形態においては、ブレークアップ、蒸発、及び熱膨張にそれぞれ対応する希釈係数K1,K2,K3を求め、基本混合気希釈速度(dλ/dt)baseの値に希釈係数K1,K2,K3を乗じることで基本混合気希釈速度(dλ/dt)baseを補正しているが、ブレークアップ、蒸発、及び熱膨張に対応する補正値をそれぞれ求め、基本混合気希釈速度(dλ/dt)baseの値に各補正値をそれぞれ加減算することで基本混合気希釈速度(dλ/dt)baseを補正してもよい。

加えて、上記実施形態においては、混合気内の液滴燃料が蒸発した時点以降において希釈係数K3を求めるに際し、混合気の膨張速度dr/dt（従って、乱れ強度Ｂ）を球形状の混合気が熱膨張すると仮定して算出しているが、楕円体、多面体、垂体等の形状を呈した混合気が熱膨張すると仮定して混合気の膨張速度（従って、乱れ強度Ｂ）を算出してもよい。

本発明の実施形態に係る内燃機関の混合気濃度推定装置を含んだ混合気温度推定装置を４気筒内燃機関（ディーゼル機関）に適用したシステム全体の概略構成図である。 或る一つの気筒のシリンダ内(筒内)に吸気マニホールドからガスが吸入され、筒内に吸入されたガスが排気マニホールドへ排出される様子を模式的に示した図である。 （ａ）は、燃料が一時に噴射された時点での同燃料の様子を模式的に示した図であり、（ｂ）は、筒内ガスと混ざり合いながら混合気となって円錐状に拡散していく燃料の様子を模式的に示した図である。 混合気希釈速度を混合気内の液滴燃料の蒸発を考慮して補正するために使用される希釈係数の値と、混合気温度との関係を示したグラフである。 球形状の混合気が熱膨張する際の様子を示した模式図である。 混合気希釈速度を混合気の熱膨張を考慮して補正するために使用される希釈係数の値と、乱れ強度（膨張速度）との関係を示したグラフである。 図１に示した混合気濃度推定装置により求められた噴射後経過時間に対する空気過剰率の変化の一例を示したグラフである。 図１に示したＣＰＵが実行する空気過剰率、及び混合気温度の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する空気過剰率、及び混合気温度の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する空気過剰率、及び混合気温度の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する空気過剰率、及び混合気温度の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する空気過剰率、及び混合気温度の算出を行うためのルーチンの第５番目の部分を示したフローチャートである。 図１に示したＣＰＵが図８に示したルーチンを実行する際に参照する指令燃料噴射量を決定するためのテーブルである。 図１に示したＣＰＵが図８に示したルーチンを実行する際に参照する基本燃料燃料噴射時期を決定するためのテーブルである。 図１に示したＣＰＵが図８に示したルーチンを実行する際に参照する基本燃料噴射圧力を決定するためのテーブルである。 図１に示したＣＰＵが図１２に示したルーチンを実行する際に参照する噴射時期補正値を決定するためのテーブルである。 図１に示したＣＰＵが図１２に示したルーチンを実行する際に参照する噴射圧力補正値を決定するためのテーブルである。 図１に示したＣＰＵが実行する混合気のブレークアップに対応する希釈係数を決定するためのルーチンを示したフローチャートである。 図１に示したＣＰＵが実行する混合気内の液滴燃料の蒸発に対応する希釈係数を決定するためのルーチンを示したフローチャートである。 図１に示したＣＰＵが実行する混合気の熱膨張に対応する希釈係数を決定するためのルーチンを示したフローチャートである。 図１に示したＣＰＵが実行する燃料噴射制御を行うためのルーチンを示したフローチャートである。

符号の説明

２１…燃料噴射弁、６０…電気制御装置、６１…ＣＰＵ、７２…吸気温センサ、７３…吸気管圧力センサ、７４…クランクポジションセンサ、７６…燃料温度センサ、７７…吸気酸素濃度センサ

Claims (7)

1. 内燃機関の燃焼室内に噴射された後において前記燃焼室内に吸入されている筒内ガス中にて時々刻々と拡散していく燃料と、前記燃料と混ざり合って混合気を形成する前記筒内ガスの一部分であってその量が前記燃料の拡散の進行に伴って時々刻々と増大していく混合気形成筒内ガスと、からなる混合気の濃度に相当する値である混合気濃度相当値の基本値を、前記燃料の噴射後において発生する前記混合気のブレークアップを考慮することなく且つ少なくとも運動量保存則を利用して、前記燃料の噴射後において逐次算出する混合気濃度算出手段と、
   前記混合気のブレークアップが発生したか否かを前記燃料の噴射後において逐次判定するブレークアップ判定手段と、
   前記ブレークアップが発生したと判定された後において、前記混合気濃度算出手段により算出される混合気濃度相当値の基本値を前記ブレークアップに基づく分だけ前記混合気の濃度の減少方向に補正して最終的な混合気濃度相当値を逐次算出する混合気濃度補正手段と、
   を備えた内燃機関の混合気濃度推定装置。
2. 請求項１に記載の内燃機関の混合気濃度推定装置において、
   前記ブレークアップ判定手段は、
   前記燃焼室内を移動していく前記混合気の位置を前記燃料の噴射後において逐次算出する混合気位置算出手段と、
   前記燃焼室内において前記混合気のブレークアップの発生位置を算出するブレークアップ位置算出手段と、
   を備え、
   前記算出された混合気の位置が前記算出されたブレークアップの発生位置に達したとき、前記ブレークアップが発生したと判定するように構成された内燃機関の混合気濃度推定装置。
3. 請求項１又は請求項２に記載の内燃機関の混合気濃度推定装置であって、
   前記混合気内の燃料が蒸発したか否かを前記燃料の噴射後において逐次判定する蒸発判定手段を備え、
   前記混合気濃度算出手段は、
   前記燃料の噴射後において発生する前記混合気の燃料の蒸発を考慮することなく前記混合気濃度相当値の基本値を逐次算出するように構成されていて、
   前記混合気濃度補正手段は、
   前記燃料が蒸発したと判定された後において、前記混合気濃度算出手段により算出される混合気濃度相当値の基本値を前記燃料の蒸発に基づく分だけ前記混合気の濃度の増大方向に補正して最終的な混合気濃度相当値を逐次算出するように構成された内燃機関の混合気濃度推定装置。
4. 請求項３に記載の内燃機関の混合気濃度推定装置において、
   前記蒸発判定手段は、
   前記混合気の温度を前記燃料の噴射後において逐次算出する混合気温度算出手段を備え、
   前記算出された混合気の温度が前記燃料の蒸発温度に達したとき、前記燃料が蒸発したと判定するように構成された内燃機関の混合気濃度推定装置。
5. 請求項４に記載の内燃機関の混合気濃度推定装置において、
   前記混合気濃度補正手段は、
   前記混合気の温度が高いほど、前記燃料の蒸発に基づく前記混合気の濃度の増大方向への補正の量を大きくするように構成された内燃機関の混合気濃度推定装置。
6. 請求項３乃至請求項５の何れか一項に記載の内燃機関の混合気濃度推定装置において、
   前記混合気濃度算出手段は、
   前記燃料の噴射後において前記混合気内の反応により発生する熱による前記混合気の膨張を考慮することなく前記混合気濃度相当値の基本値を逐次算出するように構成されていて、
   前記混合気濃度補正手段は、
   前記燃料が蒸発したと判定された後において、前記混合気濃度算出手段により算出される混合気濃度相当値の基本値を前記混合気の膨張に基づく分だけ前記混合気の濃度の減少方向に補正して最終的な混合気濃度相当値を逐次算出するように構成された内燃機関の混合気濃度推定装置。
7. 請求項６に記載の内燃機関の混合気濃度推定装置において、
   前記混合気濃度補正手段は、
   前記混合気の膨張の程度が大きいほど、前記混合気の膨張に基づく前記混合気の濃度の減少方向への補正の量を大きくするように構成された内燃機関の混合気濃度推定装置。

Images