JP4424242B2 - 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置 - Google Patents

内燃機関の混合気状態推定装置、及びエミッション発生量推定装置 Download PDF

Info

Publication number
JP4424242B2
JP4424242B2 JP2005098148A JP2005098148A JP4424242B2 JP 4424242 B2 JP4424242 B2 JP 4424242B2 JP 2005098148 A JP2005098148 A JP 2005098148A JP 2005098148 A JP2005098148 A JP 2005098148A JP 4424242 B2 JP4424242 B2 JP 4424242B2
Authority
JP
Japan
Prior art keywords
fuel
mixture
injection
air
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2005098148A
Other languages
English (en)
Other versions
JP2006274991A (ja
Inventor
知美 大西
茂樹 中山
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Toyota Motor Corp
Original Assignee
Toyota Motor Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Toyota Motor Corp filed Critical Toyota Motor Corp
Priority to JP2005098148A priority Critical patent/JP4424242B2/ja
Priority to PCT/JP2006/306769 priority patent/WO2006106893A1/ja
Priority to CN2006800106974A priority patent/CN101151453B/zh
Priority to US11/886,595 priority patent/US7848871B2/en
Priority to EP06730717.3A priority patent/EP1865178B1/en
Publication of JP2006274991A publication Critical patent/JP2006274991A/ja
Application granted granted Critical
Publication of JP4424242B2 publication Critical patent/JP4424242B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/30—Controlling fuel injection
    • F02D41/38—Controlling fuel injection of the high pressure type
    • F02D41/40—Controlling fuel injection of the high pressure type with means for controlling injection timing or duration
    • F02D41/402—Multiple injections
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D35/00—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for
    • F02D35/02—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions
    • F02D35/023—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions by determining the cylinder pressure
    • F02D35/024—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions by determining the cylinder pressure using an estimation
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D35/00—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for
    • F02D35/02—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions
    • F02D35/025—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions by determining temperatures inside the cylinder, e.g. combustion temperatures
    • F02D35/026—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions by determining temperatures inside the cylinder, e.g. combustion temperatures using an estimation
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/02—Circuit arrangements for generating control signals
    • F02D41/14—Introducing closed-loop corrections
    • F02D41/1438—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor
    • F02D41/1444—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases
    • F02D41/1454—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being an oxygen content or concentration or the air-fuel ratio
    • F02D41/1458—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being an oxygen content or concentration or the air-fuel ratio with determination means using an estimation
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/02—Circuit arrangements for generating control signals
    • F02D41/14—Introducing closed-loop corrections
    • F02D41/1438—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor
    • F02D41/1444—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases
    • F02D41/146—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being an NOx content or concentration
    • F02D41/1461—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being an NOx content or concentration of the exhaust gases emitted by the engine
    • F02D41/1462—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being an NOx content or concentration of the exhaust gases emitted by the engine with determination means using an estimation
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/02—Circuit arrangements for generating control signals
    • F02D41/14—Introducing closed-loop corrections
    • F02D41/1438—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor
    • F02D41/1444—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases
    • F02D41/1466—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being a soot concentration or content
    • F02D41/1467—Introducing closed-loop corrections using means for determining characteristics of the combustion gases; Sensors therefor characterised by the characteristics of the combustion gases the characteristics being a soot concentration or content with determination means using an estimation
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/02—Circuit arrangements for generating control signals
    • F02D41/14—Introducing closed-loop corrections
    • F02D41/1401—Introducing closed-loop corrections characterised by the control or regulation method
    • F02D2041/1433—Introducing closed-loop corrections characterised by the control or regulation method using a model or simulation of the system
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D2250/00—Engine control related to specific problems or objectives
    • F02D2250/31—Control of the fuel pressure
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D35/00—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for
    • F02D35/02—Controlling engines, dependent on conditions exterior or interior to engines, not otherwise provided for on interior conditions
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02M—SUPPLYING COMBUSTION ENGINES IN GENERAL WITH COMBUSTIBLE MIXTURES OR CONSTITUENTS THEREOF
    • F02M26/00—Engine-pertinent apparatus for adding exhaust gases to combustion-air, main fuel or fuel-air mixture, e.g. by exhaust gas recirculation [EGR] systems
    • F02M26/02—EGR systems specially adapted for supercharged engines
    • F02M26/04—EGR systems specially adapted for supercharged engines with a single turbocharger
    • F02M26/05—High pressure loops, i.e. wherein recirculated exhaust gas is taken out from the exhaust system upstream of the turbine and reintroduced into the intake system downstream of the compressor
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02M—SUPPLYING COMBUSTION ENGINES IN GENERAL WITH COMBUSTIBLE MIXTURES OR CONSTITUENTS THEREOF
    • F02M26/00—Engine-pertinent apparatus for adding exhaust gases to combustion-air, main fuel or fuel-air mixture, e.g. by exhaust gas recirculation [EGR] systems
    • F02M26/13—Arrangement or layout of EGR passages, e.g. in relation to specific engine parts or for incorporation of accessories
    • F02M26/22—Arrangement or layout of EGR passages, e.g. in relation to specific engine parts or for incorporation of accessories with coolers in the recirculation passage
    • F02M26/23—Layout, e.g. schematics
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Mechanical Engineering (AREA)
  • General Engineering & Computer Science (AREA)
  • Combined Controls Of Internal Combustion Engines (AREA)

Description

本発明は、内燃機関の燃焼室内に噴射された燃料が同燃焼室内に吸入されている筒内ガスと混ざり合うことで形成される混合気の状態（温度等）を推定する内燃機関の混合気状態推定装置、及び同混合気の燃焼に起因して発生する有害物質であるエミッションの発生量を推定するエミッション発生量推定装置に関する。
内燃機関（特に、ディーゼル機関）の燃焼室内で発生するスート（Soot、煤煙）、ＮＯｘ等のエミッションの発生量は、燃焼室内にて進行していく混合気の（特に、着火後の混合気）の状態（特に、温度、燃料濃度等）に強い相関を有している。従って、エミッションの発生量を精度良く制御して同エミッションの発生量（従って、排出量）を少なくするためには、混合気の状態を精度良く推定する必要がある。
混合気の状態を推定する技術としては種々のものが知られている。例えば、下記特許文献１に記載のディーゼル機関の燃料噴射装置は、機関冷却水温、吸気温度、吸気圧力等のような混合気温度に影響を与える各種運転状態量に基づいて混合気温度を推定するようになっている。
特開２００１−２５４６４５号公報
ところで、燃焼室内にて進行していく混合気の状態（温度等）は、ピストンの運動に応じて時々刻々と変化する筒内ガスの温度等に依存して順次決定されていく。また、混合気の状態は、燃料に対して筒内ガスが混ざり合っていく程度（従って、例えば、混合気の空気過剰率、燃料濃度等）にも依存する。この混合気の空気過剰率は、同混合気に含まれる燃料の噴射時点での噴射圧力、及び同噴射時点での筒内ガスの密度に依存して決定される値である。即ち、混合気の状態は、同混合気に含まれる燃料の噴射時点での噴射圧力、及び同噴射時点での筒内ガスの密度にも依存して順次決定されていく。
他方、燃料噴射後の経過時間に対する筒内ガスの温度等は、その燃料の噴射時点によって異なる。また、燃料噴射時点での燃料の噴射圧力、及び同噴射時点での筒内ガスの密度も、その燃料の噴射時点によって異なる。従って、燃焼室内にて進行していく混合気の状態は、同混合気に含まれる燃料の噴射時点によって異なることになる。
加えて、燃焼室内に或る量の燃料が噴射される場合、その噴射量に応じて設定された噴射期間だけ連続して燃料が噴射される。従って、燃焼室内にて進行していく混合気を微視的に見れば、同混合気の先頭部分により近い部分に含まれる燃料ほどより早期に噴射されたことになる。換言すれば、混合気が占める領域内の位置（具体的には、混合気先頭位置からの噴孔方向への距離）によって、その位置に対応する燃料の部分の噴射時点が異なる。
以上のことから、燃焼室内にて進行していく混合気を微視的に見た場合、（燃料噴射後の経過時間に対する）同混合気の状態（温度等）は、同混合気が占める領域内の位置に応じて異なる。従って、混合気の燃焼に起因して発生するエミッションの発生度合も同混合気が占める領域内の位置に応じて異なる。換言すれば、燃焼室内にて進行していく混合気の状態が不均一となるから、エミッションの発生度合も不均一になり、その結果、燃焼室内で発生するエミッションの総発生量は、係る不均一性に大きく影響され得る。
しかしながら、上記従来の装置において推定される混合気温度（即ち、混合気の状態）は、上記噴射期間を考慮することなく燃料が一時に（噴射開始時点において瞬時に）噴射されたものとして算出された値である。即ち、この混合気温度は、上述した混合気の不均一性を考慮することなく求められる値である。従って、上記従来の装置により推定される混合気温度に基づいてエミッションの総発生量を推定すると、同推定されるエミッションの総排出量に誤差が生じ、この結果、エミッションの発生量（従って、排出量）を精度良く制御することができないという問題がある。
本発明は、かかる問題に対処するためになされたものであって、その目的は、燃焼室内にて進行していく混合気の不均一性を考慮して混合気の状態を精度良く推定できる内燃機関の混合気状態推定装置、及び、係る混合気の不均一性に基づくエミッション発生度合の不均一性を考慮して燃焼室内でのエミッションの総発生量を精度良く推定できるエミッション発生量推定装置を提供することにある。
本発明に係る混合気状態推定装置の特徴は、内燃機関の燃焼室内にて所定の噴射開始時点から所定の噴射期間だけ連続して噴射される燃料を複数の部分に区分する噴射燃料区分手段と、前記区分された噴射燃料の各部分が前記所定の噴射開始時点からの時間経過に従って独立して順次噴射されていくとの仮定のもと、前記噴射燃料の各部分が前記燃焼室内に吸入されている筒内ガスと混ざり合って形成されていくそれぞれの混合気の状態を個別に推定する混合気状態推定手段とを備えたことにある。ここで、混合気の状態とは、例えば、同混合気の温度、同混合気内の燃料濃度、同混合気内の酸素濃度、同混合気内の窒素濃度等である。
また、上述したように、混合気の状態は、燃料に対して筒内ガスが混ざり合っていく程度（例えば、混合気の空気過剰率、燃料濃度等）に依存して決定されていくから、前記混合気状態推定手段は、前記区分された噴射燃料の各部分に対して前記筒内ガスが混ざり合っていく程度を表す値を個別に取得し、対応する部分に対して前記取得された前記筒内ガスが混ざり合っていく程度を表す値を使用して各部分に基づいて形成されていくそれぞれの混合気の状態を個別に推定するように構成されることが好適である。ここで、前記筒内ガスが混ざり合っていく程度を表す値とは、例えば、混合気の空気過剰率、混合気内の燃料濃度等である。
これによれば、上記区分された噴射燃料の各部分が所定の噴射開始時点からの時間経過に従って独立して順次噴射されていくとの仮定のもと、噴射燃料の各部分に基づいて形成されていくそれぞれの混合気の状態が、対応する部分の噴射時点で同対応する部分の燃料が噴射されたものとして個別に推定され得る。従って、噴射燃料の各部分に基づいて形成されていくそれぞれの混合気間の不均一性が考慮されて、それぞれの混合気の状態が、個別、且つ精度良く推定され得る。
この場合、前記噴射燃料区分手段は、前記所定の噴射期間を所定の規則に従って複数の期間に区分し、前記噴射燃料を同区分された各期間に対応して順次噴射されていくそれぞれの部分に区分するように構成されると好適である。
ここにおいて、上記区分された複数の（区分数個の）期間は、例えば、前記噴射期間を区分数で除することで得られる期間と等しい値に設定されてもよい。この場合、上記複数の期間は均一になる。この場合、上記噴射期間に亘る噴射圧力が一定に維持され得ないこと等に起因して、同各期間に対応する噴射燃料の部分の燃料量は不均一になり得る。また、上記複数の期間は、各期間に対応する噴射燃料の部分の燃料量が等しくなるように設定されてもよい。この場合、上記複数の期間は不均一になり得る。
また、上記区分数は、上記噴射期間に応じて変更してもよい。即ち、例えば、上記区分数は、上記噴射期間が長くなるほど（即ち、総燃料噴射量が大きくなるほど）大きくしてもよい。このようにして、上記構成によれば、上記噴射期間を時間に関して簡易な規則に従って複数の期間に区分するのみで噴射燃料における複数の部分を簡易に得ることができる。
また、上記本発明に係る混合気状態推定装置においては、前記混合気の着火時期を表す値を取得する着火時期取得手段を更に備え、前記噴射燃料区分手段は、前記取得された着火時期を表す値に基づいて前記噴射燃料を前記混合気の着火前に噴射される部分と同混合気の着火後に噴射される部分とに区分するように構成されている。ここにおいて、混合気の着火時期を表す値とは、例えば、混合気の着火時期そのもの、燃料噴射開始時点から混合気の着火までに要する時間（以下、「着火遅れ時間」と称呼する。）等である。
一般に、着火前に噴射された燃料に基づいて形成されていく混合気については、予混合的な燃焼が主として発生し、着火後に噴射された燃料に基づいて形成されていく混合気については、拡散的な燃焼が主として発生することが知られている。他方、予混合的燃焼と拡散的燃焼とでは反応速度（従って、燃焼に起因する発熱量））が異なる（詳細は後述する。）。
換言すれば、着火前に噴射された燃料に基づいて形成されていく混合気と、着火後に噴射された燃料に基づいて形成されていく混合気とでは、燃料噴射後の経過時間に対する温度の増加量が異なる。
上記構成は、係る観点に基づくものである。即ち、上記構成によれば、前記噴射燃料が前記混合気の着火前に噴射される部分と同混合気の着火後に噴射される部分とに区分されるから、着火前に噴射された燃料に基づいて形成されていく混合気の状態（温度等）と、着火後に噴射された燃料に基づいて形成されていく混合気の状態（温度等）とを個別に推定していくことができる。
従って、それぞれの混合気の状態推定において最適な反応速度（従って、発熱量）を設定することができ、この結果、それぞれの混合気の状態（特に、温度）がより一層精度良く推定され得る。
また、上記本発明に関連する混合気状態推定装置においては、前記区分された噴射燃料における第１番目に噴射される部分に対して前記筒内ガスが混ざり合っていく程度を表す値（例えば、上記空気過剰率等）である第１番目部分混合指標値を同第１番目に噴射される部分の前記燃焼室内での進行に対応させて取得していく混合指標値手段を更に備え、前記混合気状態推定手段は、前記第１番目に噴射される部分に基づいて形成されていく混合気の状態を前記取得されている第１番目部分混合指標値を使用して推定するとともに、前記区分された噴射燃料における第２番目以降に噴射される部分に対して前記筒内ガスが混ざり合っていく程度を表す値であって同第２番目以降に噴射される部分の前記燃焼室内での進行に対応する値が同進行に対応する前記取得されている前記第１番目部分混合指標値に基づいて決定されるとの仮定のもと、同決定された第２番目以降に噴射される部分に対して前記筒内ガスが混ざり合っていく程度を表す値を使用して、前記第２番目以降に噴射される部分に基づいて形成されていく混合気の状態を推定するように構成されてもよい。
この場合、「前記区分された噴射燃料における第２番目以降に噴射される部分に対して前記筒内ガスが混ざり合っていく程度を表す値であって同第２番目以降に噴射される部分の前記燃焼室内での進行（例えば、噴射からの経過時間、噴孔からの位置等）に対応する値」が同進行に対応する前記取得されている前記第１番目部分混合指標値と等しい値に決定されるとより好適である。
一般に、燃料に対して筒内ガスが混ざり合っていく程度を表す値（例えば、空気過剰率等）は、計算負荷の大きいべき演算等を含む非常に複雑な実験式等を使用して算出されていく。従って、係る算出には多大な計算負荷がかかることになるから、係る算出回数を減らすことが望まれている。
上記構成は、係る観点に基づくものである。即ち、上記構成によれば、「第２番目以降に噴射される部分に対して筒内ガスが混ざり合っていく程度を表す値」が上記第１番目部分混合指標値に基づいて決定されるから、「第２番目以降に噴射される部分に対して筒内ガスが混ざり合っていく程度を表す値」を上記複雑な実験式等を使用して計算する必要がなくなる。この結果、噴射燃料の各部分に基づいて形成されていくそれぞれの混合気間の不均一性が考慮されつつ、少ない計算負荷でそれぞれの混合気の状態が、簡易、且つ個別に推定され得る。
また、本発明に係るエミッション発生量推定装置は、上記本発明に係る混合気状態推定装置の何れかにより推定された対応する混合気の状態に基づいて、前記それぞれの混合気の燃焼に起因して発生する有害物質であるエミッションの発生量を個別に推定するエミッション発生量推定手段を備えている。ここにおいて、上記エミッションには、近年、排出量の低減が強く望まれているスート、及びＮＯｘ等を少なくとも含むことが好適である。
これによると、個別に推定された複数の混合気の状態のうち対応する混合気の状態に基づいてそれぞれの混合気の燃焼に起因して発生するエミッション（スート、ＮＯｘ等）の発生量が個別に推定される。従って、上述したエミッションの発生度合の不均一性が考慮されてエミッションの総発生量が精度良く推定され得る。
また、本発明に関連するエミッション発生量推定装置は、内燃機関の燃焼室内にて噴射された燃料が同燃焼室内に吸入されている筒内ガスと混ざり合って形成される混合気が着火した後において同燃焼室内にて定常火炎が発生しているとの仮定のもと、同定常火炎内の状態に基づいて、同定常火炎内における燃料が過剰な領域内の燃焼に起因して発生する有害物質であるエミッションの発生量と、同定常火炎内における酸素が過剰な領域内の燃焼に起因して発生する同エミッションの発生量とを個別に推定するエミッション発生量推定手段を備えている。
燃料の噴射期間が比較的長いと、混合気が着火した後において燃焼室内にて所謂定常火炎（或いは、定常火炎に非常に近い形態を有する火炎）が発生する場合がある。この場合、後述するように、定常火炎が占める領域を、同定常火炎内における燃料が過剰な領域（従って、定常酸素濃度が「０」となる領域）と、同定常火炎内における酸素が過剰な領域（従って、定常燃料濃度が「０」となる領域）とに区分することができる。
他方、後述するように、スートの発生量（より具体的には、スートの発生速度）の算出に使用される代表的な公知の実験式の一つには、燃料濃度の値を含むスートの生成速度（従って、増加速度）を求める項と、酸素濃度の値を含むスートの酸化速度（従って、減少速度）を求める項からなる２つの項が存在する。
従って、上記構成のように、定常火炎が占める領域を、定常酸素濃度が「０」となる領域と、定常燃料濃度が「０」となる領域とに区分すれば、それぞれの領域内におけるスートの発生量を上記代表的な実験式に従って計算する場合において、それぞれの領域について上記２つの項の何れかの値が常に「０」に維持される。従って、それぞれの領域内におけるスートの発生量の計算に際し、上記２つの項の何れかの計算を省略することができるから、スートの発生量の計算に関する計算負荷を小さくすることができる。
この場合、前記エミッション発生量推定手段は、前記それぞれの混合気、或いは前記定常火炎内のそれぞれの領域に対応する前記エミッションの発生速度を個別に算出し、同算出されたエミッションの発生速度を時間積分することで前記エミッションの発生量を個別に推定するように構成することが好適である。
一般に、エミッションの発生量の計算では、エミッション発生速度を求める公知の実験式等（例えば、上述したスートの発生速度を求める実験式等）が使用される。従って、上記構成によれば、十分に信頼の高い公知の実験式等を用いてエミッションの発生量が計算され得る。
また、この場合、前記エミッション発生量推定手段は、前記混合気の温度、或いは前記定常火炎内の前記領域内の温度が所定温度未満となったとき、対応する前記エミッションの発生速度の時間積分処理を終了するように構成することが好適である。
一般に、エミッションは、混合気の温度が反応限界温度未満になると殆ど発生しなくなる。従って、混合気の温度が対応する反応限界温度未満になった時点以降については、上記時間積分処理を継続するか否かにかかわらず同時間積分処理により得られるエミッションの総発生量は殆ど同一の値となる。従って、上記のように、混合気の温度（或いは、定常火炎内の領域内の温度）が所定温度未満（例えば、反応限界温度未満）となったとき、対応するエミッションの発生速度の時間積分処理を終了するように構成すれば、エミッションの発生量の計算に必要な計算負荷を低減することができる。
以下、本発明による内燃機関（ディーゼル機関）の混合気状態推定装置を含んだエミッション発生量推定装置の各実施形態について図面を参照しつつ説明する。
（第１実施形態）
図１は、本発明の第１実施形態に係る内燃機関のエミッション発生量推定装置を４気筒内燃機関（ディーゼル機関）１０に適用したシステム全体の概略構成を示している。このシステムは、燃料供給系統を含むエンジン本体２０、エンジン本体２０の各気筒の燃焼室（筒内）にガスを導入するための吸気系統３０、エンジン本体２０からの排ガスを放出するための排気系統４０、排気還流を行うためのＥＧＲ装置５０、及び電気制御装置６０を含んでいる。
エンジン本体２０の各気筒の上部には燃料噴射弁（噴射弁、インジェクタ）２１が配設されている。各燃料噴射弁２１は、図示しない燃料タンクと接続された燃料噴射用ポンプ２２に燃料配管２３を介して接続されている。燃料噴射用ポンプ２２は、電気制御装置６０と電気的に接続されていて、同電気制御装置６０からの駆動信号（後述する指令最終燃料噴射圧力Pcrfinに応じた指令信号）により燃料の実際の噴射圧力（吐出圧力）が同指令最終燃料噴射圧力Pcrfinになるように同燃料を昇圧するようになっている。
これにより、燃料噴射弁２１には、燃料噴射用ポンプ２２から前記指令最終燃料噴射圧力Pcrfinまで昇圧された燃料が供給されるようになっている。また、燃料噴射弁２１は、電気制御装置６０と電気的に接続されていて、同電気制御装置６０からの駆動信号（指令燃料噴射量（質量）Qfinに応じた指令信号）により噴射期間TAUだけ開弁し、これにより各気筒の燃焼室内に前記指令最終燃料噴射圧力Pcrfinにまで昇圧された燃料を前記指令燃料噴射量Qfinだけ直接噴射するようになっている。
吸気系統３０は、エンジン本体２０の各気筒の燃焼室にそれぞれ接続された吸気マニホールド３１、吸気マニホールド３１の上流側集合部に接続され同吸気マニホールド３１とともに吸気通路を構成する吸気管３２、吸気管３２内に回動可能に保持されたスロットル弁３３、電気制御装置６０からの駆動信号に応答してスロットル弁３３を回転駆動するスロットル弁アクチュエータ３３ａ、スロットル弁３３の上流において吸気管３２に順に介装されたインタクーラー３４と過給機３５のコンプレッサ３５ａ、及び吸気管３２の先端部に配設されたエアクリーナ３６とを含んでいる。
排気系統４０は、エンジン本体２０の各気筒にそれぞれ接続された排気マニホールド４１、排気マニホールド４１の下流側集合部に接続された排気管４２、排気管４２に配設された過給機３５のタービン３５ｂ、及び排気管４２に介装されたディーゼルパティキュレートフィルタ（以下、「ＤＰＮＲ」と称呼する。）４３を含んでいる。排気マニホールド４１及び排気管４２は排気通路を構成している。
ＥＧＲ装置５０は、排気ガスを還流させる通路（ＥＧＲ通路）を構成する排気還流管５１と、排気還流管５１に介装されたＥＧＲ制御弁５２と、ＥＧＲクーラー５３とを備えている。排気還流管５１はタービン３５ｂの上流側排気通路（排気マニホールド４１）とスロットル弁３３の下流側吸気通路（吸気マニホールド３１）を連通している。ＥＧＲ制御弁５２は電気制御装置６０からの駆動信号に応答し、再循環される排気ガス量（排気還流量、ＥＧＲガス流量）を変更し得るようになっている。
電気制御装置６０は、互いにバスで接続されたＣＰＵ６１、ＣＰＵ６１が実行するプログラム、テーブル（ルックアップテーブル、マップ）、及び定数等を予め記憶したＲＯＭ６２、ＣＰＵ６１が必要に応じてデータを一時的に格納するＲＡＭ６３、電源が投入された状態でデータを格納するとともに同格納したデータを電源が遮断されている間も保持するバックアップＲＡＭ６４、並びにＡＤコンバータを含むインターフェース６５等からなるマイクロコンピュータである。
インターフェース６５は、吸気管３２に配置された熱線式エアフローメータ７１、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に設けられた吸気温センサ７２、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気管圧力センサ７３、クランクポジションセンサ７４、アクセル開度センサ７５、燃料噴射用ポンプ２２の吐出口の近傍の燃料配管２３に配設された燃料温度センサ７６、気筒毎にそれぞれ配置された着火時期取得手段としての筒内圧力センサ７７、及び、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気酸素濃度センサ７８と接続されていて、これらのセンサからの信号をＣＰＵ６１に供給するようになっている。また、インターフェース６５は、燃料噴射弁２１、燃料噴射用ポンプ２２、スロットル弁アクチュエータ３３ａ、及びＥＧＲ制御弁５２と接続されていて、ＣＰＵ６１の指示に応じてこれらに駆動信号を送出するようになっている。
熱線式エアフローメータ７１は、吸気通路内を通過する吸入空気の質量流量（単位時間当りの吸入空気量、単位時間あたりの新気量）を計測し、同質量流量Ga（空気流量Ga）を表す信号を発生するようになっている。吸気温センサ７２は、エンジン１０のシリンダ（即ち、燃焼室、筒内）に吸入されるガスの温度（即ち、吸気温度）を検出し、同吸気温度Tbを表す信号を発生するようになっている。吸気管圧力センサ７３は、エンジン１０のシリンダに吸入されるガスの圧力（即ち、吸気管圧力）を検出し、同吸気管圧力Pbを表す信号を発生するようになっている。
クランクポジションセンサ７４は、各気筒の絶対クランク角度を検出し、実クランク角度CAactを表すとともにエンジン１０の回転速度であるエンジン回転速度NEをも表す信号を発生するようになっている。アクセル開度センサ７５は、アクセルペダルＡＰの操作量を検出し、アクセル操作量Accpを表す信号を発生するようになっている。燃料温度センサ７６は、燃料配管２３を通過する燃料の温度を検出し、燃料温度Tcrを表す信号を発生するようになっている。各筒内圧力センサ７７は、燃焼室内のガスの圧力（従って、上記筒内ガスの圧力）を検出し、筒内ガス圧力Paを表す信号を発生するようになっている。この筒内圧力センサ７７は、後述するように着火時点を検出するためにのみ使用される。吸気酸素濃度センサ７８は、吸気中の酸素濃度を検出し、吸気酸素濃度RO2inを表す信号を発生するようになっている。
（混合気状態の推定方法の概要）
次に、上記のように構成された混合気状態推定装置を含んだエミッション発生量推定装置（以下、「本装置」と云う。）による混合気状態の推定方法について説明する。
図２は、或る一つの気筒のシリンダ内(筒内、燃焼室内)に吸気マニホールド３１からガスが吸入され、燃焼室内に吸入されたガスが排気マニホールド４１へ排出される様子を模式的に示した図である。図２に示したように、燃焼室内に吸入されるガス（従って、筒内ガス）には、吸気管３２の先端部からスロットル弁３３を介して吸入された新気と、排気還流管５１からＥＧＲ制御弁５２を介して吸入されたＥＧＲガスが含まれる。
吸入される新気量（質量）と吸入されるＥＧＲガス量（質量）の和に対するＥＧＲガス量の割合（即ち、ＥＧＲ率）は、運転状態に応じて電気制御装置６０（ＣＰＵ６１）により適宜制御されるスロットル弁３３の開度、及びＥＧＲ制御弁５２の開度に応じて変化する。
かかる新気、及びＥＧＲガスは、吸気行程において開弁している吸気弁Vinを介してピストンの下降に伴って燃焼室内に吸入されて筒内ガスとなる。筒内ガスは、ピストンが圧縮下死点に達する時点近傍で吸気弁Vinが閉弁することにより燃焼室内に密閉され、その後の圧縮行程においてピストンの上昇に伴って圧縮される。
そして、ピストンが圧縮上死点近傍に達すると（具体的には、後述する燃料噴射開始時期（クランク角度）CAinjが到来すると）、本装置は、前記指令燃料噴射量Qfinに応じた噴射期間TAUだけ燃料噴射弁２１を開弁することで燃料を燃焼室内に直接噴射する。この結果、燃料噴射弁２１の噴孔から噴射された（液体の）燃料は、圧縮により高温になっている筒内ガスから受ける熱により直ちに燃料蒸気になるとともに、時間の経過に伴って同筒内ガスを取り込みながら混合気となって燃焼室内において円錐状に拡散していく。
上述したように、指令燃料噴射量Qfinの燃料は、実際には、上記燃料噴射開始時期CAinjから噴射期間TAUだけ連続して噴射されるが、以下においては便宜上、指令燃料噴射量Qfinの燃料が燃料噴射開始時期CAinjにて一時（瞬時）に噴射されるものとして説明を続ける。
図３（ａ）は、燃料噴射弁２１の噴孔から、噴射期間TAUに相当する指令燃料噴射量（質量）Qfinの燃料が一時に噴射された時点（即ち、噴射後経過時間ｔ＝０）での質量Qfinの燃料（燃料蒸気）の様子を模式的に示した図である。図３（ｂ）は、その後の或る時点（任意の噴射後経過時間ｔ）での図３（ａ）に示した質量Qfinの燃料蒸気の様子を模式的に示した図である。
図３（ｂ）に示すように、質量Qfinの燃料蒸気は、燃料噴射開始時期CAinj（即ち、噴射後経過時間ｔ＝０）において噴射された後、噴霧角θをもって円錐状に拡散しながら筒内ガスを順次取り込んでいく。そして、質量Qfinの燃料蒸気は、任意の噴射後経過時間ｔにおいて、噴射後経過時間ｔの関数である質量Gの筒内ガス（以下、「混合気形成筒内ガス」と云うこともある。）と混ざり合って質量（Qfin＋G）の混合気となっているものと仮定する。
本装置は、この混合気の任意の噴射後経過時間ｔにおける状態を推定する。混合気の状態としては、後述するエミッションの発生量の推定に必要となる混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、混合気内の酸素濃度[O2]mix、混合気内の窒素濃度[N2]mixが推定される。以下、先ず、係る混合気の状態推定に必要となる、任意の噴射後経過時間ｔにおける空気過剰率λの取得方法について説明する。
<空気過剰率λの取得>
噴射後経過時間ｔにおける空気過剰率λは、下記(1)式に示すように定義される。下記(1)式において、stoichは、単位質量の燃料の燃焼に必要な筒内ガスの質量（以下、「筒内ガス理論空燃比stoich」と呼ぶ。）である。筒内ガス理論空燃比stoichの値は、吸気中の酸素濃度に応じて変化すると考えられるから、上記吸気酸素濃度RO2inを引数とする所定の関数に従って取得され得る。
Figure 0004424242
このように定義される空気過剰率λは、例えば、日本機械学会論文集 25-156(1959年)，820ページ 「ディーゼル機関の噴霧到達距離に関する研究」 和栗雄太郎，藤井勝，網谷竜夫，恒屋礼次郎 （以下、「非特許文献１」と称呼する。）にて紹介された実験式である下記(2)式、及び下記(3)式に基づいて噴射後経過時間ｔの関数として求めることができる。
Figure 0004424242
Figure 0004424242
上記(3)式において、ｔは上記噴射後経過時間であり、dλ/dtは噴射後経過時間ｔの関数である燃料希釈率である。また、ｃは収縮係数、ｄは燃料噴射弁２１の噴孔径、ρfは（液体の）燃料密度、Ｌは論理希釈ガス量であって、これらの各値は全て定数である。
上記(3)式において、ΔPは有効噴射圧力であって、上記最終燃料噴射圧力Pcrfinから噴射開始時点（即ち、噴射後経過時間ｔ＝０）での筒内ガス圧力Pg0を減じた値である。筒内ガス圧力Pg0は、圧縮行程（及び膨張行程）における筒内ガスの状態が吸気弁Vinの閉弁時（即ち、筒内ガスが密閉された時点。以下、「ＩＶＣ」と呼ぶ。）以降断熱変化するとの仮定のもと、下記(4)式に従って求めることができる。
Figure 0004424242
上記(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は筒内ガスのガス定数（本例では、一定）である。
Figure 0004424242
また、上記(3)式において、ρgは噴射後経過時間ｔにおける筒内ガス密度であって、前記筒内ガスの全質量Mgを、噴射後経過時間ｔにおける上記筒内容積Vg(CA)で除することで、噴射後経過時間ｔの関数として取得することができる。
以上のことから、有効噴射圧力ΔPと噴霧角θとを上述のようにして求めれば、噴射後経過時間ｔの値と同噴射後経過時間ｔの関数である筒内ガス密度ρgの値とを使用して、上記(3)式に従って燃料希釈率dλ/dtが噴射後経過時間ｔの関数として求められる。そして、噴射後経過時間ｔ＝０から微小時間Δｔ（例えば、０．１msec）毎に求めた燃料希釈率dλ/dtの値を上記(2)式に従って時間で積分（積算）していくことで噴射後経過時間ｔにおける空気過剰率λを噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。
なお、上記(3)式から取得される燃料希釈率dλ/dtの値は常に正の値となることから上記(2)式から取得される空気過剰率λの値は噴射後経過時間ｔの増大に従って増加していく。そうすると、上記(1)式から理解できるように、混合気形成筒内ガスの質量Gが噴射後経過時間ｔの増大に従って増加していく。このことは、噴射後の燃料蒸気が円錐状に拡散していくことに伴って燃料蒸気と混ざり合う（燃料蒸気が取り込む）筒内ガス（従って、混合気形成筒内ガス）の量が増大していくことに対応している。
<混合気温度Tmixの取得>
次に、上述のように取得される空気過剰率λの値を利用して任意の噴射後経過時間ｔにおける混合気温度Tmixを取得する方法について説明する。一般に、混合気の熱エネルギー（エンタルピ）Hmixは、混合気温度Tmixを用いて下記(6)式に従って表すことができる。
Figure 0004424242
上記(6)式において、Mmixは混合気の総質量（混合気質量）、Cmixは混合気の定圧比熱である。従って、混合気のエンタルピHmix、混合気質量Mmix、及び混合気の定圧比熱Cmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていく（更新していく）ことで、下記(7)式に従って混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。以下、先ず、混合気質量Mmixの求め方について説明する。
Figure 0004424242
<<混合気質量Mmix>>
上述したように、質量Qfinの燃料蒸気は、任意の噴射後経過時間ｔにおいて質量Gの混合気形成筒内ガスと混ざり合って質量（Qfin＋G）の混合気となっているから、任意の噴射後経過時間ｔにおける混合気質量Mmixは(Qfin＋G)である。ここで、上記(1)式より「G＝stoich・λ・Qfin」と表すことができるから、混合気質量Mmixは、空気過剰率λを用いて下記(8)式にて表すことができる。
Figure 0004424242
よって、上述したように噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得され得る空気過剰率λの値を上記(8)式に順次適用していくことで、混合気質量Mmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。
<<混合気の定圧比熱Cmix>>
次に、混合気の定圧比熱Cmixの求め方について説明する。一般に、混合気の定圧比熱Cmixは、同混合気内の酸素濃度[O2]mix、及び混合気温度Tmixに大きく依存すると考えられる。ここで、混合気内の酸素濃度[O2]mixは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。従って、混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができれば、下記(9)式に従って、混合気の定圧比熱Cmixを微小時間Δｔ毎に求めることができる。
Figure 0004424242
上記(9)式において、funcCmixは、混合気の酸素濃度[O2]mix、及び混合気温度Tmixを引数とする混合気の定圧比熱Cmixを求めるための関数である。なお、上記(9)式を使用して混合気の定圧比熱Cmixを微小時間Δｔ毎に求めていく際における[O2]mix，Tmixの引数値としては、それぞれ現時点（即ち、噴射後経過時間ｔ）よりも微小時間Δｔ前の値が使用される。
<<混合気のエンタルピHmix>>
次に、混合気のエンタルピHmixの求め方について説明する。いま、噴射後経過時間（ｔ−Δt）における混合気のエンタルピHmix(ｔ−Δｔ)が既知である場合において、同噴射後経過時間(ｔ−Δｔ)から噴射後経過時間ｔまでの微小時間Δｔの間における混合気のエンタルピの増加分ΔHmixについて考える。この混合気のエンタルピの増加分ΔHmixとしては、微小時間Δｔの間に混合気に新たに取り込まれる筒内ガスの熱エネルギーΔHgと、同微小時間Δｔの間において混合気内で発生する化学反応による反応熱Hrとが挙げられる。
先ず、上記筒内ガスの熱エネルギーΔHgは、下記(10)式により表すことができる。ここにおいて、gは微小時間Δｔの間に混合気に新たに取り込まれる筒内ガスの質量である。この質量gは、噴射後経過時間ｔにおける混合気形成筒内ガスの質量から噴射後経過時間（ｔ−Δｔ）における混合気形成筒内ガスの質量を減じた値である。従って、上述した関係「G＝stoich・λ・Qfin」を利用して下記(11)式により求めることができる。(11)式において、λ(ｔ)，λ(ｔ−Δｔ)はそれぞれ、噴射後経過時間ｔ，(ｔ−Δｔ)における空気過剰率であり、上記(2)式、(3)式から求めることができる。
Figure 0004424242
Figure 0004424242
また、上記(10)式において、Tgは、噴射後経過時間ｔにおける筒内ガスの温度であり、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(12)式に従って求めることができる。下記(12)式において、上述したように、TgivcはＩＶＣにおける筒内ガス温度であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。
Figure 0004424242
また、上記(10)式において、Cgは、噴射後経過時間ｔにおける筒内ガスの定圧比熱であり、混合気の定圧比熱Cmixを求める上記(9)式と同様、下記(13)式に従って求めることができる。下記(13)式において、funcCgは、吸気中の酸素濃度[O2]in、及び筒内ガス温度Tgを引数とする筒内ガスの定圧比熱Cgを求めるための関数である。
Figure 0004424242
なお、上記(13)式を使用して筒内ガスの定圧比熱Cgを微小時間Δｔ毎に求めていく際における[O2]inの引数値としては、吸気酸素濃度センサ７８により検出される上記吸気酸素濃度RO2inが使用される。また、上記筒内ガス温度Tgの引数値としては、現時点（即ち、噴射後経過時間ｔ）における値が使用される。以上により、上記(10)式の右辺の項の全てを求めることができるから、同(10)式に従って上記筒内ガスの熱エネルギーΔHgを求めることができる。
次に、上記微小時間Δｔの間において混合気内で発生する化学反応による反応熱Hrは、下記(14)式で表すことができる。下記(14)式において、Hfは所定の定数であり、qrは上記微小時間Δｔの間において混合気内で発生する化学反応による燃料消費量である。
Figure 0004424242
上記燃料消費量qrの対象となる化学反応としては、着火反応（熱炎反応）、低温酸化反応（冷炎反応）のみならず、その他の種々の化学反応が含まれる。この燃料消費量qrは、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、混合気温度Tmixに大きく依存すると考えられるから、下記(15)式に従って表すことができる。
Figure 0004424242
上記(15)式において、funcqrは、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、及び混合気温度Tmixを引数とする上記燃料消費量qrを求めるための関数である。混合気内の燃料濃度[Fuel]mixも、上記混合気内の酸素濃度[O2]mixと同様、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。上記(15)式を使用して上記燃料消費量qrを微小時間Δｔ毎に求めていく際における[O2]mix，[Fuel]mixの引数値としては、ぞれぞれ現時点（即ち、噴射後経過時間ｔ）よりも微小時間Δｔ前の値が使用される。
また、上記(15)式の混合気温度Tmixの引数値としては化学反応前の混合気温度Tpreが使用される。この化学反応前の混合気温度Tpreは、上記(11)式にて算出される質量gの筒内ガスが混合気に新たに取り込まれた後であって、且つ、噴射後経過時間(ｔ−Δｔ)からの上記微小時間Δｔの間における化学反応が発生する前の段階における混合気温度であり、下記(16)式に従って求めることができる。
Figure 0004424242
上記(16)式において、Mmix(ｔ−Δｔ)，Cmix(ｔ−Δｔ)はそれぞれ、噴射後経過時間(ｔ−Δｔ)における混合気質量、及び混合気の定圧比熱であって、上記(8)式、及び上記(9)式によりそれぞれ取得され得る。更には、噴射後経過時間(ｔ−Δｔ)における混合気のエンタルピHmix(ｔ−Δｔ)は既知である。以上により、上記化学反応前の混合気温度Tpreを求めることができる。従って、上記(15)式の右辺の引数値を全て取得することができるから、(15)式，(14)式に従って上記化学反応による反応熱Hrを求めることができる。
以上のことから、噴射後経過時間（ｔ−Δt）における混合気のエンタルピHmix(ｔ−Δｔ)が既知である場合において、同噴射後経過時間(ｔ−Δｔ)から噴射後経過時間ｔまでの微小時間Δｔの間における混合気のエンタルピの増加分ΔHmix(＝ΔHg＋Hr)が求められるから、噴射後経過時間ｔにおける混合気のエンタルピHmix(ｔ)（＝Hmix(ｔ−Δｔ)＋ΔHmix）を求めることができる。
更には、噴射後経過時間ｔ＝０における混合気は、筒内ガスが取り込まれる前の状態（即ち、燃料蒸気のみ）であるから（図３（ａ）を参照）、この時点での混合気のエンタルピHmix(0)は、下記(17)式にて求めることができる。ここにおいて、Cfは燃料（蒸気）の定圧比熱（ここでは、定数）である。
Figure 0004424242
また、Tfは、燃料蒸気そのものの温度であり、液体の燃料が噴射直後に燃料蒸気に変化する際の単位質量当たりの潜熱Qvaporを考慮して下記(18)式に従って求めることができる。下記(18)式において、Tcrは噴射後経過時間ｔ＝０において燃料温度センサ７６により検出される液体の燃料温度である。αcrは燃料が燃料噴射用ポンプ２２の吐出口近傍から燃料噴射弁２１までの燃料配管２３を通過する際の熱損失分を考慮するための補正係数である。
Figure 0004424242
従って、噴射後経過時間ｔ＝０における混合気のエンタルピHmix(0)も求めることができる。以上より、混合気のエンタルピHmixを、噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。
このようにして、混合気のエンタルピHmix、混合気質量Mmix、及び混合気の定圧比熱Cmixがそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求められるから、上記(7)式に従って混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。
<混合気内の燃料濃度[Fuel]mixの取得>
次に、混合気内の燃料濃度（質量濃度）[Fuel]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の燃料濃度[Fuel]mixは、上記(8)式により取得される噴射後経過時間ｔにおける混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する燃料の質量」の割合である。
ここで、「噴射後経過時間ｔにおいて混合気内に存在する燃料の質量」は、噴射後経過時間ｔ＝０にて噴射された燃料量（指令燃料噴射量Qfin）から、噴射から現時点（噴射後経過時間ｔ）までの間に化学反応で消費された燃料分を減じた値である。従って、噴射後経過時間ｔにおける混合気内の燃料濃度[Fuel]mixは、下記(19)式で表すことができる。
Figure 0004424242
上記(19)式において、「Σqr」は、噴射から現時点（噴射後経過時間ｔ）まで微小時間Δｔ毎に上記(15)式に従って逐次取得・更新されていくそれぞれの燃料消費量qrの和である。このように、燃料消費量qr、及び混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(19)式に従って混合気内の燃料濃度[Fuel]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。
<混合気内の酸素濃度[O2]mixの取得>
次に、混合気内の酸素濃度（質量濃度）[O2]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の酸素濃度[O2]mixは、噴射後経過時間ｔにおける上記混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する酸素の質量」の割合である。
ここで、「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」は、上述した噴射後経過時間ｔにおける混合気形成筒内ガスの質量Gから、噴射から現時点（噴射後経過時間ｔ）までの間に化学反応で消費された筒内ガス分を減じた値である。上記微小時間Δｔにおける燃料消費量qrの燃料と反応して同微小時間Δｔの間において混合気内で化学反応により消費される筒内ガスの消費量grは、下記(20)式にて表すことができる。
Figure 0004424242
従って、「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」は、「G−Σgr」と表すことができる。ここで、「Σgr」は、噴射から現時点（噴射後経過時間ｔ）まで微小時間Δｔ毎に上記(20)式に従って逐次取得・更新されていくそれぞれの筒内ガス消費量grの和である。
「噴射後経過時間ｔにおいて混合気内に存在する酸素の質量」は、上記「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」に、筒内ガス内の酸素濃度（従って、上記吸気中の酸素濃度[O2]in）を乗じることで求めることができる。以上のことから、噴射後経過時間ｔにおける混合気内の酸素濃度[O2]mixは、下記(21)式で表すことができる。
Figure 0004424242
このように、混合気形成筒内ガスの質量G、筒内ガス消費量gr、及び混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(21)式に従って混合気内の酸素濃度[O2]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。
<混合気内の窒素濃度[N2]mixの取得>
次に、混合気内の窒素濃度（質量濃度）[N2]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の窒素濃度[N2]mixは、噴射後経過時間ｔにおける上記混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する窒素の質量」の割合である。
ここで、筒内ガス中の窒素は不活性ガスであって噴射から現時点（噴射後経過時間ｔ）までの間に混合気内で化学反応により消費されることがない。従って、「噴射後経過時間ｔにおいて混合気内に存在する窒素の質量」は、噴射後経過時間ｔにおける混合気形成筒内ガスの質量Gに、筒内ガス内の窒素濃度（即ち、吸気中の窒素濃度[N2]in）を乗じることで求めることができる。以上のことから、噴射後経過時間ｔにおける混合気内の窒素濃度[N2]mixは、下記(22)式で表すことができる。
Figure 0004424242
このように、混合気形成筒内ガスの質量G、混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(22)式に従って混合気内の窒素濃度[N2]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。以上、エミッション発生量の推定に必要となる混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、混合気内の酸素濃度[O2]mix、及び混合気内の窒素濃度[N2]mixの取得方法について説明した。
（エミッション発生量の推定方法）
次に、本装置によるエミッション発生量の推定方法について説明する。エミッション発生量としては、本例では、スート（Soot）発生量と、ＮＯ発生量とが推定される。以下、先ず、スート発生量の推定について説明する。
<スート発生量>
スート発生量は、例えば、日本機械学会論文集(Ｂ編) 48巻432号「直接噴射式ディーゼル機関の燃焼モデルと性能予測」（以下、「非特許文献２」と称呼する。）にて紹介された実験式である下記(23)式から得られる、混合気内のスート濃度[Soot]mixの変化速度（以下、「スート発生速度d[Soot]mix/dt」と称呼する。）を利用して求めることができる。
Figure 0004424242
上記(23)式において、dmsf/dtは混合気内でのスート生成に伴う混合気内のスート濃度[Soot]mixの増加速度（以下、「スートの生成速度」と称呼する。）であり、dmso/dtは混合気内で生成されたスートの酸化に伴う混合気内のスート濃度[Soot]mixの減少速度（以下、「スートの酸化速度」と称呼する。）である。これらはそれぞれ、下記(24)式、(25)式に従って表される。
Figure 0004424242
Figure 0004424242
上記(24)式、(25)式において、Af,Aoは定数、Esf,Esoは活性化エネルギー（ここでは、定数）である。また、Pgは、筒内ガスの圧力であり、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(26)式に従って求めることができる。下記(26)式において、上述したように、PgivcはＩＶＣにおける筒内ガス圧力であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。
Figure 0004424242
上記(24)式は、混合気内の燃料濃度[Fuel]mixが大きいほど、混合気温度Tmixが高いほど、筒内ガス圧力Pgが高いほど、スートの生成速度が大きくなる（従って、スートが生成され易くなる）ことを示している。また、上記(25)式は、混合気内のスート濃度[Soot]mixが大きいほど、混合気内の酸素濃度[O2]mixが大きいほど、混合気温度Tmixが高いほど、筒内ガス圧力Pgが高いほど、スートの酸化速度が大きくなる（従って、スートが酸化して消滅し易くなる）ことを示している。
上記(23)式〜(25)式に、上述したように噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次取得され得る混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、及び混合気内の酸素濃度[O2]mixを同微小時間Δｔ毎に順次適用していくことで、スート発生速度d[Soot]mix/dtを噴射後経過時間ｔ＝０から微小時間Δｔ毎に求めていくことができる。
従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に求められるスート発生速度d[Soot]mix/dtのそれぞれの値を時間で積分（積算）していくことで噴射後経過時間ｔにおける混合気内のスート濃度[Soot]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に更新していくことができる。このようにして、噴射後経過時間ｔにおける混合気内のスート濃度[Soot]mixを求めることができれば、この混合気内のスート濃度[Soot]mixの値に混合気質量Mmixの値を乗じることで噴射後経過時間ｔにおけるスート発生量を求めることができる。
<ＮＯ発生量>
次に、ＮＯ発生量の推定について説明する。ＮＯ発生量は、例えば、森北出版「燃焼工学」水谷幸夫著（以下、「非特許文献３」と称呼する。）にて紹介されている拡大ゼルドヴィッチ(Zel’dovich)機構に基づく式である下記(27)式から得られる、混合気内のＮＯ濃度[NO]mixの変化速度（以下、「ＮＯ発生速度d[NO]mix/dt」と称呼する。）を利用して求めることができる。
Figure 0004424242
上記(27)式において、Anは定数、Enは活性化エネルギー（ここでは、定数）である。また、Kfo(Tmix)は平衡定数であって混合気温度Tmixに基づいて決定される。上記(27)式は、混合気内の酸素濃度[O2]mixが大きいほど、混合気内の窒素濃度[N2]mixが大きいほど、混合気温度Tmixが高いほど、ＮＯの発生速度が大きくなる（従って、ＮＯが生成され易くなる）ことを示している。
上記(27)式に、上述したように噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次取得され得る混合気温度Tmix、混合気内の窒素濃度[N2]mix、及び混合気内の酸素濃度[O2]mixを同微小時間Δｔ毎に順次適用していくことで、ＮＯ発生速度d[NO]mix/dtを噴射後経過時間ｔ＝０から微小時間Δｔ毎に求めていくことができる。
従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に求められるＮＯ発生速度d[NO]mix/dtのそれぞれの値を時間で積分（積算）していくことで噴射後経過時間ｔにおける混合気内のＮＯ濃度[NO]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に更新していくことができる。このようにして、噴射後経過時間ｔにおける混合気内のＮＯ濃度[NO]mixを求めることができれば、この混合気内のＮＯ濃度[NO]mixの値に混合気質量Mmixの値を乗じることで噴射後経過時間ｔにおけるＮＯ発生量を求めることができる。
（噴射燃料の区分、混合気の区分）
以上のように取得される混合気の状態（温度Tmix等）、及びエミッション発生量は、図３（ａ），(ｂ)に示すように、指令燃料噴射量Qfinの燃料が燃料噴射開始時期CAinjにて一時（瞬時）に噴射されるとの仮定のもとで得られる値である。即ち、質量Qfinの燃料の全てが、燃料噴射開始時期CAinj以降、燃料噴射開始時期CAinjでの有効噴射圧力ΔP、同燃料噴射開始時期CAinjでの筒内ガス密度ρg0に基づいて決定される上記噴霧角θ（図３（ｂ）を参照）をもって「ひとつの混合気」となって円錐状に拡散していくものとして扱われる。
しかしながら、上述したように、指令燃料噴射量Qfinの燃料は、実際には、上記燃料噴射開始時期CAinjから噴射期間TAUだけ連続して噴射される。従って、燃焼室内にて進行していく混合気を微視的に見れば、混合気の末尾部分により近い部分に含まれる燃料ほど上記燃料噴射開始時期CAinjから遅れて噴射されることになる。換言すれば、混合気が占める領域内の位置によって、その位置に対応する燃料の部分の噴射時点が異なる。
他方、噴射時点が異なれば、噴射後の経過時間に対する筒内ガスの温度Tg、密度ρg等の状態が異なることに加え、噴射時点での有効噴射圧力ΔP、及び噴射時点での筒内ガス密度ρg0（従って、噴霧角θ）も異なる。換言すれば、燃焼室内にて進行していく混合気の状態は、同混合気に含まれる燃料の噴射時点によって異なる。
以上のことから、噴射後の経過時間に対する混合気の状態（温度等）は、同混合気が占める領域内の位置に応じて異なる。従って、混合気内でのエミッション発生量も同混合気が占める領域内の位置に応じて異なることになる。即ち、燃焼室内にて進行していく混合気の不均一性に起因して混合気内でのエミッション発生量も不均一になる。
係る混合気の不均一性、エミッション発生量の不均一性を考慮して混合気の状態、及びエミッション発生量を精度良く推定するためには、噴射期間TAUを複数の期間に区分して指令燃料噴射量Qfinの燃料を同区分された各期間に対応して順次噴射されていくそれぞれの部分に区分することが考えられる。そして、上記区分された燃料の各部分に基づいて形成されていくそれぞれの混合気の状態を個別に推定することが考えられる。
そこで、本例では、本装置は、図４に示すように、噴射期間TAUを「前期1/3TAU」「中期1/3TAU」「後期1/3TAU」の３つの期間に均等に区分し、前期1/3TAU、中期1/3TAU、後期1/3TAUに対応してそれぞれ質量Q(1)，Q(2)，Q(3)の燃料が順次個別に噴射されるものとする。
より具体的には、１番目噴射に係わる噴射期間「前期1/3TAU」分の質量Q(1)の燃料が燃料噴射開始時期CAinjにて一時に噴射され、１番目噴射から1/3TAU経過した後に２番目噴射に係わる噴射期間「中期1/3TAU」分の質量Q(2)の燃料が一時に噴射され、２番目噴射から1/3TAU経過した後に３番目噴射に係わる噴射期間「後期1/3TAU」分の質量Q(3)の燃料が一時に噴射されるものとする。なお、「Q(1)＋Q(2)＋Q(3)＝Qfin」の関係があるが、噴射期間TAUに亘る噴射圧力Pcrが一定に維持され得ないこと等に起因してQ(1)，Q(2)，Q(3)は互いに等しい値とはならない。
本装置は、１番目噴射に基づいて形成されていく混合気（１番目混合気）、２番目噴射に基づいて形成されていく混合気（２番目混合気）、及び３番目噴射に基づいて形成されていく混合気（３番目混合気）の「３つの混合気」をそれぞれ個別に扱い、混合気毎に、上述した手法により混合気の状態（温度Tmix等）及びエミッション発生量を個別に推定する。
これにより、混合気間での、噴射後の経過時間に対する筒内ガスの温度Tg、密度ρg等の状態の相違、噴射時点での有効噴射圧力ΔP及び筒内ガス密度ρg0（従って、噴霧角θ）の相違が考慮されて混合気の状態及びエミッション発生量が個別に推定され得る。
そして、本装置は、混合気毎に推定されたエミッション発生量を合計することでエミッションの総発生量（具体的には、スート総発生量Soot、ＮＯ総発生量ＮＯ）を推定する。これにより、上述したエミッション発生量の不均一性が考慮されてエミッションの総発生量が精度良く推定され得る。
（燃料噴射制御の概要）
本装置は、上述した混合気の状態の推定、及びエミッション発生量の推定に係わる計算を筒内ガスの量が確定するＩＶＣの直後に開始し、燃料噴射開始時期CAinjの到来前にエミッション総発生量（具体的には、スート総発生量Soot、及びＮＯ総発生量ＮＯ）の推定を完了する。
そして、本装置は、目標エミッション発生量（スート目標発生量Sootter、ＮＯ目標発生量ＮＯter）を機関の運転状態から求め、上記推定されたスート総発生量Soot、及びＮＯ総発生量ＮＯの一方が対応する目標発生量に対して十分に大きい場合、その一方の総発生量が小さくなるように燃料噴射圧力をフィードバック制御する。
具体的には、本装置は、スート総発生量Sootからスート目標発生量Sootterを減じた値が所定量よりも大きいとき、燃料噴射開始圧力を基本燃料噴射圧力Pcrbaseよりも所定量だけ高くする。これにより、スート発生量が小さくなる方向に制御される。一方、本装置は、ＮＯ総発生量ＮＯからＮＯ目標発生量ＮＯterを減じた値が所定量よりも大きいとき、燃料噴射開始圧力を基本燃料噴射圧力Pcrbaseよりも所定量だけ低くする。これにより、ＮＯ発生量が小さくなる方向に制御される。以上が、燃料噴射制御の概要である。
（実際の作動）
次に、上記のように構成された内燃機関のエミッション発生量推定装置の実際の作動について説明する。
<混合気温度等及びエミッション量の算出>
ＣＰＵ６１は、図５〜図９に一連のフローチャートにより示した混合気温度等、及びエミッション量の算出を行うためのルーチンを所定時間の経過毎に、気筒毎に、繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１はステップ５００から処理を開始し、ステップ５０５に進んで吸気弁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、及び図１１に示したテーブルMapCAinjから燃料噴射開始時期CAinjを決定する。テーブルMapCAinjは、指令燃料噴射量Qfin及びエンジン回転速度NEと燃料噴射開始時期CAinjとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。
続いて、ＣＰＵ６１はステップ５３０に進んで、指令燃料噴射量Qfin、エンジン回転速度NE、及び図１２に示したテーブルMapPcrbaseから基本燃料噴射圧力Pcrbaseを決定する。テーブルMapPcrbaseは、指令燃料噴射量Qfin及びエンジン回転速度NEと基本燃料噴射圧力Pcrbaseとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。
次に、ＣＰＵ６１はステップ５３５に進み、上記求めた燃料噴射期間TAUを３で除した値「TAU/3」と、上記求めた基本燃料噴射圧力Pcrbaseと、関数funcPcrとから、１番目噴射、２番目噴射、及び３番目噴射のそれぞれの燃料噴射圧力Pcr(1)，Pcr(2)，Pcr(3)を求める。上述のごとく、１番目噴射（質量Q(1)）は燃料噴射開始時期CAinjにて一時に実行され、２番目噴射（質量Q(2)）は１番目噴射から1/3TAU経過した後に一時に実行され、３番目噴射（質量Q(3)）は２番目噴射から1/3TAU経過した後に一時に実行されるものとする。
関数funcPcrは、噴射圧力が燃料噴射開始時期CAinjにおいて上記求められた基本燃料噴射圧力Pcrbaseに調整されている状態にて、燃料噴射開始時期CAinjから燃料噴射期間TAUに亘って燃料が連続して噴射される場合における噴射圧力の変動（低下）を考慮して、燃料噴射圧力Pcr(1)，Pcr(2)，Pcr(3)を決定する関数である。具体的には、これにより、１番目噴射の燃料噴射圧力Pcr(1)は基本燃料噴射圧力Pcrbaseと等しい値に設定され、２，３番目噴射圧力Pcr(2)，Pcr(3)は、基本燃料噴射圧力Pcrbaseよりも所定量だけ低い値に設定される。
続いて、ＣＰＵ６１はステップ５４０に進み、上記求めた燃料噴射期間TAU、値「TAU/3」と、上記求めた基本燃料噴射圧力Pcrbaseと、関数funcQとから、１番目噴射、２番目噴射、及び３番目噴射のそれぞれの燃料量（質量）Q(1)，Q(2)，Q(3)を求める。このステップ５４０が噴射燃料区分手段に相当する。
関数funcQは、噴射圧力が燃料噴射開始時期CAinjにおいて上記求められた基本燃料噴射圧力Pcrbaseに調整されている状態にて、燃料噴射開始時期CAinjから燃料噴射期間TAUに亘って燃料が連続して噴射される場合における噴射圧力の変動（低下）を考慮して、燃料量Q(1)，Q(2)，Q(3)を決定する関数である。これにより、燃料量（質量）Q(1)，Q(2)，Q(3)は、上述のごとく、「Qfin＝Q(1)＋Q(2)＋Q(3)」の関係が成立するように設定される一方で、互いに等しい値とはならない。
次に、ＣＰＵ６１はステップ５４５に進んで上記値「TAU/3」と、現時点でのエンジン回転速度NEと、上記決定された燃料噴射開始時期CAinjと、関数funcCAinjとから、１番目噴射、２番目噴射、及び３番目噴射のそれぞれの噴射時点でのクランク角度CAinj(1)，CAinj(2)，CAinj(3)を求める。
関数funcCAinjは、エンジン回転速度NEが現時点での値で一定であるものとしてクランク角度CAinj(1)，CAinj(2)，CAinj(3)を求める関数である。具体的には、これにより、１番目噴射時クランク角度CAinj(1)は上記燃料噴射開始時期CAinjに対応する値に設定され、２番目噴射時クランク角度CAinj(2)は燃料噴射開始時期CAinj（従って、１番目噴射時）から期間「TAU/3」だけ遅れた時期に対応する値に設定され、３番目噴射時クランク角度CAinj(3)は２番目噴射時から期間「TAU/3」だけ遅れた時期に対応する値に設定される。
次いで、ＣＰＵ６１はステップ５５０に進み、燃料温度センサ７６から得られる現時点での燃料温度Tcrと、上記(18)式とに基づいて燃料蒸気温度Tfを求める。続いて、ＣＰＵ６１はステップ５５５に進んで、上記吸気酸素濃度[O2]inと、[O2]inを引数とする筒内ガス理論空燃比stoichを求めるための関数funcstoichとに基づいて筒内ガス理論空燃比stoichを求める。
次に、ＣＰＵ６１はステップ５６０に進み、上記吸気酸素濃度[O2]inと、[O2]inを引数とする筒内ガス中の窒素濃度を求めるための関数func[N2]とに基づいて筒内ガス中の窒素濃度（即ち、吸気窒素濃度[N2]in）を求める。
続いて、ＣＰＵ６１はステップ５６５に進んで、現時点でのエンジン回転速度NEと、微小時間Δｔ（例えば、０．１msec）と、NE，Δｔを引数とする微小クランク角度ΔCAを求めるための関数funcΔCAとに基づいて、同微小時間Δｔに相当するクランク角度である微小クランク角度ΔCAを求める。この微小クランク角度ΔCAは、エンジン回転速度NEが現時点（即ち、ＩＶＣ直後）での値である場合における、微小時間Δｔに相当するクランク角度である。
そして、ＣＰＵ６１はステップ５７０に進み、スート総発生量Sootを初期値Soot0に設定するとともにＮＯ総発生量ＮＯを初期値ＮＯ0に設定し、続くステップ５７５にて変数iの値を「０」にクリアする。ここで、初期値Soot0は、燃料噴射前の時点から筒内ガス中に含まれている（従って、ＥＧＲガス中に予め含まれている）スート量に相当し、初期値ＮＯ0は、燃料噴射前の時点から筒内ガス中に含まれている（従って、ＥＧＲガス中に予め含まれている）ＮＯ量に相当する。
また、変数iの値は、何番目の噴射かを識別するための値であり、「i=1」は１番目噴射（従って、１番目混合気）を表し、「i=2」は２番目噴射（従って、２番目混合気）を表し、「i=3」は３番目噴射（従って、３番目混合気）を表す。
次に、ＣＰＵ６１は図６のステップ６０５へと進み、変数iの値をその時点での値（現時点では「０」）に「１」を加えた値に設定する。これにより、変数iは現時点では「１」になる。従って、以降、「i=1」である限りにおいて１番目噴射に係わる算出がなされていく。具体的には、先ず、以下のようにして、１番目噴射に関する各種初期値が決定されていく。
先ず、ＣＰＵ６１はステップ６１０に進んで、先のステップ５１５にて求めた筒内ガスの全質量Mgを、先のステップ５４５にて求めたi番目（１番目）噴射時クランク角度CAinj(i)から得られる１番目噴射時筒内容積Vg(CAinj(1))で除することで（１番目噴射時の）筒内ガス密度ρg0を求める。
続いて、ＣＰＵ６１はステップ６１５に進み、先のステップ５１０にて求めたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記i番目（１番目）噴射時筒内容積Vg(CAinj(i))と、上記(4)式に相当する式とに基づいて（１番目噴射時の）筒内ガス圧力Pg0を求める。
次に、ＣＰＵ６１はステップ６２０に進んで、先のステップ５３５にて求めたi番目（１番目）噴射圧力Pcr(i)から上記筒内ガス圧力Pg0を減じることで（１番目噴射時の）有効噴射圧力ΔPを求め、続くステップ６２５にて、上記求めた有効噴射圧力ΔPと、筒内ガス密度ρg0と、テーブルMapθとに基づいて（１番目噴射に係わる）噴霧角θを求める。これにより、噴霧角θは、i番目（１番目）噴射時点（即ち、クランク角度CAinj(1)）での有効噴射圧力ΔP及び筒内ガス密度ρg0に基づいて決定される。
次いで、ＣＰＵ６１はステップ６３０に進み、空気過剰率前回値λbを初期値「０」に設定し、続くステップ６３５にて（１番目混合気に係わる）混合気形成筒内ガス質量Gの値を初期値「０」に設定するとともに、続くステップ６４０にて（１番目混合気に係わる）燃料消費量積算値sumqr、及び筒内ガス消費量積算値sumgrを共に初期値「０」に設定する。
続いて、ＣＰＵ６１はステップ６４５に進んで、（１番目混合気に係わる）混合気のエンタルピHmixを、上記(17)式に相当する式に従って、初期値（即ち、先のステップ５４０にて求めたi番目（１番目）噴射量Q(i)と、燃料の定圧比熱Cfと、先のステップ５５０にて求めた燃料蒸気温度Tfの積）に設定する。
次に、ＣＰＵ６１はステップ６５０に進み、（１番目混合気に係わる）混合気の定圧比熱Cmixを初期値である上記燃料の定圧比熱Cfに設定し、続くステップ６５５にて（１番目混合気に係わる）混合気質量Mmixを初期値であるi番目（１番目）噴射量Q(i)に設定する。
次いで、ＣＰＵ６１はステップ６６０に進んで、（１番目混合気に係わる）混合気内のＮＯ濃度［NO］mix、スート濃度[Soot]mix、酸素濃度[O2]mixをそれぞれ初期値「０」に設定するとともに、（１番目混合気に係わる）混合気内の燃料濃度[Fuel]mixを初期値「１」に設定する。
続いて、ＣＰＵ６１はステップ６６５に進み、（１番目混合気に係わる）噴射後経過時間ｔを初期値「０」に設定するとともに、（１番目噴射に係わる）クランク角度CAを初期値であるi番目（１番目）噴射時クランク角度CAinj(i)に設定する。これにより、i番目（１番目）混合気に係わる噴射後経過時間ｔがi番目（１番目）噴射時からカウントされることになる。
次に、ＣＰＵ６１はステップステップ６７０に進んで、フラグＥＮＤsootの値、フラグＥＮＤnoの値を共に「０」に初期化する。ここで、フラグＥＮＤsootは、その値が「１」のときスート濃度[Soot]mixの値が更新中であることを示し、その値が「０」のときスート濃度[Soot]mixの更新が終了していることを示す。また、フラグＥＮＤnoは、その値が「１」のときＮＯ濃度[NO]mixの値が更新中であることを示し、その値が「０」のときＮＯ濃度[NO]mixの更新が終了していることを示す。このようにして、i番目（１番目）噴射に関する各種初期値が決定される。
次に、ＣＰＵ６１は図７のルーチンに進み、i番目（１番目）噴射に関する混合気温度の算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ７０５に進み、（１番目噴射に係わる）噴射後経過時間ｔの値（現時点では、「０」）を微小時間Δｔだけ増大・更新させるとともに、（１番目噴射に係わる）クランク角度CAの値（現時点では、「CAinj(1)」）を微小クランク角度ΔCAだけ増大・更新させる。このように、クランク角度CAの値が噴射後経過時間ｔに対応する値に維持されていく。これにより、以降、（１番目噴射に係わる）噴射後経過時間ｔ＝Δｔとなり、（１番目噴射に係わる）クランク角度CA＝CAinj(1)＋ΔCAとなる。
続いて、ＣＰＵ６１はステップ７１０に進んで、先のステップ５１５にて求めた筒内ガスの全質量Mgを、上記ステップ７０５にて更新されたクランク角度CAに対応する筒内容積Vg(CA)で除することで噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）筒内ガス密度ρgを求める。
次に、ＣＰＵ６１はステップ７１５に進んで、上記ステップ５１０にて求められたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記クランク角度CAに対応する筒内容積Vg(CA)と、上記(26)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）筒内ガス圧力Pgを求める。
次いで、ＣＰＵ６１はステップ７２０に進み、上記ステップ５１０にて求められたＩＶＣ時筒内ガス温度Tgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記クランク角度CAに対応する筒内容積Vg(CA)と、上記(12)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）筒内ガス温度Tgを求める。
次に、ＣＰＵ６１はステップ７２５に進んで、上記ステップ５１０にて求められた上記吸気酸素濃度[O2]inと、ステップ７２０にて求めた筒内ガス温度Tgと、[O2]in，Tgを引数とする筒内ガスの定圧比熱Cgを求めるための関数funcCgと、上記(13)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）筒内ガスの定圧比熱Cgを求める。
続いて、ＣＰＵ６１はステップ７３０に進み、先のステップ７１０にて求めた筒内ガス密度ρgと、先のステップ６２５にて求めた噴霧角θと、先のステップ６２０にて求めた有効噴射圧力ΔPと、先のステップ７０５にて更新したi番目（１番目）噴射に係わる噴射後経過時間ｔと、上記(3)式とに基づいてi番目（１番目）噴射に係わる燃料希釈率dλ/dtを求める。
次いで、ＣＰＵ６１はステップ７３５に進んで、上記(2)式に従って、i番目（１番目）噴射に係わる空気過剰率λを、その時点での空気過剰率前回値λb（現時点では、ステップ６３０の処理により「０」）に、上記求めた燃料希釈率dλ/dtに微小時間Δｔを乗じた値「dλ/dt・Δｔ」を加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）空気過剰率λが求められる。
次に、ＣＰＵ６１はステップ７４０に進み、先のステップ５５５にて求めた筒内ガス理論空燃比stoichと、ステップ７３５にて求めた空気過剰率λと、空気過剰率前回値λb（現時点では、ステップ６３０の処理により「０」。次回からは、後述するステップ７８５にて設定されている値）と、ステップ５４０にて設定されたi番目（１番目）噴射量Q(i)と、上記(11)式に相当する式とに基づいて（１番目混合気に係わる）微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気に新たに取り込まれた筒内ガス質量gを求める。
続いて、ＣＰＵ６１はステップ７４５に進んで、混合気形成筒内ガス質量Gを、その時点での値（現時点では、ステップ６３５の処理により「０」）に上記新たに取り込まれた筒内ガス質量gを加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気形成筒内ガス質量Gが求められる。
次に、ＣＰＵ６１はステップ７５０に進み、混合気質量Mmixを、その時点での値（現時点では、ステップ６５５の処理により１番目噴射量Q(1)）に上記新たに取り込まれた筒内ガス質量gを加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気質量Mmixが求められる。
続いて、ＣＰＵ６１はステップ７５５に進んで、化学反応前の混合気のエンタルピHpreを、その時点での混合気のエンタルピHmix（現時点では、ステップ６４５の処理により値「Q(1)・Cf・Tf」）に、上記(10)式に従って求められる「上記新たに取り込まれた筒内ガスの熱エネルギーΔHg＝g・Cg・Tg」を加えた値に設定する。
次に、ＣＰＵ６１はステップ７６０に進み、上記求めた化学反応前の混合気のエンタルピHpreを、上記ステップ７５０にて求めた混合気質量Mmixにその時点での混合気の定圧比熱Cmix（現時点では、ステップ６５０の処理により燃料の定圧比熱Cf。次回からは、後述するステップ８３０にて設定されている値）を乗じた値で除することにより上記(16)式に相当する式に従って化学反応前の混合気の温度Tpreを求める。
次いで、ＣＰＵ６１はステップ７６５に進んで、その時点での混合気内の酸素濃度[O2]min（現時点では、ステップ６６０の処理により「０」。次回からは、後述するステップ８２０にて設定されている値）と、燃料濃度[Fuel]mix（現時点では、ステップ６６０の処理により「０」。次回からは、後述するステップ８１５にて設定されている値）と、上記求めた化学反応前の混合気の温度Tpreと、上記(15)式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrを求める。
続いて、ＣＰＵ６１はステップ７７０に進み、上記求めた消費燃料量qrと、上記(14)式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による反応熱Hrを求め、続くステップ７７５にて混合気のエンタルピHmixを、上記求めた化学反応前の混合気のエンタルピHpreに上記求めた反応熱Hrを加えた値に設定・更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気のエンタルピHmixが求められる。
そして、ＣＰＵ６１はステップ７８０に進んで、上記ステップ７７５にて求めた混合気のエンタルピHmixと、上記ステップ７５０にて求めた混合気質量Mmixと、その時点での混合気の定圧比熱Cmix（現時点では、ステップ６５０の処理により燃料の定圧比熱Cf。次回からは、後述するステップ８３０にて設定されている値）と、上記(7)式とに基づいて混合気温度Tmixを算出する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj(1)＋ΔCA）における（１番目噴射に係わる）混合気温度Tmixが求められる。このステップ７８０は混合気状態推定手段に相当する。
次に、ＣＰＵ６１はステップ７８５に進んで、次回の計算の準備のため、空気過剰率前回値λbを上記ステップ７３５にて求めた空気過剰率λの値に設定する。この値は、以降、上述したステップ７３５にて使用される。このようにして、噴射後経過時間ｔにおけるi番目（１番目）噴射に関する混合気温度Tmixが算出される。
次に、ＣＰＵ６１は図８のルーチンに進み、i番目（１番目）噴射に関する各種濃度の算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ８０５に進み、上記ステップ７６５にて求めた微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrと、ステップ５５５にて求めた筒内ガス理論空燃比stoichと、上記(20)式とに基づいて、微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による筒内ガス消費量grを求める。
続いて、ＣＰＵ６１はステップ８１０に進んで、（１番目混合気に係わる）燃料消費量積算値sumqrを、その時点での値（現時点では、ステップ６４０の処理により「０」）にステップ６７５にて求めた上記燃料消費量qrを加えた値に設定・更新し、（１番目混合気に係わる）筒内ガス消費量積算値sumgrを、その時点での値（現時点では、ステップ６４０の処理により「０」）にステップ８０５にて求めた上記筒内ガス消費量grを加えた値に設定・更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）燃料消費量積算値sumqr、及び筒内ガス消費量積算値sumgrが求められる。
次いで、ＣＰＵ６１はステップ８１５に進み、ステップ５４０にて求めたi番目（１番目）噴射量Q(i)と、上記求めた燃料消費量積算値sumqrと、ステップ７５０にて求めた混合気質量Mmixと、上記(19)式に相当する式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気内の燃料濃度[Fuel]mixを求める。
次に、ＣＰＵ６１はステップ８２０に進んで、ステップ７４５にて求めた混合気形成筒内ガス質量Gと、上記求めた筒内ガス消費量積算値sumgrと、ステップ５１０にて設定された吸気酸素濃度[O2]inと、ステップ７５０にて求めた混合気質量Mmixと、上記(21)式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気内の酸素濃度[O2]mixを求める。
続いて、ＣＰＵ６１はステップ８２５に進み、ステップ７４５にて求めた混合気形成筒内ガス質量Gと、ステップ５６０にて設定された吸気窒素濃度[N2]inと、ステップ７５０にて求めた混合気質量Mmixと、上記(22)式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気内の窒素濃度[N2]mixを求める。このステップ８１５〜８２５も混合気状態推定手段に相当する。
次に、ＣＰＵ６１はステップ８３０に進んで、上記ステップ８２０にて求めた混合気内の酸素濃度[O2]mixと、ステップ７８０にて求めた混合気温度Tmixと、上記(9)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）混合気の定圧比熱Cmixを求める。この値は、以降、ステップ７６０、７８０にて使用される。
次いで、ＣＰＵ６１はステップ８３５に進み、フラグＥＮＤsootの値が「０」であるか否かを判定し、「Ｎｏ」と判定する場合、後述する混合気内のスート濃度[Soot]mixの更新（スート発生速度d[soot]mix/dtの時間積分処理。ステップ８５５）を行うことなくステップ８７０に直ちに進む。
現時点では、先のステップ６７０の処理によりフラグＥＮＤsootの値は「０」になっている。従って、ＣＰＵ６１はステップ８３５にて「Ｙｅｓ」と判定してステップ８４０に進み、ステップ８１５にて求めた混合気内の燃料濃度[Fuel]mixと、ステップ７１５にて求めた筒内ガス圧力Pgと、ステップ７８０にて求めた混合気温度Tmixと、上記(24)式とに基づいてスートの生成速度dmsf/dtを求める。
続いて、ＣＰＵ６１はステップ８４５に進んで、その時点での混合気内のスート濃度[Soot]mix（現時点では、ステップ６６０の処理により「０」。次回からは、後述するステップ８５５にて設定されている値）と、ステップ８２０にて求めた混合気内の酸素濃度[O2]mixと、ステップ７１５にて求めた筒内ガス圧力Pgと、ステップ７８０にて求めた混合気温度Tmixと、上記(25)式とに基づいてスートの酸化速度dmso/dtを求める。
次に、ＣＰＵ６１はステップ８５０に進み、上記求めたスートの生成速度dmsf/dtと、上記求めたスートの酸化速度dmso/dtと、上記(23)式とに基づいてスート発生速度d[soot]mix/dtを求め、続くステップ８５５にて、混合気内のスート濃度[Soot]mixを、その時点での値（現時点では、ステップ６６０の処理により「０」）に上記求めたスート発生速度d[soot]mix/dtに微小時間Δｔを乗じた値「d[Soot]mix/dt・Δｔ」を加えた値に設定・更新する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj(1)＋ΔCA）における（１番目噴射に係わる）混合気内のスート濃度[Soot]mixが求められる。このステップ８５５がエミッション発生量推定手段に相当する。
続いて、ＣＰＵ６１はステップ８６０に進んで、クランク角度CAが圧縮上死点（以下、「ＴＤＣ」と称呼する。）以降であって、且つ、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）ステップ７８０で求められた混合気温度Tmixがスートの反応限界温度TminSootより低いか否かを判定する。なお、スートの反応限界温度TminSootとは、混合気温度Tmixがその温度未満になるとスートが殆ど発生しなくなる温度である。
ここで、「Ｙｅｓ」と判定する場合、ＣＰＵ６１はステップ８６５に進んでフラグＥＮＤsootの値を「０」から「１」に変更する。これにより、以降、ＣＰＵ６１はステップ８３５に進んだとき「Ｎｏ」と判定するようになり、上述したように、混合気内のスート濃度[Soot]mixの更新（ステップ８５５）が実行されなくなる。これにより、スート濃度[Soot]mixの更新に係わる不必要な計算が省略でき、ＣＰＵ６１の計算負荷を低減することができる。
現時点はＩＶＣの直後であるから、ＴＤＣ前である。従って、ＣＰＵ６１はステップ８６０にて「Ｎｏ」と判定してステップ８７０に直ちに進む。即ち、フラグＥＮＤsootの値が「０」に維持される。
ＣＰＵ６１はステップ８７０に進むと、フラグＥＮＤnoの値が「０」であるか否かを判定する。ここで、「Ｎｏ」と判定する場合、ＣＰＵ６１は後述する混合気内のＮＯ濃度[NO]mixの更新（ＮＯ発生速度d[NO]mix/dtの時間積分処理。ステップ８８０）を行うことなく図９のステップ９０５に直ちに進む。
現時点では、先のステップ６７０の処理によりフラグＥＮＤnoの値は「０」になっている。従って、ＣＰＵ６１はステップ８７０にて「Ｙｅｓ」と判定してステップ８７５に進み、ステップ８２５にて求めた混合気内の窒素濃度[N2]mixと、ステップ８２０にて求めた混合気内の酸素濃度[O2]mixと、ステップ７８０にて求めた混合気温度Tmixと、上記(27)式とに基づいてＮＯ発生速度d[NO]mix/dtを求める。
続いて、ＣＰＵ６１はステップ８８０に進み、混合気内のＮＯ濃度[Soot]mixを、その時点での値（現時点では、ステップ６６０の処理により「０」）に上記求めたＮＯ発生速度d[NO]mix/dtに微小時間Δｔを乗じた値「d[NO]mix/dt・Δｔ」を加えた値に設定・更新する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj(1)＋ΔCA）における（１番目噴射に係わる）混合気内のＮＯ濃度[NO]mixが求められる。このステップ８８０もエミッション発生量推定手段に相当する。
続いて、ＣＰＵ６１はステップ８８５に進んで、クランク角度CAがＴＤＣ以降であって、且つ、噴射後経過時間ｔ（従って、クランク角度CA）における（１番目噴射に係わる）ステップ７８０で求められた混合気温度TmixがＮＯの反応限界温度TminNOより低いか否かを判定する。なお、ＮＯの反応限界温度TminNOとは、混合気温度Tmixがその温度未満になるとＮＯが殆ど発生しなくなる温度である。
ここで、「Ｙｅｓ」と判定する場合、ＣＰＵ６１はステップ８９０に進んでフラグＥＮＤnoの値を「０」から「１」に変更する。これにより、以降、ＣＰＵ６１はステップ８７０に進んだとき「Ｎｏ」と判定するようになり、上述したように、混合気内のＮＯ濃度[NO]mixの更新（ステップ８８０）が実行されなくなる。これにより、ＮＯ濃度[NO]mixの更新に係わる不必要な計算が省略でき、ＣＰＵ６１の計算負荷を低減することができる。
上述と同様、現時点はＩＶＣの直後であるから、ＴＤＣ前である。従って、現時点では、ＣＰＵ６１はステップ８８５にて「Ｎｏ」と判定して図９のルーチンに直ちに進む。即ち、フラグＥＮＤnoの値が「０」に維持される。以上により、変数i=1についての（従って、１番目噴射に係わる）噴射後経過時間ｔ＝Δｔ（クランク角度CA＝CAinj(1)＋ΔCA）における、混合気（１番目混合気）の空気過剰率λ、１番目混合気の状態（温度Tmix等）、及び１番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix）が算出される。
ＣＰＵ６１は図９のステップ９０５に進むと、フラグＥＮＤsootの値、及びフラグＥＮＤnoの値が共に「１」であるか、又は、クランク角度CAがＴＤＣ以降の所定の終了判定クランク角度CAendに一致したか否かを判定する。現時点では、上述のごとく、フラグＥＮＤsootの値、及びフラグＥＮＤnoの値が共に「０」であり、且つ、クランク角度CAは、１番目噴射時クランク角度CAinj(1)に微小クランク角度ΔCAを加えた値（従って、ＴＤＣ前）であるから上記終了判定クランク角度CAendに達していない。
従って、現時点では、ＣＰＵ６１はステップ９０５にて「Ｎｏ」と判定して図７のステップ７０５に戻り、（１番目噴射に係わる）噴射後経過時間ｔ（現時点では、「1・Δｔ」）を微小時間Δｔだけ増大・更新させるとともに、（１番目噴射に係わる）クランク角度CA（現時点では、「CAinj(1)＋ΔCA」）を微小クランク角度ΔCAだけ増大・更新させた後、上述した図７のステップ７１０〜図９のステップ９０５の処理を再び実行する。
これにより、変数i=1についての（従って、１番目噴射に係わる）噴射後経過時間ｔ＝２・Δｔ（クランク角度CA＝CAinj(1)＋２・ΔCA）における、混合気（１番目混合気）の空気過剰率λ（ステップ７３５を参照。）、１番目混合気の状態（温度Tmix等。ステップ７８０を参照。）、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）１番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix。ステップ８５５、８８０を参照。）が算出される。
そして、図９のステップ９０５の判定において「Ｎｏ」と判定される毎に、図７のステップ７０５〜図９のステップ９０５の処理が繰り返し実行されていく。これにより、図９のステップ９０５の判定において「Ｎｏ」と判定される限りにおいて、変数i＝１に維持された状態で、１番目混合気の空気過剰率λ、１番目混合気の状態（温度Tmix等）、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）１番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix）が１番目噴射時から微小時間Δｔ毎に（即ち、CAinj(1)から微小クランク角度ΔCA毎に）求められていく。
以降、ＴＤＣ後の膨張行程における筒内容積の増大等に伴って混合気温度Tmixが低下することで上述したステップ８６０の条件が成立した場合、フラグＥＮＤsootの値が「０」から「１」に変更されて、それ以降、上述のごとく、スート濃度[Soot]mixの更新（ステップ８５５）が実行されなくなる。また、膨張行程における筒内容積の増大等に伴って上述したステップ８８５の条件が成立した場合、フラグＥＮＤnoの値が「０」から「１」に変更されて、それ以降、上述のごとく、ＮＯ濃度[NO]mixの更新（ステップ８８０）が実行されなくなる。
そして、フラグＥＮＤsootの値とフラグＥＮＤnoの値が共に「１」となった場合、又は、（フラグＥＮＤsootの値とフラグＥＮＤnoの値が共に「１」となっていなくても）クランク角度CAが上記終了判定クランク角度CAendに達した場合、ＣＰＵ６１は図９のステップ９０５に進んだとき「Ｙｅｓ」と判定してステップ９１０以降に進み、変数i=1についての（従って、１番目噴射に係わる）計算を終了するための処理を行う。
即ち、ＣＰＵ６１はステップ９１０に進むと、先のステップ８５５の処理にて更新されている現時点での混合気内のスート濃度[Soot]mixの値に、先のステップ７５０にて更新されている現時点での混合気質量Mmixを乗じることでi番目（１番目）混合気内のスート発生量Soot(i)を求めるとともに、先のステップ８８０の処理にて更新されている現時点での混合気内のＮＯ濃度[NO]mixの値に、同混合気質量Mmixを乗じることでi番目（１番目）混合気内のＮＯ発生量NO(i)を求める。
次いで、ＣＰＵ６１はステップ９１５に進み、スート総発生量Sootを、その時点での値（現時点では、ステップ５７０の処理により初期値Soot0）に上記求めたi番目（１番目）混合気内のスート発生量Soot(i)を加えた値（Soot0＋Soot(1)）に更新するとともに、ＮＯ総発生量ＮＯを、その時点での値（現時点では、ステップ５７０の処理により初期値ＮＯ0）に上記求めたi番目（１番目）混合気内のＮＯ発生量NO(i)を加えた値（ＮＯ0＋NO(1)）に更新する。以上により、１番目噴射（従って、１番目混合気）に関する計算が終了する。
そして、ＣＰＵ６１はステップ９２０に進んで変数iの値が「３」であるか否かを判定する。現時点では、変数i＝1であるから、ＣＰＵ６１はステップ９２０にて「Ｎｏ」と判定して図６のステップ６０５に戻り、変数iの値を「１」だけ増大させる。これにより、以降、変数iが「２」に設定されるから、２番目噴射（従って、２番目混合気）に係わる算出がなされていく。
即ち、先ず、上述した図６のステップ６１０〜６７０にて、２番目噴射に関する各種初期値が設定される。具体的には、ステップ６１０にて、先のステップ５４５にて求めた２番目噴射時クランク角度CAinj(2)から得られる２番目噴射時筒内容積Vg(CAinj(2))を使用して（２番目噴射時の）筒内ガス密度ρg0が求められる。
続いて、ステップ６１５にて、上記２番目噴射時筒内容積Vg(CAinj(2))を使用して（２番目噴射時の）筒内ガス圧力Pg0が求められる。この結果、ステップ６２０にて、２番目噴射圧力Pcr(2)から上記筒内ガス圧力Pg0を減じることで（２番目噴射時の）有効噴射圧力ΔPが求められ、続くステップ６２５にて、上記求めた有効噴射圧力ΔPと、筒内ガス密度ρg0と、テーブルMapθとに基づいて（２番目噴射に係わる）噴霧角θが求められる。これにより、噴霧角θは、２番目噴射時点（即ち、クランク角度CAinj(2)）での有効噴射圧力ΔP及び筒内ガス密度ρg0に基づいて決定される。
また、ステップ６４５にて、（２番目混合気に係わる）混合気のエンタルピHmixが、初期値（即ち、先のステップ５４０にて求めた２番目噴射量Q(2)と、燃料の定圧比熱Cfと、先のステップ５５０にて求めた燃料蒸気温度Tfの積）に設定され、ステップ６５５にて、（２番目混合気に係わる）混合気質量Mmixが初期値である２番目噴射量Q(2)に設定される。
そして、ステップ６６５にて、（２番目混合気に係わる）噴射後経過時間ｔが初期値「０」に設定されるとともに、（２番目噴射に係わる）クランク角度CAが初期値である２番目噴射時クランク角度CAinj(2)に設定される。これにより、２番目混合気に係わる噴射後経過時間ｔが２番目噴射時からカウントされることになる。
このようにして、図６のステップ６１０〜６７０にて２番目噴射に関する各種初期値が設定されると、次に、上述した図７のステップ７０５〜図９のステップ９０５の処理が実行される。これにより、図９のステップ９０５の判定において「Ｎｏ」と判定される限りにおいて、変数i＝２に維持された状態で、２番目混合気の空気過剰率λ、２番目混合気の状態（温度Tmix等）、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）２番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix）が２番目噴射時から微小時間Δｔ毎に（即ち、CAinj(2)から微小クランク角度ΔCA毎）に求められていく。
そして、上述したステップ９０５の条件が成立すると、ステップ９１０にて２番目混合気内のスート発生量Soot(2)とＮＯ発生量NO(2)が求められ、ステップ９１５にて、スート総発生量Sootが、現時点での値（初期値Soot0＋Soot(1)）に上記求めた２番目混合気内のスート発生量Soot(2)を加えた値（Soot0＋Soot(1)＋Soot(2)）に更新されるとともに、ＮＯ総発生量ＮＯが、現時点での値（初期値ＮＯ0＋NO(1)）に上記求めた２番目混合気内のＮＯ発生量NO(2)を加えた値（ＮＯ0＋NO(1)＋NO(2)）に更新される。以上により、２番目噴射（従って、２番目混合気）に関する計算が終了する。
そして、ステップ９２０にて変数iの値が「３」であるか否かが判定される。現時点では、変数i＝２であるから、「Ｎｏ」と判定されて図６のステップ６０５に再び戻り、変数iの値が「１」だけ増大させられる。これにより、以降、変数iが「３」に設定されるから、３番目噴射（従って、３番目混合気）に係わる算出がなされていく。
即ち、先ず、上述した図６のステップ６１０〜６７０にて、３番目噴射に関する各種初期値が設定される。具体的には、ステップ６１０〜６２５にて（３番目噴射に係わる）噴霧角θが、３番目噴射時点（即ち、クランク角度CAinj(3)）での有効噴射圧力ΔP及び筒内ガス密度ρg0に基づいて決定される。
また、ステップ６４５にて、（３番目混合気に係わる）混合気のエンタルピHmixが、初期値（即ち、先のステップ５４０にて求めた３番目噴射量Q(3)と、燃料の定圧比熱Cfと、先のステップ５５０にて求めた燃料蒸気温度Tfの積）に設定され、ステップ６５５にて、（３番目混合気に係わる）混合気質量Mmixが初期値である３番目噴射量Q(3)に設定される。
そして、ステップ６６５にて、（３番目混合気に係わる）噴射後経過時間ｔが初期値「０」に設定されるとともに、（３番目噴射に係わる）クランク角度CAが初期値である３番目）噴射時クランク角度CAinj(3)に設定される。これにより、３番目混合気に係わる噴射後経過時間ｔが３番目噴射時からカウントされることになる。
このようにして、図６のステップ６１０〜６７０にて３番目噴射に関する各種初期値が設定されると、次に、上述した図７のステップ７０５〜図９のステップ９０５の処理が実行される。これにより、図９のステップ９０５の判定において「Ｎｏ」と判定される限りにおいて、変数i＝３に維持された状態で、３番目混合気の空気過剰率λ、３番目混合気の状態（温度Tmix等）、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）３番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix）が３番目噴射時から微小時間Δｔ毎に（即ち、CAinj(3)から微小クランク角度ΔCA毎）に求められていく。
そして、上述したステップ９０５の条件が成立すると、ステップ９１０にて３番目混合気内のスート発生量Soot(3)とＮＯ発生量NO(3)が求められ、ステップ９１５にて、スート総発生量Sootが、現時点での値（初期値Soot0＋Soot(1)＋Soot(2)）に上記求めた３番目混合気内のスート発生量Soot(3)を加えた値（Soot0＋Soot(1)＋Soot(2)＋Soot(3)）に更新されるとともに、ＮＯ総発生量ＮＯが、現時点での値（初期値ＮＯ0＋NO(1)＋NO(2)）に上記求めた３番目混合気内のＮＯ発生量NO(3)を加えた値（ＮＯ0＋NO(1)＋NO(2)＋NO(3)）に更新される。
以上により、３番目噴射（従って、３番目混合気）に関する計算が終了する。また、これにより、スート総発生量Sootが各混合気内のスート発生量の和として求められるとともに、ＮＯ総発生量ＮＯが各混合気内のＮＯ発生量の和として求められる。即ち、エミッション発生量の推定が終了する。ＣＰＵ６１は、以上説明した混合気状態の推定、及びエミッション発生量の推定をＩＶＣの直後（即ち、燃料噴射開始時期CAinjの前の時点）において完了する。
そして、ステップ９２０にて変数iの値が「３」であるか否かが判定される。現時点では、変数i＝３であるから、「Ｙｅｓ」と判定される。この場合、ＣＰＵ６１はステップ９２５に進み、現時点（即ち、ＩＶＣの直後の時点）におけるエンジン回転速度NEと、上記ステップ５２０にて決定されている指令燃料噴射量Qfinと、テーブルMapSootter，MapＮＯterとに基づいてスート目標発生量Sootter、及びＮＯ目標発生量ＮＯterを求める。
次いで、ＣＰＵ６１はステップ９３０に進んで、スート発生量偏差Δsootを、ステップ９１５にて求められているスート総発生量Sootから上記求めたスート目標発生量Sootterを減じた値に設定するとともに、ＮＯ発生量偏差ΔＮＯを、ステップ９１５にて求められているＮＯ総発生量ＮＯから上記求めたＮＯ目標発生量ＮＯterを減じた値に設定する。
続いて、ＣＰＵ６１はステップ９３５に進み、上記求めたスート発生量偏差Δsootが基準値Sootrefより大きいか否かを判定し、「Ｙｅｓ」と判定する場合、ステップ９４０に進んで最終燃料噴射圧力Pcrfinを、上記ステップ５３０にて決定されている基本燃料噴射圧力Pcrbaseに所定値ΔPcrを加えた値に設定する。これにより、スート発生量が小さくなる方向に燃料噴射圧力が補正されることになる。
一方、ステップ９３５にて「Ｎｏ」と判定する場合、ＣＰＵ６１はステップ９４５に進み、上記求めたＮＯ発生量偏差ΔＮＯが基準値ＮＯrefより大きいか否かを判定し、「Ｙｅｓ」と判定する場合、ステップ９５０に進んで最終燃料噴射圧力Pcrfinを、上記ステップ５３０にて決定されている基本燃料噴射圧力Pcrbaseから上記所定値ΔPcrを減じた値に設定する。これにより、ＮＯ発生量が小さくなる方向に燃料噴射圧力が補正されることになる。
他方、ステップ９４５にて「Ｎｏ」と判定する場合（即ち、ΔSoot≦Sootref、且つ、ΔＮＯ≦ＮＯrefの場合）、ＣＰＵ６１はステップ９５５に進んで、最終燃料噴射圧力Pcrfinを上記ステップ５３０にて決定されている基本燃料噴射圧力Pcrbaseと等しい値に設定する。即ち、この場合、燃料噴射圧力が補正されない。
そして、ＣＰＵ６１はステップ９６０に進むと、燃料噴射圧力が上記設定された最終燃料噴射圧力Pcrfinになるように燃料噴射用ポンプ２２（の駆動回路）に対して制御指示を行い、ステップ５９５に進んで図５〜図９の一連の本ルーチンを一旦終了する。以降、ＣＰＵ６１は、次のＩＶＣが到来するまでの間、ステップ５０５に進む毎に「Ｎｏ」と判定し続ける。
この結果、本ルーチンの実行により、ＩＶＣが到来する毎に、燃料噴射形態（噴射量、噴射圧力、噴射時期）が決定されるとともに混合気の状態、及びエミッション総発生量が直ちに推定され、係るエミッション総発生量の推定結果に基づいて噴射圧力が補正されていく。
また、ＣＰＵ６１は、図１３にフローチャートにより示した燃料噴射制御を行うためのルーチンを所定時間の経過毎に、気筒毎に、繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１はステップ１３００から処理を開始し、ステップ１３０５に進んで実クランク角度CAactが先のステップ５２５にて決定されている燃料噴射開始時期CAinjに一致したか否かを判定し、「Ｎｏ」と判定する場合、ステップ１３９５に直ちに進んで本ルーチンを一旦終了する。
いま、実クランク角度CAactが上記燃料噴射開始時期CAinjに一致したものとすると、ＣＰＵ６１はステップ１３１０に進んで、対応する燃料噴射弁２１に対してステップ５２０にて決定されている指令燃料噴射量Qfinの燃料の噴射指示（具体的には、燃料噴射期間TAUに亘る開弁指示）を行い、ステップ１３９５に進んで本ルーチンを一旦終了する。これにより、指令燃料噴射量Qfinの燃料が先のステップ９４０、９５０、９５５の何れかにて設定されている最終燃料噴射圧力Pcrfinをもって噴射される。
以上、説明したように、本発明による混合気状態推定装置、及びエミッション発生量推定装置の第１実施形態によれば、噴射期間TAUを「前期1/3TAU」「中期1/3TAU」「後期1/3TAU」の３つの期間に均等に区分し、「前期1/3TAU」に対応する１番目噴射（質量Q(1)）が燃料噴射開始時期CAinjにて一時に実行され、１番目噴射から1/3TAU経過した後に「中期1/3TAU」に対応する２番目噴射（質量Q(2)）が一時に実行され、２番目噴射から1/3TAU経過した後に「後期1/3TAU」に対応する３番目噴射（質量Q(3)）が一時に実行されるものとする。そして、１番目噴射に基づく１番目混合気、２番目噴射に基づく２番目混合気、及び３番目噴射に基づく３番目混合気の「３つの混合気」をそれぞれ個別に扱い、混合気毎に、噴射後経過時間ｔに対する混合気の空気過剰率λ（筒内ガスが混ざり合っていく程度を表す値）、混合気の状態（温度Tmix等）及び混合気内でのエミッション発生量（スート発生量、及びＮＯ発生量）を個別に推定する。
これにより、混合気間での、噴射後経過時間ｔに対する空気過剰率λの相違、噴射後経過時間ｔに対する筒内ガスの温度Tg、密度ρg等の状態の相違、噴射時点での有効噴射圧力ΔP及び筒内ガス密度ρg0（従って、噴霧角θ）の相違が考慮されて混合気の状態及びエミッション発生量が個別に推定され得る。そして、混合気毎に推定されたエミッション発生量を合計することでエミッションの総発生量（具体的には、スート総発生量Soot、ＮＯ総発生量ＮＯ）が推定される。これにより、上述した混合気の不均一性、エミッション発生量の不均一性が考慮されてエミッションの総発生量が精度良く推定され得る。
本発明は上記第１実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記第１実施形態では、指令燃料噴射量Qfin（従って、燃料噴射期間TAU）にかかわらず、噴射燃料（従って、混合気）を３つに区分しているが、指令燃料噴射量Qfinに応じて区分数を変更してもよい。この場合、指令燃料噴射量Qfinが大きいほど区分数を大きくすることが好適である。
また、上記第１実施形態においては、噴射期間TAUを複数の期間に均等に区分しているが、区分されたそれぞれの期間に対応する噴射燃料量が等しくなるように噴射期間TAUを複数の期間に区分してもよい。この場合、複数の期間は不均一になる。
また、上記第１実施形態においては、上記(15)式（図７のステップ７６５）にて算出される燃料消費量qrの対象となる化学反応として、着火反応（熱炎反応）、及び低温酸化反応（冷炎反応）のみならず、その他の種々の化学反応が含まれているが、着火反応、及び低温酸化反応に比してその他の種々の化学反応に係わる燃料消費量は十分に小さいことを考慮して、上記燃料消費量qrの対象となる化学反応として着火反応、及び低温酸化反応のみを扱うように構成してもよい。
この場合、燃料消費量qrを求めるための関数funcqrは、化学反応前の混合気温度Tpreが低温酸化反応が発生する温度範囲内にある場合に低温酸化反応により消費される燃料量の値を出力し、同化学反応前の混合気温度Tpreが着火反応が発生する温度範囲内にある場合には着火反応により消費される燃料量の値を出力し、同化学反応前の混合気温度Tpreがいずれの温度範囲内にもない場合には「０」を出力するように構成される。
また、上記第１実施形態においては、クランク角度CAがＴＤＣ以降であって、且つ、混合気温度Tmixがエミッションの反応限界温度より低くなった時点以降、ＣＰＵ６１の計算負荷低減のためにエミッション発生速度の時間積分処理が終了せしめられるように構成されているが、これに加え、クランク角度CAがＴＤＣ以前であって、且つ、混合気温度Tmixがエミッションの反応限界温度を超えるまでの間も、エミッション発生速度の時間積分処理を行わないように構成してもよい。これにより、エミッション発生量の算出に係わる不必要な計算が省略でき、より一層ＣＰＵ６１の計算負荷を低減することができる。
（第２実施形態）
次に、本発明の第２実施形態に係る内燃機関の混合気状態推定装置を含んだエミッション発生量推定装置について説明する。この第２実施形態は、燃料噴射期間TAUの途中において混合気が着火する場合を想定し、噴射燃料を混合気の着火前に噴射される部分と、同混合気の着火後に噴射される部分とに区分する点のみで上記第１実施形態と異なっている。従って、以下、係る相違点を中心に説明する。
（第２実施形態における噴射燃料の区分、混合気の区分）
燃料噴射期間TAUが比較的長い場合、燃料噴射期間TAUの途中において（即ち、燃料噴射継続中において）燃焼室内を進行していく混合気が着火する場合が多い。この場合、着火前に噴射された燃料に基づいて形成されていく混合気（以下、「着火前噴射に基づく混合気」と称呼する。）については、噴射から着火までの期間が比較的長いから、着火時点では既に混合気が十分に分散している。従って、予混合圧縮着火燃焼に類似した予混合的な燃焼が支配的になる場合が多い。
一方、着火後に噴射された燃料に基づいて形成されていく混合気（以下、「着火後噴射に基づく混合気」と称呼する。）については、噴射後直ちに着火することになるから、混合気が十分に分散していない状態で拡散しながら着火する。従って、拡散燃焼に類似した拡散的な燃焼が支配的になる場合が多い。
他方、予混合的燃焼と拡散的燃焼とでは燃料の反応速度（燃焼速度）が異なる。即ち、予混合的燃焼では、燃料と酸素とが十分に混ざり合った状態で着火が発生するから、燃料と容易に化学反応し得る酸素が十分に存在し得る。この結果、燃料の反応速度（燃焼速度）が比較的大きくなる。
これに対し、拡散的燃焼では、燃料と酸素が十分に混ざり合っていない状態で着火が発生するから、燃料と容易に化学反応し得る酸素が不足する。この結果、燃料の反応速度（燃焼速度）が比較的小さくなる。
更には、燃料の反応速度（燃焼速度）が異なると、上記(15)式にて表される微小時間Δｔの間において混合気内で発生する化学反応による燃料消費量qr（従って、反応熱Hr）も異なる。即ち、着火前噴射に基づく混合気と、着火後噴射に基づく混合気とでは、上記(15)式の関数funcqrの引数である、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、混合気温度Tmixが同じであっても、燃料消費量qrが異なることになる。
係る着火前噴射に基づく混合気と着火後噴射に基づく混合気の間の燃料消費量qr（従って、反応熱Hr）の相違を考慮して混合気の状態、及びエミッション発生量を精度良く推定するためには、上記(15)式に相当する関数funcqrを、着火前噴射に基づく混合気用（関数funcqrpre）と着火後噴射に基づく混合気用（関数funcqrpost）の２種類準備し、噴射燃料を混合気の着火前に噴射される部分と同混合気の着火後に噴射される部分とに区分することが考えられる。そして、着火前噴射に基づく混合気の状態を上記着火前噴射に基づく混合気用の関数funcqrpreを利用して得られる燃料消費量qrに基づいて推定するとともに、着火後噴射に基づく混合気の状態を上記着火後噴射に基づく混合気用の関数funcqrpostを利用して得られる燃料消費量qrに基づいて推定することが考えられる。
そこで、第２実施形態は、着火遅れ時間（噴射開始から着火までの時間）Tdelayを推定し、図１４に示すように、噴射期間TAUを「前期Tdelay」「後期(TAU−Tdelay)」の２つの期間に区分するとともに、前期Tdelay、後期(TAU−Tdelay)に対応してそれぞれ質量Q(1)，Q(2)の燃料が順次個別に噴射されるものとする。着火遅れ時間Tdelayの推定については後述する。
より具体的には、着火前噴射（１番目噴射）に係わる噴射期間「前期Tdelay」分の質量Q(1)の燃料が燃料噴射開始時期CAinjにて一時に噴射され、着火前噴射からTdelay経過した後に着火後噴射（２番目噴射）に係わる噴射期間「後期(TAU−Tdelay)」分の質量Q(2)の燃料が一時に噴射されるものとする。なお、「Q(1)＋Q(2)＝Qfin」の関係がある。
第２実施形態は、着火前噴射に基づく混合気（１番目混合気）の状態を求める際には上記着火前噴射に基づく混合気用の関数funcqrpreを利用して得られる燃料消費量qrを使用し、着火後噴射に基づく混合気（２番目混合気）の状態を求める際には上記着火後噴射に基づく混合気用の関数funcqrpostを利用して得られる燃料消費量qrを使用する。
そして、第２実施形態は、着火前噴射に基づく混合気と着火後噴射に基づく混合気とをそれぞれ個別に扱い、混合気毎に、上述した第１実施形態と同様の手法により混合気の状態（温度Tmix等）及びエミッション発生量を個別に推定する。そして、第２実施形態は、混合気毎に推定されたエミッション発生量を合計することでエミッションの総発生量（具体的には、スート総発生量Soot、ＮＯ総発生量ＮＯ）を推定する。
これにより、混合気間での、噴射後の経過時間に対する筒内ガスの温度Tg、密度ρg等の状態の相違、噴射時点での有効噴射圧力ΔP及び筒内ガス密度ρg0（従って、噴霧角θ）の相違に加え、混合気間での上述した燃料消費量qr（従って、反応熱Hr）の相違が考慮されて混合気の状態及びエミッション発生量が個別に推定され得る。即ち、上述した混合気の不均一性、エミッション発生量の不均一性に加えて上述した燃料消費量qr(従って、反応熱Hr)の相違が考慮されてエミッションの総発生量が精度良く推定され得る。
（第２実施形態の実際の作動）
以下、第２実施形態に係る混合気状態推定装置を含んだエミッション発生量推定装置の実際の作動について説明する。この装置のＣＰＵ６１は、第１実施形態のＣＰＵ６１が実行する図５〜図９、及び図１３に示したルーチンのうち、図６、図８、及び図１３に示したルーチンをそのまま実行するとともに、図５、図７、及び図９のルーチンをそれぞれ以下のように一部だけ変更して実行する。以下、主として図５、図７、及び図９のルーチンの変更部分について説明する。
第２実施形態のＣＰＵ６１は、図５に示したルーチンのステップ５３５〜５４５に代えて図１５に示したステップ１５０５〜１５２０を実行する。即ち、ＣＰＵ６１はステップ５３０からステップ１５０５に進むと、現時点から最も近い過去において燃焼が発生した気筒の筒内圧力センサ７７の出力経緯からその燃焼に係わる着火遅れ時間を取得し、同取得した着火遅れ時間を今回の燃焼についての着火遅れ時間Tdelayとする。
ここで、着火遅れ時間（噴射開始から着火までの時間）の取得は、着火発生時点にて筒内圧力が急増することに着目して筒内圧力センサ７７の出力経緯から着火時期を特定することで取得し得る。このようにして取得された着火遅れ時間を今回の燃焼についての着火遅れ時間Tdelayとするのは、今回の燃焼についての着火遅れ時間Tdelayは現時点から最も近い過去における燃焼に係わる着火遅れ時間に近い値になるであろうとの予測に基づく。
次に、ＣＰＵ６１はステップ１５１０に進んで、上記求めた着火遅れ時間Tdelayと、図５のステップ５３０にて求めた基本燃料噴射圧力Pcrbaseと、上記関数funcPcrとから、１番目噴射（着火前噴射）、及び２番目噴射（着火後噴射）のそれぞれの燃料噴射圧力Pcr(1)，Pcr(2)を求める。上述のごとく、着火前噴射（質量Q(1)）は燃料噴射開始時期CAinjにて一時に実行され、着火後噴射（質量Q(2)）は着火前噴射から着火遅れ時間Tdelay経過した後に一時に実行されるものとする。これにより、第１実施形態と同様、着火前噴射の燃料噴射圧力Pcr(1)は基本燃料噴射圧力Pcrbaseと等しい値に設定され、着火後噴射の燃料噴射圧力Pcr(2)は、基本燃料噴射圧力Pcrbaseよりも所定量だけ低い値に設定される。
続いて、ＣＰＵ６１はステップ１５１５に進み、図５のステップ５２０にて求めた燃料噴射期間TAU、着火遅れ時間Tdelayと、上記基本燃料噴射圧力Pcrbaseと、上記関数funcQとから、着火前噴射、及び着火後噴射のそれぞれの燃料量（質量）Q(1)，Q(2)を求める。これにより、燃料量（質量）Q(1)，Q(2)は、上述のごとく、「Qfin＝Q(1)＋Q(2)」の関係が成立するように設定される。このステップ１５１５が噴射燃料区分手段に相当する。
次いで、ＣＰＵ６１はステップ１５２０に進み、上記着火遅れ時間Tdelayと、現時点でのエンジン回転速度NEと、図５のステップ５２５にて決定された燃料噴射開始時期CAinjと、上記関数funcCAinjとから、着火前噴射、及び着火後噴射のそれぞれの噴射時点でのクランク角度CAinj(1)，CAinj(2)を求める。これにより、着火前噴射時クランク角度CAinj(1)は上記燃料噴射開始時期CAinjに対応する値に設定され、着火後噴射時クランク角度CAinj(2)は燃料噴射開始時期CAinj（従って、着火前噴射時）から着火遅れ時間Tdelayだけ遅れた時期に対応する値に設定される。そして、ＣＰＵ６１は図５のステップ５５０へと進む。
以上、ステップ１５０５〜１５２０の処理は、変数i＝１が着火前噴射（従って、着火前噴射に基づく混合気）に対応し、変数i＝２が着火後噴射（従って、着火後噴射に基づく混合気）に対応することを意味している。
また、第２実施形態のＣＰＵ６１は、図７に示したルーチンのステップ７６５に代えて図１６に示したステップ１６０５〜１６１５を実行する。即ち、ＣＰＵ６１はステップ７６０からステップ１６０５に進むと、変数iの値が「１」であるか否かを判定する。
いま、変数i＝１であるものとすると（即ち、着火前噴射に係わる算出がなされているものとすると）、ＣＰＵ６１はステップ１６０５にて「Ｙｅｓ」と判定してステップ１６１０に進み、その時点での混合気内の酸素濃度[O2]minと、燃料濃度[Fuel]mixと、図７のステップ７６０にて求めた化学反応前の混合気の温度Tpreと、上記(15)式に相当する着火前噴射に基づく混合気用の関数Funcqrpreとに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrを求める。
一方、変数i＝２であるものとすると（即ち、着火後噴射に係わる算出がなされているものとすると）、ＣＰＵ６１はステップ１６０５にて「Ｎｏ」と判定してステップ１６１５に進み、その時点での混合気内の酸素濃度[O2]minと、燃料濃度[Fuel]mixと、図７のステップ７６０にて求めた化学反応前の混合気の温度Tpreと、上記(15)式に相当する着火後噴射に基づく混合気用の関数Funcqrpostとに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrを求める。そして、ＣＰＵ６１は図７のステップ７７０へと進む。
また、第２実施形態のＣＰＵ６１は、図９に示したルーチンのステップ９２０に代えて図１７に示したステップ１７０５を実行する。即ち、ＣＰＵ６１はステップ９１５からステップ１７０５に進むと、変数iが「２」であるか否かを判定し、変数iが「２」でない場合（即ち、i＝１である場合）、図６のステップ６０５に戻る一方、変数iが「２」である場合、ステップ９２５以降へと進む。
これにより、ステップ１７０５にて「Ｙｅｓ」と判定された時点でエミッション発生量の推定が完了し、スート総発生量Sootが、初期値Soot0に、着火前噴射に基づく混合気内のスート発生量Soot(1)と、着火後噴射に基づく混合気内のスート発生量Soot(2)とを加えた値として求められる。同様に、ＮＯ総発生量ＮＯが、初期値ＮＯ0に、着火前噴射に基づく混合気内のＮＯ発生量NO(1)と、着火後噴射に基づく混合気内のＮＯ発生量NO(2)とを加えた値として求められる。
以上、説明したように、本発明による混合気状態推定装置、及びエミッション発生量推定装置の第２実施形態によれば、噴射期間TAUを「前期Tdelay」「後期(TAU−Tdelay)」の２つの期間に区分し、「前期Tdelay」に対応する着火前噴射（質量Q(1)）が燃料噴射開始時期CAinjにて一時に実行され、着火前噴射から着火遅れ時間Tdelay経過した後に「後期(TAU−Tdelay)」に対応する着火後噴射（質量Q(2)）が一時に実行されるものとする。
そして、着火前噴射に基づく１番目混合気、及び着火後噴射に基づく２番目混合気の「２つの混合気」をそれぞれ個別に扱い、混合気毎に、噴射後経過時間ｔに対する混合気の空気過剰率λ、混合気の状態（温度Tmix等）及び混合気内でのエミッション発生量（スート発生量、及びＮＯ発生量）を個別に推定する。そして、混合気毎に推定されたエミッション発生量を合計することでエミッションの総発生量が推定される。
加えて、着火前噴射に基づく混合気（１番目混合気）の状態を求める際には上記着火前噴射に基づく混合気用の関数funcqrpreを利用して得られる燃料消費量qrを使用し、着火後噴射に基づく混合気（２番目混合気）の状態を求める際には上記着火後噴射に基づく混合気用の関数funcqrpostを利用して得られる燃料消費量qrを使用する。
これにより、上述した混合気の不均一性、エミッション発生量の不均一性のみならず、上述した「着火前噴射に基づく混合気と着火後噴射に基づく混合気の間の燃料消費量qrの相違」が考慮されてエミッションの総発生量が精度良く推定され得る。
本発明は上記第２実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記第２実施形態では、噴射燃料を、混合気の着火前に噴射される部分と、同混合気の着火後に噴射される部分との「２つの部分」に区分するように構成されているが、混合気の着火前に噴射される部分を更に複数の（例えば、Ｍ個の）部分に区分し、混合気の着火後に噴射される部分も更に複数の（例えば、Ｎ個の）部分に区分することで噴射燃料を合計(Ｍ＋Ｎ)個の部分に区分するように構成してもよい。
この場合、上記Ｍ個の部分に基づくＭ個の混合気の状態を求める際には、上記着火前噴射に基づく混合気用の関数funcqrpreを利用して得られる燃料消費量qrが使用され、上記Ｎ個の部分に基づくＮ個の混合気の状態を求める際には、上記着火後噴射に基づく混合気用の関数funcqrpostを利用して得られる燃料消費量qrが使用される。
これによると、上述した「着火前噴射に基づく混合気と着火後噴射に基づく混合気の間の燃料消費量qrの相違」が考慮されつつ、上述した混合気の不均一性、エミッション発生量の不均一性がより一層明確に考慮されてエミッションの総発生量がより一層精度良く推定され得る。
（第３実施形態）
次に、本発明の第３実施形態に係る内燃機関の混合気状態推定装置を含んだエミッション発生量推定装置について説明する。この第３実施形態は、燃料噴射期間TAUを微小時間Δｔ（例えば、0.1msec）毎に（多数の（n個の）期間に）区分して噴射燃料を同区分された各期間に対応して噴射されるそれぞれの部分（n個の部分）に区分する点、並びに、２番目以降に噴射される部分に基づいて形成されていく混合気の空気過剰率を上記(2)式、(3)式を利用することなく１番目に噴射される部分に基づいて形成されていく混合気の空気過剰率に基づいて決定する点で、上記第１、第２実施形態と異なっている。以下、係る相違点を中心に説明する。なお、以下、i番目（i：n以下の自然数）に噴射される部分に基づいて形成されていく混合気を「i番目混合気」と称呼する。
上記第１、第２実施形態では、（２つ、又は３つに区分された）混合気毎に、且つ、微小時間Δｔ（例えば、0.1msec）毎に、上記(2)式、(3)式を利用して空気過剰率λが算出されている（図７のステップ７３０、７３５を参照）。これは、混合気毎に、且つ、微小時間Δｔ毎に、混合気に新たに取り込まれる筒内ガス量gを求める必要があるからである（図７のステップ７４０を参照）。
ここで、(3)式の右辺には、時々刻々と変化する変数ρg，tに関するべき演算が含まれている。べき演算には多大な計算負荷が伴う。従って、上記(2)式、(3)式を利用した空気過剰率λの算出には多大な計算負荷が伴うから、係る算出回数をできるだけ減らすことが望ましい。
しかしながら、噴射燃料を多数の（n個の）部分に区分するこの第３実施形態で上記第１、第２実施形態と同様に、混合気毎、且つ、微小時間Δｔ毎に上記(2)式、(3)式を利用した空気過剰率λの算出を行うものとすると、係る算出回数が非常に大きな値となってしまい、計算負荷が膨大となる。
そこで、第３実施形態では、２番目以降に噴射される部分に基づいて形成されていくそれぞれの混合気の噴射後経過時間ｔ＝０からの微小時間Δｔ毎の空気過剰率は、１番目混合気の噴射後経過時間ｔ＝０からの微小時間Δｔ毎の空気過剰率と等しいと仮定される。係る仮定により、２番目以降に噴射される部分に基づいて形成されていくそれぞれの混合気について、微小時間Δｔ毎に混合気に新たに取り込まれる筒内ガス量gを上記(2)式、(3)式を利用することなく算出することができる。この結果、２番目以降に噴射される部分に基づいて形成されていくそれぞれの混合気についての上記(2)式、(3)式を利用した空気過剰率λの算出が一切不要となる。以下、この点について図１８を参照しながら説明する。
図１８は、n個に区分された噴射燃料のそれぞれの部分（q(1),q(2),・・・,q(n−1),q(n)）に基づいて形成されていくn個の混合気の時間経過に対する様子（噴射開始から噴射終了までの様子）を模式的に示した図である。
図１８に示すように、第３実施形態では、１番目噴射に係わる噴射期間「１番目の微小時間Δｔ」分の質量q(1)の燃料が燃料噴射開始時期CAinjにて一時に噴射され、２番目噴射に係わる噴射期間「２番目の微小時間Δｔ」分の質量q(2)の燃料が燃料噴射開始時期CAinjからΔｔ経過した後に一時に噴射され、・・・、n番目噴射に係わる噴射期間「n番目の微小時間Δｔ」分の質量q(n)の燃料が燃料噴射開始時期CAinjから(n−1)・Δｔ経過した後に一時に噴射されるものとする。なお、指令燃料噴射量Qfinとq(i)（i：n以下の自然数）には、下記(28)式に示す関係が成立するが、燃料噴射期間TAUに亘る噴射圧力の変動等に起因してq(i)（i：n以下の自然数）は互いに等しい値とはならない。
Figure 0004424242
更には、図１８に示すように、i番目混合気についてi番目混合気に係わる噴射後経過時間ｔ＝(k−1)・Δｔ
〜 k・Δｔの間にi番目混合気に新たに取り込まれる筒内ガス量を「g(i,k)」と表すものとする（i：n以下の自然数、k：自然数。以下、同様）。
先ず、１番目噴射（従って、１番目混合気）のみについて考える。１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおける１番目混合気の空気過剰率λ(k)（第１番目部分混合指標値）は、上記第１、第２実施形態と同様、上記(2)式、(3)式を利用して順次求めることができる。
従って、１番目混合気に新たに取り込まれる筒内ガス量g(1,k)（図１８において斜線で示した部分に対応する）は、上記(11)式に相当する下記(29)式に従って求めることができる。なお、λ(0)＝0とする。即ち、１番目混合気に新たに取り込まれる筒内ガス量g(1,k)だけは、上記第１、第２実施形態と同様、上記(2)式、(3)式を利用して求められる。
Figure 0004424242
ここで、図１８、及び上記(29)式から理解できるように、１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔの時点で１番目混合気に取り込まれている筒内ガスの総量sumg(k)は、下記(30)式にて表される。
Figure 0004424242
例えば、１番目混合気に係わる噴射後経過時間ｔ＝２・Δｔの時点で１番目混合気に取り込まれている筒内ガスの総量sumg(2)＝g(1,1)＋g(1,2)であり、ｔ＝３・Δｔの時点ではsumg(3)＝g(1,1)＋g(1,2)＋g(1,3)となる。従って、１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔの空気過剰率λ(k)は、下記(31)式に従って表すことができる。
Figure 0004424242
従って、上記(31)から下記(32)式が導出され得るから、値「λ(k)−λ(k−1)」を、既知の値である「g(1,k)」と「q(1)」とを利用して求めることができる。
Figure 0004424242
ここで、上記仮定によれば、２番目以降のi番目混合気（i：２以上n以下の自然数）に係わる噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気の空気過剰率は、上記求められた「１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおける１番目混合気の空気過剰率λ(k)」と等しい。従って、２番目以降のi番目混合気についてi番目混合気に係わる噴射後経過時間ｔ＝(k−1)・Δｔ
〜 k・Δｔの間にi番目混合気に新たに取り込まれる筒内ガス量g(i,k)は、上述した仮定を利用すれば、上述した「１番目混合気に新たに取り込まれる筒内ガス量g(1,k)」を表す上記(29)式と同様に、下記(33)式にて表すことができる。
Figure 0004424242
この上記(33)式に上記(32)式を代入すると、下記(34)式を得ることができる。この下記(34)式によれば、２番目以降のi番目混合気に新たに取り込まれる筒内ガス量g(i,k)を、計算負荷の大きい上記(2)式、(3)式を利用することなく、既知の値である「g(1,k)」と「q(i)」と「q(1)」とを利用して簡易に求めることができる。即ち、２番目以降のi番目混合気に新たに取り込まれる筒内ガス量g(i,k)を求めるために、２番目以降のi番目混合気についての上記(2)式、(3)式を利用した空気過剰率λを算出する必要がなくなる。
Figure 0004424242
このように、第３実施形態では、２番目以降のi番目の混合気の噴射後経過時間ｔ＝０からの微小時間Δｔ毎の空気過剰率が、１番目混合気の噴射後経過時間ｔ＝０からの微小時間Δｔ毎の空気過剰率λ(k)と等しいと仮定することで、上記(2)式、(3)式を利用した空気過剰率λの算出を必要とする混合気を１番目混合気のみとすることができる。従って、計算負荷の大きい上記(2)式、(3)式を利用した空気過剰率λの算出回数を減らすことができ、ＣＰＵ６１の計算負荷を減らすことができる。
なお、この第３実施形態では、第１、第２実施形態と同様に混合気間での噴射後経過時間ｔに対する筒内ガスの状態（温度Tg、圧力Pg等）の相違は考慮され得る一方、混合気間での噴射時点での有効噴射圧力ΔP及び筒内ガス密度ρg0（従って、噴霧角θ）の相違、即ち、混合気間での噴射後経過時間ｔに対する空気過剰率λの相違が考慮されることなく、混合気の状態及びエミッション発生量が個別に推定されることになる。従って、第１、第２実施形態に比して上述した混合気の不均一性、及びエミッション発生量の不均一性が考慮される程度は低くなる。
（第３実施形態の実際の作動）
以下、第３実施形態に係る混合気状態推定装置を含んだエミッション発生量推定装置の実際の作動について説明する。この装置のＣＰＵ６１は、第１実施形態のＣＰＵ６１が実行する図５〜図９に示した一連のルーチン、及び図１３に示したルーチンのうち、図１３に示したルーチンのみをそのまま実行するとともに、図５〜図９に示した一連のルーチンに代えて図５〜図９のルーチンにそれぞれ対応する図１９〜図２３にフローチャートにより示した一連のルーチンを実行する。加えて、この装置のＣＰＵ６１は、図２４〜図２７にフローチャートにより示した一連のルーチンを追加的に実行する。
なお、図１９〜図２７に示したルーチンのステップであって図５〜図９に示したルーチンのステップと同一のものについては、図５〜図９に示したルーチンの対応するステップの番号と同一の番号を付すことでそれらの説明を省略する。以下、第３実施形態に特有の図１９〜図２７に示した各ルーチンについて説明する。
第３実施形態のＣＰＵ６１は、図５〜図９に示した一連のルーチンに対応する図１９〜図２３に示した一連のルーチンを所定時間の経過毎に繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１は図１９のステップ１９００から処理を開始し、ステップ５０５に進んで、「Ｙｅｓ」と判定する場合（即ち、ＩＶＣが到来した場合）、ステップ５１０〜５３０の処理を順に実行し、その後、ステップ１９０５に進んで、ステップ５２０にて求めた燃料噴射期間TAUを微小時間Δｔ（例えば、0.1msec）で除することで区分数nを求める。
次いで、ＣＰＵ６１はステップ１９１０に進み、図５のステップ５４０と同様にして、上記求めた燃料噴射期間TAU、微小時間Δｔと、上記求めた基本燃料噴射圧力Pcrbaseと、上記関数funcQとから、１番目噴射からn番目噴射までのそれぞれの燃料量（質量）q(1)，q(2)，・・・,q(n)を求める。これにより、燃料量（質量）q(1)，q(2),・・・,q(n)は、上述のごとく、「Qfin＝q(1)＋q(2)＋・・・＋q(n)」の関係が成立するように設定される一方で、互いに等しい値とはならない。このステップ１９１０は噴射燃料区分手段に相当する。
続いて、ＣＰＵ６１はステップ１９１５に進んで、図５のステップ５４５と同様にして、上記微小時間Δｔと、現時点でのエンジン回転速度NEと、ステップ５２５にて決定された燃料噴射開始時期CAinjと、上記関数funcCAinjとから、１番目噴射からn番目噴射までのそれぞれの噴射時点でのクランク角度CAinj(1)，CAinj(2)，・・・,CAinj(n)を求める。
これにより、１番目噴射時クランク角度CAinj(1)は上記燃料噴射開始時期CAinjに対応する値に設定され、２番目噴射時クランク角度CAinj(2)は燃料噴射開始時期CAinjからΔｔだけ遅れた時期に対応する値に設定され、・・・、n番目噴射時クランク角度CAinj(n)は燃料噴射開始時期CAinjから(n−1)・Δｔだけ遅れた時期に対応する値に設定される。
次に、ＣＰＵ６１はステップ５５０〜５７０の処理を順に実行し、その後、図２０のルーチンに進んで１番目噴射のみに関する各種初期値の決定を行うための処理を行う。具体的には、ＣＰＵ６１はステップ２００５に進んで、図６のステップ６１０と同様にして、ステップ５１５にて求めた筒内ガスの全質量Mgを、先のステップ１９１５にて求めた１番目噴射時クランク角度CAinj(1)から得られる１番目噴射時筒内容積Vg(CAinj(1))で除することで、１番目噴射時の筒内ガス密度ρg0を求める。
続いて、ＣＰＵ６１はステップ２０１０に進み、図６のステップ６１５と同様にして、ステップ５１０にて求めたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記１番目噴射時筒内容積Vg(CAinj(1))と、上記(4)式に相当する式とに基づいて１番目噴射時の筒内ガス圧力Pg0を求める。
次に、ＣＰＵ６１はステップ２０１５に進んで、ステップ５３０にて求めた基本燃料噴射圧力Pcrbaseから上記１番目噴射時の筒内ガス圧力Pg0を減じることで１番目噴射時の有効噴射圧力ΔPを求め、続くステップ６２５にて、上記求めた有効噴射圧力ΔPと、筒内ガス密度ρg0と、上記テーブルMapθとに基づいて１番目噴射に係わる噴霧角θを求める。これにより、噴霧角θは、１番目噴射時点（即ち、クランク角度CAinj(1)）での有効噴射圧力ΔP及び筒内ガス密度ρg0に基づいて決定される。
続いて、ＣＰＵ６１はステップ２０２０に進み、上述したように、計算の便宜上、１番目混合気に係わる噴射後経過時間ｔ＝０における１番目混合気の空気過剰率λ(0)（この値は、後述する図２１のステップ２１１０にて使用される。）を「０」に設定する。
次に、ＣＰＵ６１はステップ６３５、６４０の処理を順に実行し、その後、ステップ２０２５に進んで、図６のステップ６４５と同様にして、１番目混合気に係わる混合気のエンタルピHmixを、上記(17)式に相当する式に従って、初期値（即ち、先のステップ１９１０にて求めた１番目噴射量q(1)と、燃料の定圧比熱Cfと、ステップ５５０にて求めた燃料蒸気温度Tfの積）に設定する。
次いで、ＣＰＵ６１はステップ６５０の処理を実行し、その後、ステップ２０３０に進み、１番目混合気に係わる混合気質量Mmixを初期値である１番目噴射量q(1)に設定する。続いて、ＣＰＵ６１はステップ６６０の処理を実行し、その後、ステップ２０３５に進んで、図６のステップ６６５と同様、１番目混合気に係わる噴射後経過時間ｔを初期値「０」に設定する。また、ＣＰＵ６１は、１番目噴射に係わるクランク角度CAを初期値である１番目噴射時クランク角度CAinj(1)に設定する。これにより、１番目混合気に係わる噴射後経過時間ｔが１番目噴射時からカウントされることになる。加えて、このステップ２０３５では、ＣＰＵ６１は変数kの値を「０」に設定する。この変数kの値は、噴射後経過時間ｔが「k・Δｔ」であることを表す。
次に、ＣＰＵ６１はステップ６７０の処理を実行し、その後、図２１のルーチンに進んで、１番目噴射に関する混合気温度の算出のための処理を開始する。具体的には、ＣＰＵ６１は、先ず、ステップ２１０５に進み、図７のステップ７０５と同様に、１番目混合気に係わる噴射後経過時間ｔをΔｔだけ進めるとともに、１番目噴射に係わるクランク角度CAを図１９のステップ５６５で求めたΔCAだけ進める。
加えて、このステップ２１０５では、ＣＰＵ６１は変数kの値を「１」だけインクリメントする。これにより、１番目噴射に係わるクランク角度CAの値、及び変数kの値が、１番目混合気に係わる噴射後経過時間ｔに対応する値に維持されていく。
次に、ＣＰＵ６１はステップ７１０〜７３０の処理を順に実行し、その後、ステップ２１１０に進み、図７のステップ７３５と同様、１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおける１番目混合気の空気過剰率λ(k)（第１番目部分混合気指標値）を、上記(2)式、(3)式に従って、１番目混合気に係わる噴射後経過時間ｔ＝(k−1)・Δｔにおける１番目混合気の空気過剰率λ(k−1)（k＝1の場合、λ(k−1)＝λ(0)＝0）に、ステップ７３０にて求めた燃料希釈率dλ/dtにΔｔを乗じた値「dλ/dt・Δｔ」を加えることで求める。このステップ２１１０が混合指標値取得手段に相当する。
続いて、ＣＰＵ６１はステップ２１１５に進んで、ステップ２１１０にて既に求められているλ(k)，λ(k-1)と、図１９のステップ１９１０に求められているq(1)と、上記(29)式とに従って、１番目混合気に係わる噴射後経過時間ｔ＝(k−1)・Δｔ
〜 k・Δｔの間に１番目混合気に新たに取り込まれる筒内ガス量g(1,k)を求める。
次いで、ＣＰＵ６１はステップ２１２０に進み、図７のステップ７４５と同様、混合気形成筒内ガス質量Gを、その時点での値（初期値は、図２０のステップ６３５により「０」に設定されている。）に上記求めた１番目混合気に新たに取り込まれる筒内ガス質量g(1,k)を加えた値に更新する。これにより、噴射後経過時間ｔ＝k・Δｔにおける１番目混合気に係わる混合気形成筒内ガス質量Gが求められる。
次に、ＣＰＵ６１はステップ２１２５に進んで、図７のステップ７５０と同様、混合気質量Mmixを、その時点での値（初期値は、図２０のステップ２０３０により「q(1)」に設定されている。）に上記１番目混合気に新たに取り込まれる筒内ガス質量g(1,k)を加えた値に更新する。これにより、噴射後経過時間ｔ＝k・Δｔにおける１番目混合気に係わる混合気質量Mmixが求められる。
続いて、ＣＰＵ６１はステップ２１３０に進み、図７のステップ７５５と同様、化学反応前の混合気のエンタルピHpreを、その時点での混合気のエンタルピHmix（初期値は、図２０のステップ２０２５の処理により値「q(1)・Cf・Tf」）に、「上記１番目混合気に新たに取り込まれる筒内ガスの熱エネルギーΔHg＝g(1,k)・Cg・Tg」を加えた値に設定する。
次に、ＣＰＵ６１はステップ７６０〜７８０の処理を順に実行する。これにより、ステップ７８０において、噴射後経過時間ｔ＝k・Δｔ（従って、クランク角度CA＝CAinj(1)＋k・ΔCA）における１番目混合気の温度Tmixが求められる。
次いで、ＣＰＵ６１は図２２のルーチンに進み、１番目噴射に関する各種濃度の算出のための処理を開始する。この図２２のルーチンは、図８のルーチンのステップ８１５をステップ２２０５に代えた点のみにおいて図８のルーチンと異なる。このステップ２２０５では、ＣＰＵ６１は、図１９のステップ１９１０にて求めた１番目噴射量q(1)と、ステップ８１０にて求めた燃料消費量積算値sumqrと、図２１のステップ２１２５にて求めた混合気質量Mmixと、上記(19)式に相当する式とに基づいて、噴射後経過時間ｔ＝k・Δｔにおける１番目混合気内の燃料濃度[Fuel]mixを求める。
そして、ＣＰＵ６１は図２３のルーチンに進み、ステップ９０５の判定において「Ｎｏ」と判定する毎に、図２１のステップ２１０５〜図２３のステップ９０５の処理を繰り返し実行する。これにより、ステップ９０５の判定において「Ｎｏ」と判定される毎に、図２１のステップ２１０５にて、変数kの値が「１」だけインクリメントされ、１番目混合気に係わる噴射後経過時間ｔがΔｔだけ進められるとともに、１番目噴射に係わるクランク角度CAがΔCAだけ進められる。
即ち、ステップ９０５の判定において「Ｎｏ」と判定される限りにおいて、１番目混合気の空気過剰率λ(k)、１番目混合気に新たに取り込まれる筒内ガス質量g(1,k)、１番目混合気の状態（温度Tmix等）、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）１番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix）が１番目噴射時から微小時間Δｔ毎に（即ち、CAinj(1)から微小クランク角度ΔCA毎に）求められていく。
そして、ステップ９０５の条件が成立すると、ＣＰＵ６１は図２３のステップ９０５に進んだとき「Ｙｅｓ」と判定してステップ２３０５以降に進み、１番目混合気に係わる計算を終了するための処理を行う。
即ち、ＣＰＵ６１はステップ２３０５に進むと、図９のステップ９１０と同様に、１番目混合気内のスート発生量Soot(1)、及び１番目混合気内のＮＯ発生量NO(1)を求める。次いで、ＣＰＵ６１はステップ２３１０に進み、図９のステップ９１５と同様、スート総発生量Sootを、その時点での値（現時点では、図１９のステップ５７０の処理により初期値Soot0）に上記求めた１番目混合気内のスート発生量Soot(1)を加えた値（Soot0＋Soot(1)）に更新するとともに、ＮＯ総発生量ＮＯを、その時点での値（現時点では、図１９のステップ５７０の処理により初期値ＮＯ0）に上記求めた１番目混合気内のＮＯ発生量NO(1)を加えた値（ＮＯ0＋NO(1)）に更新する。以上により、１番目噴射（従って、１番目混合気）に関する計算が終了する。
次に、ＣＰＵ６１はステップ２３１５に進んで、変数iを「１」に設定する。ここで、変数iは、第１、第２実施形態と同様、何番目の噴射か（従って、何番目の混合気か）を識別するための値である。
続いて、ＣＰＵ６１はステップ２３２０に進み、変数iの値を「１」だけインクリメントし、続くステップ２３２５を経由して後述する図２４〜図２７に示した一連の「i(2≦i≦n)番目噴射に関するエミッション量の算出」ルーチンを実行することにより、i番目混合気内のエミッション発生量である、スート発生量Soot(i)、及びＮＯ発生量NO(i)を求める（2≦i≦n）。ここで、値nは、図１９のステップ１９０５にて求められた噴射燃料の区分数である。
次いで、ＣＰＵ６１はステップ２３３０に進んで、スート総発生量Sootを、その時点での値（i＝２の場合、Soot0＋Soot(1)）に上記求めたi番目混合気内のスート発生量Soot(i)を加えた値に更新するとともに、ＮＯ総発生量ＮＯを、その時点での値（i＝２の場合、ＮＯ0＋NO(1)）に上記求めたi番目混合気内のＮＯ発生量NO(i)を加えた値に更新する。これにより、スート総発生量Soot＝Soot0＋Soot(1)＋・・・＋Soot(i)となり、ＮＯ総発生量ＮＯ＝ＮＯ0＋NO(1)＋・・・＋NO(i)となる。
そして、ＣＰＵ６１はステップ２３３５に進んで、変数iの値が上記区分数nに一致したか否かを判定し、「Ｎｏ」と判定する場合、ステップ２３２０に戻る。即ち、ステップ２３２０の繰り返し実行により変数iが区分数nに達するまでの間、ステップ２３２０〜２３３０の処理が繰り返し実行される。この結果、変数iの値が「１」ずつインクリメントされていくとともに、ステップ２３３０にてスート総発生量Sootの値、及びＮＯ総発生量ＮＯが更新されていく。
そして、変数iの値が区分数nに達すると、ＣＰＵ６１はステップ２３３５に進んだとき「Ｙｅｓ」と判定してステップ９２５に進む。この時点で、エミッション発生量の推定が終了して、スート総発生量Sootが「Soot0＋Soot(1)＋・・・＋Soot(n)」に確定し、ＮＯ総発生量ＮＯが「ＮＯ0＋NO(1)＋・・・＋NO(n)」に確定する。ＣＰＵ６１は、以上説明した混合気状態の推定、及びエミッション発生量の推定をＩＶＣの直後（即ち、燃料噴射開始時期CAinjの前の時点）において完了する。
以下、図２４〜図２７に示した一連の「i(2≦i≦n)番目噴射に関するエミッション量の算出」ルーチンについて説明する。このルーチンは、i（2≦i≦n）番目混合気のスート発生量Soot(i)、及びＮＯ発生量NO(i)を求めるルーチンであって、図２４〜図２７のルーチンはそれぞれ、図６〜図９のルーチンに対応している。
ＣＰＵ６１は図２３のステップ２３２０を実行した後、ステップ２３２５を経由して、先ず、i(2≦i≦n)番目噴射に関する各種初期値を決定するために図２４のルーチンから実行する。具体的には、ＣＰＵ６１は、先ず、ステップ６３５、６４０の処理を順に実行し、その後、ステップ２４０５に進んで、図２０のステップ２０２５と同様にして、i番目混合気に係わる混合気のエンタルピHmixを、初期値（即ち、先のステップ１９１０にて求めたi番目噴射量q(i)と、燃料の定圧比熱Cfと、ステップ５５０にて求めた燃料蒸気温度Tfの積）に設定する。
次いで、ＣＰＵ６１はステップ６５０の処理を実行し、その後、ステップ２４１０に進み、i番目混合気に係わる混合気質量Mmixを初期値である上記i番目噴射量q(i)に設定する。続いて、ＣＰＵ６１はステップ６６０の処理を実行し、その後、ステップ２４１５に進んで、図２０のステップ２０３５と同様、i番目噴射に係わるクランク角度CAを初期値であるi番目噴射時クランク角度CAinj(i)に設定する。加えて、このステップ２４１５では、ＣＰＵ６１は変数kの値を「０」に設定する。この変数kの値は、i番目噴射からの経過時間が「k・Δｔ」であること（即ち、i番目噴射に係わるクランク角度CAが「CAinj(i)＋k・ΔCA」であること）を表す。
次に、ＣＰＵ６１はステップ６７０の処理を実行し、その後、図２５のルーチンに進んで、i（2≦i≦n）番目噴射に関する混合気温度の算出のための処理を開始する。具体的には、ＣＰＵ６１は、先ず、ステップ２５０５に進み、図２１のステップ２１０５と同様に、i番目噴射に係わるクランク角度CAを図１９のステップ５６５で求めたΔCAだけ進める。加えて、このステップ２５０５では、ＣＰＵ６１は変数kの値を「１」だけインクリメントする。これにより、i番目噴射に係わるクランク角度CAの値と変数kの値とが、互いに対応する値に維持されていく。
次に、ＣＰＵ６１はステップ７１５〜７２５の処理を順に実行し、その後、ステップ２５１０に進み、図２１のステップ２１１５にて求めた「１番目混合気に新たに取り込まれる筒内ガス量g(1,k)」と、図１９のステップ１９１０にて求めたq(i)と、同ステップ１９１０にて求めたq(1)と、上記(34)式とに基づいて、上記(2)式、(3)式を利用することなく、i（2≦i≦n）番目混合気に係わる噴射後経過時間ｔ＝(k−1)・Δｔ
〜 k・Δｔの間にi番目混合気に新たに取り込まれる筒内ガス量g(i,k)を求める。
次いで、ＣＰＵ６１はステップ２５１５に進み、図２１のステップ２１２０と同様、混合気形成筒内ガス質量Gを、その時点での値（初期値は、図２４のステップ６３５により「０」に設定されている。）に上記求めたi番目混合気に新たに取り込まれる筒内ガス質量g(i,k)を加えた値に更新する。これにより、噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気に係わる混合気形成筒内ガス質量Gが求められる。
次に、ＣＰＵ６１はステップ２５２０に進んで、図２１のステップ２１２５と同様、混合気質量Mmixを、その時点での値（初期値は、図２４のステップ２４１０により「q(i)」に設定されている。）に上記i番目混合気に新たに取り込まれる筒内ガス質量g(i,k)を加えた値に更新する。これにより、噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気に係わる混合気質量Mmixが求められる。
続いて、ＣＰＵ６１はステップ２５２５に進み、図２１のステップ２１３０と同様、化学反応前の混合気のエンタルピHpreを、その時点での混合気のエンタルピHmix（初期値は、図２４のステップ２４０５の処理により値「q(i)・Cf・Tf」）に、「上記i番目混合気に新たに取り込まれる筒内ガスの熱エネルギーΔHg＝g(i,k)・Cg・Tg」を加えた値に設定する。
次に、ＣＰＵ６１はステップ７６０〜７８０の処理を順に実行する。これにより、ステップ７８０において、噴射後経過時間ｔ＝k・Δｔ（従って、クランク角度CA＝CAinj(i)＋k・ΔCA）におけるi（2≦i≦n）番目混合気の温度Tmixが求められる。
次いで、ＣＰＵ６１は図２６のルーチンに進み、i番目噴射に関する各種濃度の算出のための処理を開始する。この図２６のルーチンは、図８のルーチンのステップ８１５をステップ２６０５に代えた点のみにおいて図８のルーチンと異なる。このステップ２６０５では、ＣＰＵ６１は、図１９のステップ１９１０にて求めたi番目噴射量q(i)と、ステップ８１０にて求めた燃料消費量積算値sumqrと、図２５のステップ２５２０にて求めた混合気質量Mmixと、上記(19)式に相当する式とに基づいて、噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気内の燃料濃度[Fuel]mixを求める。
そして、ＣＰＵ６１は図２７のルーチンに進み、ステップ９０５の判定において「Ｎｏ」と判定する毎に、図２５のステップ２５０５〜図２７のステップ９０５の処理を繰り返し実行する。これにより、ステップ９０５の判定において「Ｎｏ」と判定される毎に、図２５のステップ２５０５にて、変数kの値が「１」だけインクリメントされ、i番目噴射に係わるクランク角度CAがΔCAだけ進められる。
即ち、ステップ９０５の判定において「Ｎｏ」と判定される限りにおいて、i番目混合気に新たに取り込まれる筒内ガス質量g(i,k)、i番目混合気の状態（温度Tmix等）、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）i番目混合気に係わるエミッション濃度（[Soot]mix，[NO]mix）がi番目噴射時から微小時間Δｔ毎に（即ち、CAinj(i)から微小クランク角度ΔCA毎に）求められていく。
そして、ステップ９０５の条件が成立すると、ＣＰＵ６１は図２７のステップ９０５に進んだとき「Ｙｅｓ」と判定してステップ９１０に進んで、i番目混合気内のスート発生量Soot(i)、及びi番目混合気内のＮＯ発生量NO(i)を求める。i番目混合気内のスート発生量Soot(i)は、図２６のステップ８５５の処理にて更新されている現時点でのi番目混合気内のスート濃度[Soot]mixの値に、図２５のステップ２５２０にて更新されている現時点での混合気質量Mmixを乗じることで求められる。また、i番目混合気内のＮＯ発生量NO(i)は、図２６のステップ８８０の処理にて更新されている現時点でのi番目混合気内のＮＯ濃度[NO]mixの値に、上記混合気質量Mmixを乗じることで求められる。
そして、ＣＰＵ６１はステップ２４９５を経由して（即ち、図２４〜図２７に示した一連のルーチンの実行を終了して）上述した図２３のステップ２３３０に進む。このようにして、上述した図２３のステップ２３２０〜２３３０の処理が繰り返し実行される毎に、変数i（2≦i≦n）の値が「１」ずつインクリメントされていくとともに、図２７のステップ９１０にてi番目混合気内のスート発生量Soot(i)、及びi番目混合気内のＮＯ発生量NO(i)が順次求められていく。以上、図２４〜図２７に示した一連のルーチンについて説明した。
上述したように、変数iの値が区分数nに達すると、ＣＰＵ６１は図２３のステップ２３３５に進んだとき「Ｙｅｓ」と判定してステップ９２５〜９６０を実行し、図２３のステップ２３３０にて求められているスート総発生量Sootの値、及びＮＯ総発生量ＮＯの値（従って、エミッション総発生量の推定結果）に基づいて噴射圧力を補正し、ステップ１９９５に進んで図１９〜図２３の一連の本ルーチンを一旦終了する。以降、ＣＰＵ６１は、次のＩＶＣが到来するまでの間、図１９のステップ５０５に進む毎に「Ｎｏ」と判定し続ける。
以上、説明したように、本発明による混合気状態推定装置、及びエミッション発生量推定装置の第３実施形態によれば、噴射期間TAUが多数のn（＝TAU/Δｔ）個の期間に区分され、i（1≦i≦n）番目噴射に係わる噴射期間「i番目の微小時間Δｔ」分の質量q(i)の燃料が燃料噴射開始時期CAinjから(i−1)・Δｔ経過した後に一時に噴射されるものとする。１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおける１番目混合気の空気過剰率λ(k)（第１番目部分混合指標値）は、上記第１、第２実施形態と同様、上記(2)式、(3)式を利用して求められ、１番目混合気の状態（温度Tmix等）は係るλ(k)を利用して推定される。
そして、２番目以降のi番目混合気（i：２以上n以下の自然数）に係わる噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気の空気過剰率は、上記「１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおける１番目混合気の空気過剰率λ(k)」と等しいとの仮定のもと、上記(2)式、(3)式を利用することなく、２番目以降のi番目混合気の状態（温度Tmix等）が推定される。
従って、上記(2)式、(3)式を利用した空気過剰率λの算出を必要とする混合気を１番目混合気のみとすることができる。この結果、計算負荷の大きい上記(2)式、(3)式を利用した空気過剰率λの算出回数を減らすことができるから、ＣＰＵ６１の計算負荷を減らすことができる。
本発明は上記第３実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記第３実施形態では、２番目以降のi番目混合気（i：２以上n以下の自然数）に係わる噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気の空気過剰率は、上記「１番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおける１番目混合気の空気過剰率λ(k)」と等しいと仮定されているが、同空気過剰率λ(k)に所定の係数を乗じた値になると仮定してもよい。
即ち、２番目以降のi番目混合気に係わる噴射後経過時間ｔ＝k・Δｔにおけるi番目混合気の空気過剰率をλ(i,k)（i:２以上n以下の自然数）とすると、例えば、「λ(i,k)＝ｈ(i)・λ(k)」と設定される。ここで、ｈ(i)は変数iに応じて決定される係数である。
この場合、i（2≦i≦n）番目混合気に係わる噴射後経過時間ｔ＝(k−1)・Δｔ 〜
k・Δｔの間にi番目混合気に新たに取り込まれる筒内ガス量g(i,k)は、図２５のステップ２５１０に記載された上記(34)式ではなく、上記(33)式に相当する式である「g(i,k)＝stoich・(λ(i,k)−λ(i,k−1))・q(i)」に従って取得され得る。
（第４実施形態）
次に、本発明の第４実施形態に係る内燃機関のエミッション発生量推定装置について説明する。この第４実施形態は、混合気の着火後において燃焼室内で所謂定常火炎が発生している場合を想定し、係る定常火炎に特有の特性を利用して混合気の状態（温度等）、及びエミッション発生量を推定する点で、上記第１〜第３実施形態と異なっている。以下、係る相違点を中心に説明する。
噴射開始時点以降における噴射された燃料（従って、混合気）の燃料噴射弁２１の噴孔からの到達距離（以下、「混合気到達距離X」と称呼する。）は、例えば、上記非特許文献１にて紹介された実験式である下記(35)式、(36)式に従って噴射後経過時間ｔの関数として求めることができる。下記(36)式において、dX/dtは噴射後経過時間ｔの関数である混合気移動速度である。なお、下記(36)式の右辺に示される各種値は、上記(3)式の右辺に示されるものと同一である。
Figure 0004424242
Figure 0004424242
ここで、上記(3)式の両辺を上記(36)式の両辺で除することにより、dλ/dXを、ρgとtanθと各種定数とを用いて表すことができる。いま、ρgとtanθが噴射後において一定であるとすると、dλ/dXは或る定数（正の値）になる。更に、空気過剰率λは混合気到達距離X＝０（従って、噴射時点）において「０」である。
従って、この場合、混合気到達距離Xと空気過剰率λとの関係は図２８に示すような線形関係となる。即ち、混合気到達距離Xが「０」から増加するにつれて空気過剰率λが「０」から増加し、X＝X0のとき、λ＝１となる。
ところで、燃料の噴射期間TAUが比較的長いと、図２８に示すように、混合気の着火後において燃焼室内にて所謂定常火炎（或いは、定常火炎に非常に近い形態を有する火炎）が発生する場合がある。
定常火炎が発生している場合、定常火炎内における燃料が過剰な領域（即ち、空気過剰率λ＜１となる領域。以下、「リッチ領域」と称呼する。）では、混合気内に取り込まれた筒内ガス中の酸素が不足していることにより酸素が全て燃焼により消費されている。即ち、定常火炎のリッチ領域内（図２８では、0≦X≦X0に対応する。）では、定常状態における酸素濃度（以下、「定常酸素濃度[O2]mixsteady」と称呼する。）が「０」になっている。なお、定常火炎のリッチ領域内では、酸素が全て燃焼により消費されても燃料が残存しているから、定常状態における燃料濃度（以下、「定常燃料濃度[Fuel]mixsteady」と称呼する。）は「０」より大きい値になっている。
一方、定常火炎が発生している場合、定常火炎内における酸素が過剰な領域（即ち、空気過剰率λ＞１となる領域。以下、「リーン領域」と称呼する。）では、燃料が不足していることにより燃料が全て燃焼により消費されている。即ち、定常火炎のリーン領域内（図２８では、X0≦Xに対応する。）では、定常燃料濃度[Fuel]mixsteadyが「０」になっている。なお、定常火炎のリーン領域内では、燃料が全て燃焼により消費されても酸素が残存しているから、定常酸素濃度[O2]mixsteadyは「０」より大きい値になっている。
他方、上述したように、混合気内のスート発生速度d[Soot]mix/dtを求める上記(23)式には、上記(24)式に示すように混合気内の燃料濃度[Fuel]mixの値を因数に含むスートの生成速度dmsf/dtを求める項と、上記(25)式に示すように混合気内の酸素濃度[O2]mixの値を因数に含むスートの酸化速度dmso/dtを求める項とが存在する。
従って、定常火炎のリッチ領域内におけるスート発生速度d[Soot]mix/dtを上記(23)式に従って計算する場合、上記(25)式に示すスートの酸化速度dmso/dtを求める項の値が常に「０」に維持される。換言すれば、この場合、上記スートの酸化速度dmso/dtを求める項の計算を省略することができる。
同様に、定常火炎のリーン領域内におけるスート発生速度d[Soot]mix/dtを上記(23)式に従って計算する場合、上記(24)式に示すスートの生成速度dmsf/dtを求める項の値が常に「０」に維持される。換言すれば、この場合、上記スートの生成速度dmsf/dtを求める項の計算を省略することができる。
加えて、上記(24)式、(25)式の右辺には、時々刻々と変化する変数Pg，Tmixに関するべき演算が含まれている。べき演算には多大な計算負荷が伴う。従って、(24)式、(25)式の算出回数をできるだけ減らすことが望ましい。
そこで、定常火炎が占める領域内（即ち、混合気内）におけるスート発生速度d[Soot]mix/dtを上記(23)式に基づいて推定する第４実施形態は、定常火炎のリッチ領域内（λ＜１）では、上記(23)式においてdmso/dtの項を省略することで得られる式である「d[Soot]mix/dt＝dmsf/dt」に従ってスート発生速度d[Soot]mix/dtを求める。
同様に、第４実施形態は、定常火炎のリーン領域内（λ＞１。実際には、λ≧１）では、上記(23)式においてdmsf/dtの項を省略することで得られる式である「d[Soot]mix/dt＝−dmso/dt」に従ってスート発生速度d[Soot]mix/dtを求める。これらにより、(24)式、(25)式の算出回数を減らすことができ、この結果、ＣＰＵ６１の計算負荷を減らすことができる。
以下、定常火炎が占める領域内（即ち、混合気内）におけるスート発生速度d[Soot]mix/dtを上記(23)式を利用して求める際に必要となる、混合気温度（具体的には、定常状態における燃焼による混合気温度上昇量。以下、「定常混合気温度上昇量ΔTmixsteady」と称呼する。）、上記定常燃料濃度[Fuel]mixsteady、及び、上記定常酸素濃度[O2]mixsteady、並びに、定常火炎が占める領域内（即ち、混合気内）におけるＮＯ発生速度d[NO]mix/dtを上記(27)式を利用して求める際に必要となる、定常状態における窒素濃度（以下、「定常窒素濃度[N2]mixsteady」と称呼する。）を求める手法について順に説明する。
以下、係る説明に必要となる空気過剰率λを、上記筒内ガス理論空燃比stoichと、燃料量（質量）Qと、混合気形成筒内ガス質量Gとを用いて下記(37)式に従って定義する。
Figure 0004424242
<定常混合気温度上昇量ΔTmixsteady>
定常混合気温度上昇量ΔTmixsteadyは、下記(38)式にて表すことができる。下記(38)式において、Qreacは混合気内で燃焼により消費される燃料量（質量）であり、Cgは筒内ガスの定圧比熱であり、Cfは燃料の定圧比熱である。HfはQreacを反応熱Hr(＝Hf・Qreac)に変換するための係数である。この(38)式は、「混合気内の質量Qreacの燃料の燃焼により発生した反応熱Hrは、質量Qの燃料と質量Gの筒内ガスとを含む混合気の温度を定常混合気温度上昇量ΔTmixsteadyだけ上昇させるために消費される」との仮定に基づくものである。
Figure 0004424242
ここで、定常火炎のリッチ領域内（λ＜１）では、混合気内で燃焼により消費される燃料量Qreacは、混合気内における筒内ガス中の全ての酸素と反応して消費される分の燃料量となるから、下記(39)式にて求めることができる。
Figure 0004424242
以上より、上記(37)式〜(39)式からQ,G,Qreacを消去して整理すると、定常火炎のリッチ領域内（λ＜１）における定常混合気温度上昇量ΔTmixsteadyは、空気過剰率λの関数として下記(40)式にて表すことができる。
Figure 0004424242
他方、定常火炎のリーン領域内（λ＞１。実際には、λ≧１）では、混合気内の燃料が燃焼により全て消費されるから、混合気内で燃焼により消費される燃料量Qreacは「Q」と等しくなる。以上より、「Qreac＝Q」の関係と、上記(37)式及び(38)式とからQ,G,Qreacを消去して整理すると、定常火炎のリーン領域内（λ≧１）における定常混合気温度上昇量ΔTmixsteadyは、空気過剰率λの関数として下記(41)式にて表すことができる。
Figure 0004424242
図２９は、上記(40)式、及び(41)式にて表される空気過剰率λと定常混合気温度上昇量ΔTmixsteadyとの関係を示したグラフである。なお、混合気到達距離Xが図２８に示すように空気過剰率λと線形関係にあることを考慮すると、図２９は、混合気到達距離Xと定常混合気温度上昇量ΔTmixsteadyとの関係を示したグラフであるということもできる。図２９に示すように、定常混合気温度上昇量ΔTmixsteadyは、λ＝１（即ち、X＝X0。図２８を参照。）にて最大値ΔTmix0＝Hf/(stoich・Cg＋Cf)を採る。
<定常燃料濃度[Fuel]mixsteady>
上述したように、定常火炎のリーン領域内では、定常燃料濃度[Fuel]mixsteadyは、空気過剰率λ（従って、混合気到達距離X）にかかわらず「０」に維持される。一方、定常火炎のリッチ領域内（λ＜１）では、混合気内に残存している燃料量（質量）は、上記(39)式を考慮すると、(Q−(G/stoich))と表すことができる。従って、混合気質量(Q＋G)に対する「混合気内に残存している燃料量（質量）」の割合である定常燃料濃度[Fuel]mixsteadyは、下記(42)式にて表すことができる。
Figure 0004424242
従って、上記(37)式、及び(42)式からQ,Gを消去して整理すると、定常火炎のリッチ領域内（λ＜１）における定常燃料濃度[Fuel]mixsteadyは、空気過剰率λの関数として下記(43)式にて表すことができる。
Figure 0004424242
図３０は、上記(43)式にて表される空気過剰率λと定常燃料濃度[Fuel]mixsteadyとの関係を示したグラフである。なお、混合気到達距離Xが空気過剰率λと線形関係にあることを考慮すると、図３０は、混合気到達距離Xと定常燃料濃度[Fuel]mixsteadyとの関係を示したグラフであるということもできる。
<定常酸素濃度[O2]mixsteady>
上述したように、定常火炎のリッチ領域内では、定常酸素濃度[O2]mixsteadyは、空気過剰率λ（従って、混合気到達距離X）にかかわらず「０」に維持される。一方、定常火炎のリーン領域内（λ＞１。実際には、λ≧１）では、混合気内に残存している酸素量（質量）は、上記「Qreac＝Q」の関係式を考慮すると、(G−(Q・stoich))・[O2]in
と表すことができる。ここで、[02]inは吸気酸素濃度（従って、筒内ガス内の酸素濃度）である。従って、混合気質量(Q＋G)に対する「混合気内に残存している酸素量（質量）」の割合である定常酸素濃度[02]mixsteadyは、下記(44)式にて表すことができる。
Figure 0004424242
従って、上記(37)式、及び(44)式からQ,Gを消去して整理すると、定常火炎のリーン領域内（λ≧１）における定常酸素濃度[O2]mixsteadyは、空気過剰率λの関数として下記(45)式にて表すことができる。
Figure 0004424242
図３１は、上記(45)式にて表される空気過剰率λと定常酸素濃度[O2]mixsteadyとの関係を示したグラフである。なお、混合気到達距離Xが空気過剰率λと線形関係にあることを考慮すると、図３１は、混合気到達距離Xと定常酸素濃度[O2]mixsteadyとの関係を示したグラフであるということもできる。
<定常窒素濃度[N2]mixsteady>
上述したように、筒内ガス中の窒素は不活性ガスであるから混合気内で化学反応により消費されることがない。従って、混合気内に残存している窒素量（質量）は、G・[N2]in と表すことができる。ここで、[N2]inは吸気窒素濃度（従って、筒内ガス内の窒素濃度）である。従って、混合気質量(Q＋G)に対する「混合気内に残存している窒素量（質量）」の割合である定常窒素濃度[N2]mixsteadyは、下記(46)式にて表すことができる。
Figure 0004424242
従って、上記(37)式、及び(46)式からQ,Gを消去して整理すると、定常窒素濃度[N2]mixsteadyは、空気過剰率λの関数として下記(47)式にて表すことができる。
Figure 0004424242
図３２は、上記(47)式にて表される空気過剰率λと定常窒素濃度[N2]mixsteadyとの関係を示したグラフである。なお、混合気到達距離Xが空気過剰率λと線形関係にあることを考慮すると、図３２は、混合気到達距離Xと定常窒素濃度[N2]mixsteadyとの関係を示したグラフであるということもできる。
以上説明したように、定常火炎が占める領域内（即ち、混合気内）においては、定常混合気温度上昇量ΔTmixsteady、定常燃料濃度[Fuel]mixsteady、定常酸素濃度[O2]mixsteady、及び定常窒素濃度[N2]mixsteadyを、空気過剰率λの関数でそれぞれ表すことができる。
他方、空気過剰率λは、上記(2)式、(3)式により噴射後経過時間ｔの関数funcλ(t)として表すことができる。図３３は、関数funcλ(t)から得られる噴射後経過時間ｔと空気過剰率λとの関係を示したグラフである。なお、t＝t0のとき、λ＝１となるものとする。
以上のことから、図２９〜図３２に示したそれぞれの関係に対して上記関数funcλ(t)から得られる図３３に示した関係を適用することにより、定常火炎が占める領域内（即ち、混合気内）においては、定常混合気温度上昇量ΔTmixsteady、定常燃料濃度[Fuel]mixsteady、定常酸素濃度[O2]mixsteady、及び定常窒素濃度[N2]mixsteadyをそれぞれ、噴射後経過時間ｔの関数である、funcΔTmixsteady(t)，func[Fuel]mixsteady(t)，func[02]mixsteady(t)，func[N2]mixsteady(t)で表すことができる。
図３４〜図３７は、funcΔTmixsteady(t)，func[Fuel]mixsteady(t)，func[02]mixsteady(t)，func[N2]mixsteady(t)で表される関係をそれぞれ示したグラフである。第４実施形態は、係る噴射後経過時間ｔの関数であるfuncΔTmixsteady(t)，func[Fuel]mixsteady(t)，func[02]mixsteady(t)，func[N2]mixsteady(t)を利用して、定常火炎が占める領域内（即ち、混合気内）におけるスート発生速度d[Soot]mix/dt、及びＮＯ発生速度d[NO]mix/dtを求める。
（第４実施形態の実際の作動）
以下、第４実施形態に係るエミッション発生量推定装置の実際の作動について説明する。この装置のＣＰＵ６１は、第１実施形態のＣＰＵ６１が実行する図５〜図９に示した一連のルーチン、及び図１３に示したルーチンのうち、図１３に示したルーチンのみをそのまま実行するとともに、図５〜図９に示した一連のルーチンに代えて図３８〜図４１にフローチャートにより示した一連のルーチンを実行する。
なお、図３８〜図４１に示したルーチンのステップであって図５〜図９に示したルーチンのステップと同一のものについては、図５〜図９に示したルーチンの対応するステップの番号と同一の番号を付すことでそれらの説明を省略する。以下、第４実施形態に特有の図３８〜図４１に示した各ルーチンについて説明する。なお、図３８のルーチンは図５のルーチンに、図３９のルーチンは図６のルーチンに、図４０のルーチンは図７、及び図８のルーチンに、図４１のルーチンは図９のルーチンに対応している。
第４実施形態のＣＰＵ６１は、図５〜図９に示した一連のルーチンに対応する図３８〜図４１に示した一連のルーチンを所定時間の経過毎に繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１は図３８のステップ３８００から処理を開始し、ステップ５０５に進んで、「Ｙｅｓ」と判定する場合（即ち、ＩＶＣが到来した場合）、ステップ５１０〜５３０の処理を順に実行し、その後、ステップ５５０〜５６０の処理を順に実行する。
次に、ＣＰＵ６１は、図３９のルーチンに進んで各種関数の決定を行うための処理を行う。具体的には、ＣＰＵ６１はステップ３９０５に進んで、定常火炎が発生している間の筒内ガス密度ρgの平均値である定常火炎時平均筒内ガス密度ρgaveを求める。定常火炎時平均筒内ガス密度ρgaveは、筒内ガスの全質量Mgに依存して決定されるから、図３８のステップ５１５にて求めた筒内ガスの全質量Mgと、Mgを引数とするρgaveを求める関数funcρgaveとから求められる。
続いて、ＣＰＵ６１はステップ３９１０に進み、定常火炎が発生している間の筒内ガス圧力Pgの平均値である定常火炎時平均筒内ガス圧力Pgaveを求める。定常火炎時平均筒内ガス圧力Pgaveは、ＩＶＣ時筒内ガス圧力Pgivc、ＩＶＣ時クランク角度CAivc、指令燃料噴射量Qfinに依存して決定されるから、図３８のステップ５１０にて求めたＩＶＣ時筒内ガス圧力Pgivc及びＩＶＣ時クランク角度CAivcと、図３８のステップ５２０にて求めた指令燃料噴射量Qfinと、Pgivc,CAivc,Qfinを引数とするPgaveを求める関数funcPgaveとから求められる。
次いで、ＣＰＵ６１はステップ３９１５に進んで、図３８のステップ５３０にて求めた基本燃料噴射圧力Pcrbaseから上記定常火炎時平均筒内ガス圧力Pgaveを減じることで定常火炎時平均有効噴射圧力ΔPaveを求め、続くステップ３９２０にて、上記求めた定常火炎時平均有効噴射圧力ΔPaveと、上記求めた定常火炎時平均筒内ガス密度ρgaveと、上記テーブルMapθとに基づいて定常火炎時平均噴霧角θaveを求める。これにより、定常火炎時平均噴霧角θaveは、ΔPave及びρgaveに基づいて決定される。
続いて、ＣＰＵ６１はステップ３９２５に進み、上記(3)式の有効噴射圧力ΔP、筒内ガス密度ρg、及び噴霧角θとして、上記求めた定常火炎時平均有効噴射圧力ΔPave、定常火炎時平均筒内ガス密度ρgave、及び定常火炎時平均噴霧角θaveをそれぞれ使用することで、上記(2)式、(3)式を利用して、噴射後経過時間ｔと空気過剰率λとの関係を規定する上記関数funcλ(t)（図３３を参照）を決定する。
次に、ＣＰＵ６１はステップ３９３０に進んで、上記求めた関数funcλ(t)と、上記(40)式、及び(41)式にて表される空気過剰率λを引数とするΔTmixsteadyを求めるための関数funcΔTmixsteady(λ)とから、噴射後経過時間ｔを引数とする定常混合気温度上昇量ΔTmixsteadyを求めるための関数funcΔTmixsteady(t)を決定する。
続いて、ＣＰＵ６１はステップ３９３５に進んで、上記求めた関数funcλ(t)と、上記(43)式にて表される空気過剰率λを引数とする[Fuel]mixsteadyを求めるための関数func[Fuel]mixsteady(λ)とから、噴射後経過時間ｔを引数とする定常燃料濃度[Fuel]mixsteadyを求めるための関数func[Fuel]mixsteady(t)を決定する。
次いで、ＣＰＵ６１はステップ３９４０に進んで、上記求めた関数funcλ(t)と、上記(45)式にて表される空気過剰率λを引数とする[O2]mixsteadyを求めるための関数func[O2]mixsteady(λ)とから、噴射後経過時間ｔを引数とする定常酸素濃度[O2]mixsteadyを求めるための関数func[O2]mixsteady(t)を決定する。
次に、ＣＰＵ６１はステップ３９４５に進んで、上記求めた関数funcλ(t)と、上記(47)式にて表される空気過剰率λを引数とする[N2]mixsteadyを求めるための関数func[N2]mixsteady(λ)とから、噴射後経過時間ｔを引数とする定常窒素濃度[N2]mixsteadyを求めるための関数func[N2]mixsteady(t)を決定する。
続いて、ＣＰＵ６１はステップ３９５０に進んで、図６のステップ６６０と同様、混合気内のＮＯ濃度［NO］mix、及びスート濃度[Soot]mixをそれぞれ初期値「０」に設定し、続くステップ３９５５にて、図６のステップ６６５と同様、噴射後経過時間ｔを初期値「０」に設定する。
次に、ＣＰＵ６１はステップ６７０の処理を実行し、その後、図４０のルーチンに進み、定常火炎が占める領域内（従って、混合気内）における混合気温度、各種濃度等の算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ４００２に進み、噴射後経過時間ｔの値（初期値は、図３９のステップ３９５５の処理により「０」）を微小時間Δｔ（例えば、0.1msec）だけ進め、続くステップ４００４にて、現時点での噴射後経過時間ｔの値と、図３９のステップ３９２５にて決定した関数funcλ(t)とに基づいて噴射後経過時間ｔに対応する空気過剰率λを求める。
次いで、ＣＰＵ６１はステップ４００６に進み、現時点での噴射後経過時間ｔの値と、図３９のステップ３９３０にて決定した関数funcΔTmixsteady(t)とに基づいて噴射後経過時間ｔに対応する定常混合気温度上昇量ΔTmixsteadyを求め、同求めたΔTmixsteadyに図３８のステップ５５０にて求めた燃料蒸気温度Tfを加えて、噴射後経過時間ｔに対応する定常混合気温度Tmixsteadyを求める。
続いて、ＣＰＵ６１はステップ４００８に進んで、現時点での噴射後経過時間ｔの値と、図３９のステップ３９３５にて決定した関数func[Fuel]mixsteady(t)とに基づいて噴射後経過時間ｔに対応する定常燃料濃度[Fuel]mixsteadyを求める。
次に、ＣＰＵ６１はステップ４０１０に進んで、現時点での噴射後経過時間ｔの値と、図３９のステップ３９４０にて決定した関数func[O2]mixsteady(t)とに基づいて噴射後経過時間ｔに対応する定常酸素濃度[O2]mixsteadyを求める。
次いで、ＣＰＵ６１はステップ４０１２に進んで、現時点での噴射後経過時間ｔの値と、図３９のステップ３９４５にて決定した関数func[N2]mixsteady(t)とに基づいて噴射後経過時間ｔに対応する定常窒素濃度[N2]mixsteadyを求める。
次に、ＣＰＵ６１はステップ８３５の判定を行う。いま、「Ｙｅｓ」と判定するものとすると、ＣＰＵ６１はステップ４０１６に進み、ステップ４００４にて更新された現時点での噴射後経過時間ｔに対応する空気過剰率λが「１」より小さいか否か（即ち、リッチ領域か否か）を判定する。
ここで、「Ｙｅｓ」と判定する場合（即ち、リッチ領域である場合）、ＣＰＵ６１はステップ４０１８に進んで、ステップ４００８にて求めた現時点での噴射後経過時間ｔに対応する定常燃料濃度[Fuel]mixsteadyと、図３９のステップ３９１０にて求めた定常火炎時平均筒内ガス圧力Pgaveと、ステップ４００６にて求めた現時点での噴射後経過時間ｔに対応する定常混合気温度Tmixsteadyと、上記(24)式とに基づいてスートの生成速度dmsf/dtを求める。続いて、ＣＰＵ６１はステップ４０２０に進み、上記(23)式から「dmso/dt」の項を削除したステップ４０２０内に記載の式と、上記求めたスートの生成速度dmsf/dtとから、現時点での噴射後経過時間ｔに対応するスート発生速度d[soot]mix/dtを求める。即ち、この場合、上記(25)式の計算が省略される。
一方、ステップ４０１６にて「Ｎｏ」と判定する場合（即ち、リーン領域である場合）、ＣＰＵ６１はステップ４０２２に進んで、その時点での混合気内のスート濃度[Soot]mix（初期値は、図３９のステップ３９５０の処理により「０」）と、ステップ４０１０にて求めた現時点での噴射後経過時間ｔに対応する定常酸素濃度[O2]mixsteadyと、図３９のステップ３９１０にて求めた定常火炎時平均筒内ガス圧力Pgaveと、ステップ４００６にて求めた現時点での噴射後経過時間ｔに対応する定常混合気温度Tmixsteadyと、上記(25)式とに基づいてスートの酸化速度dmso/dtを求める。続いて、ＣＰＵ６１はステップ４０２４に進み、上記(23)式から「dmsf/dt」の項を削除したステップ４０２４内に記載の式と、上記求めたスートの酸化速度dmso/dtとから、現時点での噴射後経過時間ｔに対応するスート発生速度d[soot]mix/dtを求める。即ち、この場合、上記(24)式の計算が省略される。
上記ステップ４０２０、或いはステップ４０２４にて現時点での噴射後経過時間ｔに対応するスート発生速度d[soot]mix/dtが求められると、ＣＰＵ６１はステップ８５５の処理を実行し、その後、ステップ４０２８に進んで、現時点での噴射後経過時間ｔに対応する空気過剰率λが「１」より大きく、且つ、上記求めた現時点での噴射後経過時間ｔに対応する定常混合気温度Tmixsteadyが上記スートの反応限界温度TminSootより低いか否かを判定する。
そして、図８のステップ８６０の判定と同様、ＣＰＵ６１は、ステップ４０２８にて「Ｙｅｓ」と判定する場合、ステップ８６５の処理を実行してからステップ８７０に進み、ステップ４０２８にて「Ｎｏ」と判定する場合、ステップ８７０に直ちに進む。
ＣＰＵ６１は、ステップ８７０に進むと、同ステップ８７０の判定を行う。いま、「Ｙｅｓ」と判定するものとすると、ＣＰＵ６１はステップ４０３４に進み、ステップ４００６にて求めた現時点での噴射後経過時間ｔに対応する定常混合気温度Tmixsteadyと、ステップ４０１０にて求めた現時点での噴射後経過時間ｔに対応する定常酸素濃度[O2]mixsteadyと、ステップ４０１２にて求めた現時点での噴射後経過時間ｔに対応する定常窒素濃度[N2]mixsteadyと、上記(27)式とに基づいて現時点での噴射後経過時間ｔに対応するＮＯ発生速度d[NO]mix/dtを求める。
続いて、ＣＰＵ６１はステップ８８０の処理を実行し、その後、ステップ４０３８に進んで、現時点での噴射後経過時間ｔに対応する空気過剰率λが「１」より大きく、且つ、上記求めた現時点での噴射後経過時間ｔに対応する定常混合気温度Tmixsteadyが上記ＮＯの反応限界温度TminNOより低いか否かを判定する。
そして、図８のステップ８８５の判定と同様、ＣＰＵ６１は、ステップ４０３８にて「Ｙｅｓ」と判定する場合、ステップ８９０の処理を実行してから図４１のルーチンに進み、ステップ４０３８にて「Ｎｏ」と判定する場合、図４１のルーチンに直ちに進む。
ＣＰＵ６１は図４１のルーチンに進むと、図９のステップ９０５に対応するステップ４１０５の判定を行う。ＣＰＵ６１は、ステップ４１０５の判定において「Ｎｏ」と判定する毎に、図４０のステップ４００２〜図４１のステップ４１０５の処理を繰り返し実行する。これにより、ステップ４１０５の判定において「Ｎｏ」と判定される毎に、図４０のステップ４００２にて、噴射後経過時間ｔがΔｔだけ進められる。
即ち、ステップ４１０５の判定において「Ｎｏ」と判定される限りにおいて、空気過剰率λ、及び（フラグＥＮＤsoot＝フラグＥＮＤno＝０である限りにおいて）エミッション濃度（[Soot]mix，[NO]mix）が、噴射後経過時間ｔの微小時間Δｔ毎の進行に対応して求められていく。
そして、ステップ４１０５の条件が成立すると、ＣＰＵ６１は図４１のステップ４１０５に進んだとき「Ｙｅｓ」と判定してステップ４１１０以降に進み、エミッション発生量に係わる計算を終了するための処理を行う。
即ち、ＣＰＵ６１はステップ４１１０に進むと、現時点（即ち、ステップ４１０５の条件が成立した時点）での空気過剰率λと、図３８のステップ５２０にて求めた指令燃料噴射量Qfinと、ステップ４１１０内に記載の式とに基づいて、指令燃料噴射量Qfinの燃料を含んだ定常火炎に係る混合気であって空気過剰率が値λとなる混合気の質量（定常混合気質量Mmixsteady）を求める。
次に、ＣＰＵ６１はステップ４１１５に進んで、図４０のステップ８５５にて更新されている現時点での混合気内のスート濃度[Soot]mixの値に上記定常混合気質量Mmixsteadyを乗じることで定常火炎が占める領域内でのスート発生量を求め、同求めたスート発生量に初期値Soot0を加えることでスート総発生量Sootを求める。
同様に、ＣＰＵ６１はステップ４１１５にて、図４０のステップ８８０にて更新されている現時点での混合気内のＮＯ濃度[NO]mixの値に上記定常混合気質量Mmixsteadyを乗じることで定常火炎が占める領域内でのＮＯ発生量を求め、同求めたＮＯ発生量に初期値ＮＯ0を加えることでＮＯ総発生量ＮＯを求める。
そして、ＣＰＵ６１はステップ９２５〜９６０を順に実行し、図４１のステップ４１１５にて求められているスート総発生量Sootの値、及びＮＯ総発生量ＮＯの値（従って、エミッション総発生量の推定結果）に基づいて噴射圧力を補正し、ステップ３８９５に進んで図３８〜図４１の一連の本ルーチンを一旦終了する。以降、ＣＰＵ６１は、次のＩＶＣが到来するまでの間、図３８のステップ５０５に進む毎に「Ｎｏ」と判定し続ける。
以上、説明したように、本発明によるエミッション発生量推定装置の第４実施形態は、混合気の着火後において燃焼室内で所謂定常火炎が発生している場合を想定している。そして、第４実施形態は、定常火炎のリッチ領域内（λ＜１）では、定常状態における酸素濃度（定常酸素濃度[O2]mixsteady）が「０」になることを利用して、上記(23)式においてdmso/dtの項を省略することで得られる式である「d[Soot]mix/dt＝dmsf/dt」に従ってスート発生速度d[Soot]mix/dtを求める。
同様に、第４実施形態は、定常火炎のリーン領域内（λ≧１）では、定常状態における燃料濃度（定常燃料濃度[Fuel]mixsteady）が「０」になることを利用して、上記(23)式においてdmsf/dtの項を省略することで得られる式である「d[Soot]mix/dt＝−dmso/dt」に従ってスート発生速度d[Soot]mix/dtを求める。これらにより、(24)式、(25)式の算出回数を減らすことができ、この結果、ＣＰＵ６１の計算負荷を減らすことができる。
本発明の第１実施形態に係る内燃機関の混合気状態推定装置を含んだエミッション発生量推定装置を４気筒内燃機関（ディーゼル機関）に適用したシステム全体の概略構成図である。 或る一つの気筒のシリンダ内(筒内)に吸気マニホールドからガスが吸入され、筒内に吸入されたガスが排気マニホールドへ排出される様子を模式的に示した図である。 （ａ）は、燃料が一時に噴射された時点での同燃料の様子を模式的に示した図であり、（ｂ）は、筒内ガスと混ざり合いながら混合気となって円錐状に拡散していく燃料の様子を模式的に示した図である。 噴射期間を前期・中期・後期の３つの期間に均等に区分し、各期間に対応して順次噴射されるそれぞれの燃料の部分に基づいて形成されていくそれぞれの混合気の様子を模式的に示した図である。 図１に示したＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第５番目の部分を示したフローチャートである。 図１に示したＣＰＵが図５に示したルーチンを実行する際に参照する指令燃料噴射量を決定するためのテーブルである。 図１に示したＣＰＵが図５に示したルーチンを実行する際に参照する燃料燃料噴射時期を決定するためのテーブルである。 図１に示したＣＰＵが図５に示したルーチンを実行する際に参照する基本燃料噴射圧力を決定するためのテーブルである。 図１に示したＣＰＵが実行する燃料噴射制御を行うためのルーチンを示したフローチャートである。 噴射期間を着火前・着火後の２つの期間に区分し、各期間に対応して順次噴射されるそれぞれの燃料の部分に基づいて形成されていくそれぞれの混合気の様子を模式的に示した図である。 第２実施形態のＣＰＵが図５に示したルーチンを実行する際における図５に示したルーチンの変更部分を示したフローチャートである。 第２実施形態のＣＰＵが図７に示したルーチンを実行する際における図７に示したルーチンの変更部分を示したフローチャートである。 第２実施形態のＣＰＵが図９に示したルーチンを実行する際における図９に示したルーチンの変更部分を示したフローチャートである。 噴射期間をn個の期間に区分し、n個の期間に対応して順次噴射されるn個の燃料の部分に基づいて形成されていくn個の混合気の時間経過に対する様子（噴射開始から噴射終了までの様子）を模式的に示した図である。 第３実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第５番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行するi(2≦i≦n)番目噴射に関するエミッション量の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行するi(2≦i≦n)番目噴射に関するエミッション量の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行するi(2≦i≦n)番目噴射に関するエミッション量の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。 第３実施形態のＣＰＵが実行するi(2≦i≦n)番目噴射に関するエミッション量の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。 燃焼室内にて定常火炎が発生している場合における、リッチ領域及びリーン領域と混合気到達距離との関係を示した図である。 定常火炎が占める領域内における、空気過剰率（混合気到達距離）と定常混合気温度上昇量との関係を示したグラフである。 定常火炎が占める領域内における、空気過剰率（混合気到達距離）と定常燃料濃度との関係を示したグラフである。 定常火炎が占める領域内における、空気過剰率（混合気到達距離）と定常酸素濃度との関係を示したグラフである。 定常火炎が占める領域内における、空気過剰率（混合気到達距離）と定常窒素濃度との関係を示したグラフである。 定常火炎が占める領域内における、噴射後経過時間と空気過剰率との関係を示したグラフである。 定常火炎が占める領域内における、噴射後経過時間と定常混合気温度上昇量との関係を示したグラフである。 定常火炎が占める領域内における、噴射後経過時間と定常燃料濃度との関係を示したグラフである。 定常火炎が占める領域内における、噴射後経過時間と定常酸素濃度との関係を示したグラフである。 定常火炎が占める領域内における、噴射後経過時間と定常窒素濃度との関係を示したグラフである。 第４実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。 第４実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。 第４実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。 第４実施形態のＣＰＵが実行する混合気温度等、及びエミッション量の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。
符号の説明
２１…燃料噴射弁、６０…電気制御装置、６１…ＣＰＵ、７２…吸気温センサ、７３…吸気管圧力センサ、７４…クランクポジションセンサ、７６…燃料温度センサ、７７…筒内圧力センサ、７８…吸気酸素濃度センサ

Claims (7)

  1. 内燃機関の燃焼室内にて所定の噴射開始時点から所定の噴射期間だけ連続して噴射される燃料を複数の部分に区分する噴射燃料区分手段と、
    前記区分された噴射燃料の各部分が前記所定の噴射開始時点からの時間経過に従って独立して順次噴射されていくとの仮定のもと、前記噴射燃料の各部分が前記燃焼室内に吸入されている筒内ガスと混ざり合って形成されていくそれぞれの混合気の状態を個別に推定する混合気状態推定手段と、
    前記混合気の着火時期を表す値を取得する着火時期取得手段と、
    を備え、
    前記噴射燃料区分手段は、
    前記取得された着火時期を表す値に基づいて前記噴射燃料を前記混合気の着火前に噴射される部分と同混合気の着火後に噴射される部分とに区分するように構成された内燃機関の混合気状態推定装置。
  2. 請求項１に記載の内燃機関の混合気状態推定装置において、
    前記混合気状態推定手段は、
    前記区分された噴射燃料の各部分に対して前記筒内ガスが混ざり合っていく程度を表す値を個別に取得し、対応する部分に対して前記取得された前記筒内ガスが混ざり合っていく程度を表す値を使用して各部分に基づいて形成されていくそれぞれの混合気の状態を個別に推定するように構成された内燃機関の混合気状態推定装置。
  3. 請求項１又は請求項２に記載の内燃機関の混合気状態推定装置において、
    前記混合気状態推定手段により推定される混合気の状態は、同混合気の温度、同混合気内の燃料濃度、同混合気内の酸素濃度、同混合気内の窒素濃度の少なくとも一つを含んでいる内燃機関の混合気状態推定装置。
  4. 請求項１乃至請求項３の何れか一項に記載の内燃機関の混合気状態推定装置により推定された対応する混合気の状態に基づいて、前記それぞれの混合気の燃焼に起因して発生する有害物質であるエミッションの発生量を個別に推定するエミッション発生量推定手段を備えた内燃機関のエミッション発生量推定装置。
  5. 請求項４に記載の内燃機関のエミッション発生量推定装置において、
    前記エミッション発生量推定手段は、
    前記エミッションの発生量として、少なくともスート、及びＮＯｘの発生量を推定するように構成された内燃機関のエミッション発生量推定装置。
  6. 請求項４又は請求項５に記載の内燃機関のエミッション発生量推定装置において、
    前記エミッション発生量推定手段は、
    前記それぞれの混合気に対応する前記エミッションの発生速度を個別に算出し、同算出されたエミッションの発生速度を時間積分することで前記エミッションの発生量を個別に推定するように構成された内燃機関のエミッション発生量推定装置。
  7. 請求項６に記載の内燃機関のエミッション発生量推定装置において、
    前記エミッション発生量推定手段は、
    前記混合気の温度が所定温度未満となったとき、対応する前記エミッションの発生速度の時間積分処理を終了する内燃機関のエミッション発生量推定装置。
JP2005098148A 2005-03-30 2005-03-30 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置 Expired - Fee Related JP4424242B2 (ja)

Priority Applications (5)

Application Number Priority Date Filing Date Title
JP2005098148A JP4424242B2 (ja) 2005-03-30 2005-03-30 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置
PCT/JP2006/306769 WO2006106893A1 (ja) 2005-03-30 2006-03-24 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置
CN2006800106974A CN101151453B (zh) 2005-03-30 2006-03-24 用于内燃机的气体混合物状态估计设备和排放物生成量估计设备
US11/886,595 US7848871B2 (en) 2005-03-30 2006-03-24 Gas-mixture-state estimation apparatus for internal combustion engine and emission-generation-amount estimation apparatus
EP06730717.3A EP1865178B1 (en) 2005-03-30 2006-03-24 Mixture gas state estimating device of internal combustion engine, and emission amount estimating device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2005098148A JP4424242B2 (ja) 2005-03-30 2005-03-30 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置

Publications (2)

Publication Number Publication Date
JP2006274991A JP2006274991A (ja) 2006-10-12
JP4424242B2 true JP4424242B2 (ja) 2010-03-03

Family

ID=37073446

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2005098148A Expired - Fee Related JP4424242B2 (ja) 2005-03-30 2005-03-30 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置

Country Status (5)

Country Link
US (1) US7848871B2 (ja)
EP (1) EP1865178B1 (ja)
JP (1) JP4424242B2 (ja)
CN (1) CN101151453B (ja)
WO (1) WO2006106893A1 (ja)

Families Citing this family (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8899209B2 (en) 2010-10-08 2014-12-02 Ford Global Technologies, Llc System and method for compensating cetane
US8051829B2 (en) * 2010-10-08 2011-11-08 Ford Global Technologies, Llc Method for controlling low temperature combustion
US9172217B2 (en) 2010-11-23 2015-10-27 Woodward, Inc. Pre-chamber spark plug with tubular electrode and method of manufacturing same
US8584648B2 (en) 2010-11-23 2013-11-19 Woodward, Inc. Controlled spark ignited flame kernel flow
US9476347B2 (en) 2010-11-23 2016-10-25 Woodward, Inc. Controlled spark ignited flame kernel flow in fuel-fed prechambers
JP5936369B2 (ja) * 2012-01-24 2016-06-22 三菱重工業株式会社 内燃機関の制御装置
US8949002B2 (en) 2012-02-21 2015-02-03 Ford Global Technologies, Llc System and method for injecting fuel
JP5907123B2 (ja) * 2012-07-13 2016-04-20 井関農機株式会社 スート堆積演算表示装置
US9551287B2 (en) 2012-10-30 2017-01-24 Mitsubishi Heavy Industries Ltd. Internal combustion engine control apparatus and control method thereof
US9856848B2 (en) 2013-01-08 2018-01-02 Woodward, Inc. Quiescent chamber hot gas igniter
JP5979012B2 (ja) * 2013-01-11 2016-08-24 トヨタ自動車株式会社 Ｐｍ排出量推定装置
US9765682B2 (en) 2013-06-10 2017-09-19 Woodward, Inc. Multi-chamber igniter
US9683497B2 (en) * 2013-10-25 2017-06-20 Ford Global Technologies, Llc Methods and systems for adjusting engine airflow based on output from an oxygen sensor
JP2015105572A (ja) * 2013-11-28 2015-06-08 株式会社デンソー 混合気の混合比推定装置
JP6057022B2 (ja) * 2014-05-27 2017-01-11 日産自動車株式会社 ディーゼルエンジンの制御装置および制御方法
AT516532B1 (de) * 2014-11-24 2019-10-15 Innio Jenbacher Gmbh & Co Og Verfahren zum Starten einer mit einem Brennstoff-Luft-Gemisch betriebenen Brennkraftmaschine
JP6222138B2 (ja) 2015-03-03 2017-11-01 トヨタ自動車株式会社 内燃機関のエミッション推定装置
US9653886B2 (en) 2015-03-20 2017-05-16 Woodward, Inc. Cap shielded ignition system
US9840963B2 (en) 2015-03-20 2017-12-12 Woodward, Inc. Parallel prechamber ignition system
FR3039214B1 (fr) * 2015-07-22 2019-04-05 Continental Automotive France Procede de gestion de l'injection dans un moteur a combustion interne
US9890689B2 (en) * 2015-10-29 2018-02-13 Woodward, Inc. Gaseous fuel combustion
JP6414156B2 (ja) * 2016-07-22 2018-10-31 トヨタ自動車株式会社 内燃機関の制御装置
US11118518B2 (en) * 2019-07-23 2021-09-14 Caterpillar Inc. Method and system for aftertreatment control
US11415041B2 (en) 2019-09-16 2022-08-16 Woodward, Inc. Flame triggered and controlled volumetric ignition
CN113867131A (zh) * 2021-09-30 2021-12-31 潍柴动力股份有限公司 一种egr控制的方法、装置及电子设备
CH719186A2 (de) * 2021-12-01 2023-06-15 Liebherr Machines Bulle Sa Verfahren zum Betrieb einer Verbrennungskraftmaschine mit einem gasförmigen Kraftstoff sowie Verbrennungskraftmaschine.

Family Cites Families (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4708117A (en) * 1986-04-14 1987-11-24 Colt Industries Inc. Multi-point fuel injection apparatus
JPH06108903A (ja) 1992-09-28 1994-04-19 Unisia Jecs Corp 内燃機関の燃焼制御装置
JP3383874B2 (ja) 1993-07-19 2003-03-10 克享 小西 ディーゼル機関の燃焼シミュレーション方法
JPH0783097A (ja) * 1993-09-13 1995-03-28 Honda Motor Co Ltd 内燃機関の空燃比検出方法
JPH07166922A (ja) 1993-12-13 1995-06-27 Nippon Soken Inc 内燃機関の燃料噴射制御装置
JP3562016B2 (ja) * 1994-09-06 2004-09-08 マツダ株式会社 自動車用リーンバーンエンジン
US5713328A (en) * 1997-03-31 1998-02-03 Ford Global Technologies, Inc. Spark ignited internal combustion engine with multiple event fuel injection
JP3591283B2 (ja) * 1998-01-29 2004-11-17 日産自動車株式会社 エンジンの排気浄化装置
DE19823021B4 (de) * 1998-05-22 2004-08-12 Fev Motorentechnik Gmbh Verfahren zum Betrieb einer fremdgezündeten Kolbenbrennkraftmaschine mit geregeltem Abgaskatalysator und elektromagnetisch betätigten Gaswechselventilen
US6390059B1 (en) * 1998-06-22 2002-05-21 Hitachi, Ltd. Cylinder-injection type internal combustion engine, method of controlling the engine, and fuel injection nozzle
JP3817977B2 (ja) * 1999-07-06 2006-09-06 株式会社日立製作所 圧縮着火式エンジンの制御方法
JP3684934B2 (ja) * 1999-08-30 2005-08-17 三菱自動車工業株式会社 内燃機関の排気浄化装置
JP3680259B2 (ja) 2000-03-08 2005-08-10 トヨタ自動車株式会社 ディーゼル機関の燃料噴射装置
JP2002054491A (ja) 2000-08-11 2002-02-20 Nissan Motor Co Ltd エンジンの燃焼室内温度状態推定装置
US6360531B1 (en) * 2000-08-29 2002-03-26 Ford Global Technologies, Inc. System and method for reducing vehicle emissions
JP4706134B2 (ja) 2001-06-15 2011-06-22 トヨタ自動車株式会社 内燃機関の制御装置
DE50110758D1 (de) * 2001-09-25 2006-09-28 Ford Global Tech Llc Vorrichtung und Verfahren zur Regeneration einer Abgasbehandlungseinrichtung
JP3893967B2 (ja) * 2001-12-18 2007-03-14 日産自動車株式会社 ディーゼルエンジンの制御装置
JP4089244B2 (ja) * 2002-03-01 2008-05-28 株式会社デンソー 内燃機関用噴射量制御装置
JP4161690B2 (ja) * 2002-11-20 2008-10-08 株式会社デンソー 蓄圧式燃料噴射装置
JP3956137B2 (ja) * 2003-09-18 2007-08-08 トヨタ自動車株式会社 内燃機関の混合気温度推定方法
JP2005127169A (ja) * 2003-10-22 2005-05-19 Hitachi Ltd 内燃機関の制御方法
JP4251123B2 (ja) 2003-11-04 2009-04-08 株式会社デンソー 内燃機関
JP3982508B2 (ja) * 2004-02-25 2007-09-26 トヨタ自動車株式会社 内燃機関の混合気温度推定装置
JP4027902B2 (ja) * 2004-03-24 2007-12-26 株式会社豊田中央研究所 内燃機関の混合気着火時期推定装置、及び内燃機関の制御装置
JP4692179B2 (ja) 2005-09-21 2011-06-01 トヨタ自動車株式会社 内燃機関の混合気濃度推定装置

Also Published As

Publication number Publication date
EP1865178A1 (en) 2007-12-12
EP1865178B1 (en) 2020-08-19
JP2006274991A (ja) 2006-10-12
US7848871B2 (en) 2010-12-07
CN101151453B (zh) 2010-08-18
EP1865178A4 (en) 2015-03-04
US20090055076A1 (en) 2009-02-26
WO2006106893A1 (ja) 2006-10-12
CN101151453A (zh) 2008-03-26

Similar Documents

Publication Publication Date Title
JP4424242B2 (ja) 内燃機関の混合気状態推定装置、及びエミッション発生量推定装置
JP4603951B2 (ja) 内燃機関のすす発生量推定装置
US7281368B2 (en) Nox discharge quantity estimation method for internal combustion engine
US7337766B2 (en) Gas-mixture-ignition-time estimation apparatus for internal combustion engine, and control apparatus for internal combustion engine
US20070255484A1 (en) Control Apparatus for Internal Combustion Engine
EP1408331A2 (en) Real-time nitrogen oxides (NOx) estimation process
EP1529941B1 (en) NOx generation quantity estimation method for internal combustion engine
CN102159804B (zh) 内燃机的炭烟排放量估计装置
US9926873B2 (en) Internal combustion engine control apparatus
US10221793B2 (en) Controller for internal combustion engine
US20160123264A1 (en) Control device for internal combustion engine
WO2010140263A1 (ja) 内燃機関のすす排出量推定装置
EP1544443B1 (en) Method of estimating the in cylinder temperature after combustion
EP2757238A1 (en) Control device for internal combustion engine
JP4692179B2 (ja) 内燃機関の混合気濃度推定装置
JP5077491B2 (ja) 内燃機関の制御装置
EP2610469A1 (en) Device for estimating diffuse combustion start time and device for controlling diffuse combustion start time for internal combustion engine
JP2007092583A (ja) 内燃機関の燃料噴射時期決定装置
JP2005180219A (ja) 内燃機関のＮＯｘ発生量推定モデルの学習方法
JP4775303B2 (ja) 内燃機関の燃料噴射時期決定装置

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20070907

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20090909

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20091026

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20091117

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20091130

R151 Written notification of patent or utility model registration

Ref document number: 4424242

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20121218

Year of fee payment: 3

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20131218

Year of fee payment: 4

LAPS Cancellation because of no payment of annual fees