(57)【特許請求の範囲】
【請求項１】
内燃機関の燃焼室内において燃料の反応に起因して発生するすすの発生量を推定するすす発生量推定手段を備えた内燃機関のすす発生量推定装置において、
前記すす発生量推定手段は、
前記燃料から前記すすが発生する反応過程についての反応モデルを用いて前記すすの発生速度を算出するすす発生速度算出手段を備え、
前記算出されたすすの発生速度に基づいて前記すすの発生量を推定するように構成され、
前記反応モデルとして、前記燃料から前記すすの前駆物質を経由して前記すすが発生する反応過程についてのモデルが用いられ、
前記すす発生量推定手段は、
前記燃料から前記前駆物質が発生する速度を算出する前駆物質発生速度算出手段を更に備え、
前記すす発生速度算出手段は、
前記算出された前駆物質の発生速度に基づいて前記すすの発生速度を算出するように構成され、
前記前駆物質発生速度算出手段は、
前記燃料から前記前駆物質が生成される速度を算出する前駆物質生成速度算出手段と、
前記生成された前駆物質が熱分解により分解される速度を算出する前駆物質分解速度算出手段と、
前記前駆物質から前記すすが生成される速度を算出するすす生成速度算出手段と、
を備え、
前記算出された前駆物質の生成速度、前駆物質の分解速度、及びすすの生成速度に基づいて前記前駆物質の発生速度を算出するように構成され、
前記燃料と前記燃焼室内に吸入されている筒内ガスとが混ざり合って形成されている混合気の温度が所定温度以上のとき、前記反応モデルとして、前記燃料から前記すすの前駆物質を経由して前記すすが発生する反応過程についてのモデルに代えて、前記燃料から前記すすが直接発生する反応過程についてのモデルが用いられる、内燃機関のすす発生量推定装置。
【請求項２】
請求項１に記載の内燃機関のすす発生量推定装置において、
前記すす発生速度算出手段は、
前記生成されたすすが酸化される速度を算出するすす酸化速度算出手段を備え、
前記混合気の温度が前記所定温度未満のとき、前記すす生成速度算出手段により算出されたすすの生成速度、及び前記すす酸化速度算出手段により算出されたすすの酸化速度に基づいて前記すすの発生速度を算出するように構成された内燃機関のすす発生量推定装置。
【請求項３】
請求項１又は請求項２に記載の内燃機関のすす発生量推定装置において、
前記すす生成速度算出手段は、
前記前駆物質同士の衝突に起因して前記すすが生成される速度である第１すす生成速度を算出する第１すす生成速度算出手段と、
前記前駆物質と前記すすの衝突に起因して新たな前記すすが生成される速度である第２すす生成速度を算出する第２すす生成速度算出手段と、
を備え、
前記混合気の温度が前記所定温度未満のとき、前記算出された第１すす生成速度、及び第２すす生成速度に基づいて前記すす生成速度を算出するように構成された内燃機関のすす発生量推定装置。
【請求項４】
請求項３に記載の内燃機関のすす発生量推定装置において、
前記第１すす生成速度算出手段は、
前記前駆物質同士の衝突に起因して前記すすが生成される確率を考慮して前記第１すす生成速度を算出するように構成され、
前記第２すす生成速度算出手段は、
前記前駆物質と前記すすの衝突に起因して新たな前記すすが生成される確率を考慮して前記第２すす生成速度を算出するように構成された内燃機関のすす発生量推定装置。
【請求項５】
請求項４に記載の内燃機関のすす発生量推定装置において、
前記第１すす生成速度算出手段は、
前記前駆物質同士の衝突数、及び／又は、衝突した前記前駆物質同士が凝集する確率を考慮して前記第１すす生成速度を算出するように構成され、
前記第２すす生成速度算出手段は、
前記前駆物質と前記すすの衝突数、及び／又は、衝突した前記前駆物質と前記すすとが凝集する確率を考慮して前記第２すす生成速度を算出するように構成された内燃機関のすす発生量推定装置。
【請求項６】
請求項５に記載の内燃機関のすす発生量推定装置において、
前記第１すす生成速度算出手段は、
前記燃料と前記燃焼室内に吸入されている筒内ガスとが混ざり合って形成されている混合気中の前記前駆物質の濃度と、前記混合気の温度とに基づいて前記前駆物質同士の衝突数を算出するように構成され、
前記第２すす生成速度算出手段は、
前記混合気中の前記前駆物質の濃度と、前記混合気中のすすの濃度と、前記混合気の温度とに基づいて前記前駆物質と前記すすの衝突数を算出するように構成された内燃機関のすす発生量推定装置。
【請求項７】
請求項５又は請求項６に記載の内燃機関のすす発生量推定装置において、
前記衝突した前記前駆物質と前記すすとが凝集する確率は、前記衝突した前記前駆物質同士が凝集する確率よりも大きい内燃機関のすす発生量推定装置。
【請求項８】
請求項１乃至請求項７の何れか一項に記載の内燃機関のすす発生量推定装置において、
前記前駆物質は、ＰＡＨ、芳香族中間生成物、飽和重合体、及びＣ２Ｈ２の少なくとも一つを含む単量体であり、前記すすは、前記単量体の２分子以上が凝集してなる重合体である内燃機関のすす発生量推定装置。
【請求項９】
請求項１に記載の内燃機関のすす発生量推定装置において、
前記すす発生量推定手段は、
前記混合気の温度が前記所定温度以上のとき、前記混合気の温度に基づいて前記燃料から前記すすが生成される速度の値を求める関数であって前記混合気の温度が高くなるほど減少する値を採る関数を用いて前記すすの生成速度を算出する高温時すす生成速度算出手段を更に備え、
前記すす発生速度算出手段は、前記混合気の温度が前記所定温度以上のとき、前記高温時すす生成速度算出手段により算出された前記すすの生成速度に基づいて前記すすの発生速度を算出するように構成された内燃機関のすす発生量推定装置。
【請求項１０】
請求項１乃至請求項９の何れか一項に記載の内燃機関のすす発生量推定装置において、
前記すす発生量推定手段は、
前記燃料と前記燃焼室内に吸入されている筒内ガスとが混ざり合って形成されている混合気の当量比が所定値以上の場合にのみ、前記すすの発生量の推定に係る計算を行うように構成された内燃機関のすす発生量推定装置。
【請求項１１】
請求項１０に記載の内燃機関のすす発生量推定装置において、
前記すす発生量推定手段は、
前記混合気の当量比が前記所定値以上であって、且つ、前記混合気の当量比と同混合気の温度との関係が所定の範囲内にある場合にのみ、前記すすの発生量の推定に係る計算を行うように構成された内燃機関のすす発生量推定装置。
【請求項１２】
請求項１乃至請求項１１の何れか一項に記載の内燃機関のすす発生量推定装置において、
前記すす発生量推定手段は、
前記燃料の性状を考慮して前記すすの発生量を推定するように構成された内燃機関のすす発生量推定装置。
【請求項１３】
請求項１２に記載の内燃機関のすす発生量推定装置において、
前記すす発生量推定手段は、
前記燃料の性状として、前記燃料に含まれる芳香族炭化水素の成分の含有率を考慮して前記すすの発生量を推定するように構成された内燃機関のすす発生量推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、内燃機関の燃焼室内において燃料の反応に起因して発生するすす（カーボン微粒子。以下、「Ｓｏｏｔ」とも称呼する。）の発生量を推定するすす発生量推定装置に関する。
【背景技術】
【０００２】
内燃機関（特に、ディーゼル機関）の燃焼室内で発生する粒子状物質（パティキュレート・マター（ＰＭ））を構成する主たる成分の一つはＳｏｏｔである。このＳｏｏｔの発生量を精度良く制御してＳｏｏｔの発生量を少なくするためには、Ｓｏｏｔの発生量を精度良く推定する必要がある。
【０００３】
例えば、下記特許文献１に記載の内燃機関の排気浄化装置は、エンジン回転速度、アクセル開度、冷却水の温度等を引数とするマップを利用して機関の運転状態に応じて変化するＰＭ発生量を推定するようになっている。このように、ＰＭ発生量が推定できれば、ＰＭ内のＳｏｏｔの発生量も推定することができる。
【特許文献１】特開２００４−４４４５７号公報
【０００４】
ところで、近年、燃料からＳｏｏｔが発生する反応過程におけるＳｏｏｔの生成メカニズムについて多くの研究がなされている。係るＳｏｏｔの生成メカニズムに基づく反応モデルを構築し、同反応モデルを用いてＳｏｏｔの発生量を推定すれば、同Ｓｏｏｔの発生量がより一層精度良く推定され得ることが大いに期待される。
【発明の開示】
【０００５】
本発明の目的は、Ｓｏｏｔの生成メカニズムに基づく反応モデルを用いて燃焼室内で発生するＳｏｏｔの発生量を精度良く推定できるすす発生量推定装置を提供することにある。
【０００６】
本発明に係るすす発生量推定装置の特徴は、内燃機関の燃焼室内において燃料の反応に起因して発生するすすの発生量を推定するすす発生量推定手段が、前記燃料から前記すすが発生する反応過程についての反応モデル（を表す式）を用いて前記すすの発生速度を算出するすす発生速度算出手段を備え、前記算出されたすすの発生速度に基づいて前記すすの発生量を推定するように構成されたことにある。
【０００７】
これによれば、燃料からＳｏｏｔが発生する反応過程についてのＳｏｏｔの生成メカニズムに基づく反応モデルを表す式を用いてＳｏｏｔの発生速度が算出される。Ｓｏｏｔの発生速度とは、例えば、混合気内のＳｏｏｔ量の変化速度、混合気内のＳｏｏｔ濃度の変化速度等である。
【０００８】
そして、係るＳｏｏｔの発生速度に基づいて（具体的には、Ｓｏｏｔの発生速度を時間積分（積算）することで）Ｓｏｏｔの発生量が推定される。Ｓｏｏｔの発生量とは、例えば、混合気内のＳｏｏｔ量、混合気内のＳｏｏｔ濃度等である。このように、Ｓｏｏｔの発生量がＳｏｏｔの生成メカニズムに基づく反応モデルを表す式を用いて推定されるから、Ｓｏｏｔの発生量が精度良く推定され得る。
【０００９】
ところで、近年におけるＳｏｏｔの生成メカニズムについての研究によれば、燃料からＳｏｏｔの前駆物質（precursor）を経由してＳｏｏｔが発生する反応過程が明らかになりつつある。従って、係る前駆物質が考慮された反応過程についての反応モデルを構築し、同反応モデルを表す式を用いてＳｏｏｔの発生量を推定すれば、同Ｓｏｏｔの発生量がより一層精度良く推定され得ることが期待される。
【００１０】
即ち、上記本発明に係るすす発生量推定装置においては、前記反応モデルとして、前記燃料から前記すすの前駆物質を経由して前記すすが発生する反応過程についてのモデルが用いられることが好適である。ここにおいて、前記前駆物質は、ＰＡＨ、芳香族中間生成物、飽和重合体、及びＣ２Ｈ２の少なくとも一つを含む単量体と定義され、Ｓｏｏｔは、単量体である前記前駆物質の２分子以上が凝集してなる重合体と定義される。前記前駆物質は、凝集により後にＳｏｏｔに成長し得る単量体の全てを含んでいる。
【００１１】
この場合、前記すす発生量推定手段は、前記燃料から前記前駆物質が発生する速度を算出する前駆物質発生速度算出手段を更に備え、前記すす発生速度算出手段は、前記算出された前駆物質の発生速度に基づいて前記すすの発生速度を算出するように構成されることが好適である。ここにおいて、前駆物質の発生速度とは、例えば、混合気内の前駆物質量の変化速度、混合気内の前駆物質濃度の変化速度等である。
【００１２】
これによれば、例えば、前駆物質発生速度算出手段により算出された前駆物質の発生速度に基づいて（具体的には、前記前駆物質の発生速度を時間積分（積算）することで）前駆物質の発生量が算出され得る。ここで、前駆物質の発生量とは、例えば、混合気内の前駆物質量、混合気内の前駆物質濃度等である。そして、係る前駆物質の発生量に基づいて上記Ｓｏｏｔの発生速度が算出され、これにより、上記Ｓｏｏｔの発生量が推定される。
【００１３】
このように、燃料からＳｏｏｔが発生する反応過程においてＳｏｏｔの前駆物質を経由させた反応モデルを用いることで、燃料からＳｏｏｔが発生する反応過程において実験的・経験的に知られている「Ｓｏｏｔ発生（反応）の遅れ」を精度良く模擬することができる。加えて、後述するように「前駆物質の熱分解による同前駆物質の分解速度」を考慮して上記前駆物質の発生速度を算出することで、実験的・経験的に知られている「混合気温度が或る温度を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」をも精度良く模擬することができる。
【００１４】
上記本発明に係るすす発生量推定装置において、前記前駆物質発生速度算出手段が備えられている場合、前記前駆物質発生速度算出手段は、前記燃料から前記前駆物質が生成される速度を算出する前駆物質生成速度算出手段と、前記生成された前駆物質が熱分解により分解される速度を算出する前駆物質分解速度算出手段と、前記前駆物質から前記すすが生成される速度を算出するすす生成速度算出手段とを備え、前記算出された前駆物質の生成速度、前駆物質の分解速度、及びすすの生成速度に基づいて前記前駆物質の発生速度を算出するように構成されることが好ましい。
【００１５】
燃料から前駆物質が生成される反応は、前駆物質発生速度を増大させる方向に働く。一方、生成された前駆物質は熱分解により分解され得る。係る前駆物質の熱分解反応、並びに生成された前駆物質からＳｏｏｔが生成される反応は、前駆物質発生速度を減少させる方向に働く。即ち、上記算出された前駆物質の生成速度、前駆物質の分解速度、及びＳｏｏｔの生成速度に基づいて前駆物質発生速度を算出するように構成すれば、これら３つの反応が全て考慮されて前駆物質発生速度が算出されるから、同前駆物質発生速度が精度良く算出され得る。
【００１６】
或いは、前記前駆物質発生速度算出手段は、前記燃料から前記前駆物質が生成される速度から前記生成された前駆物質が熱分解により分解される速度を減じた前駆物質の正味の生成速度を算出する前駆物質正味生成速度算出手段と、前記前駆物質から前記すすが生成される速度を算出するすす生成速度算出手段とを備え、前記算出された前駆物質の正味生成速度、及びすすの生成速度に基づいて前記前駆物質の発生速度を算出するように構成されてもよい。
【００１７】
これによれば、上記前駆物質の生成速度、及び上記前駆物質の分解速度を別個独立に算出することなく、上記前駆物質の生成速度から上記前駆物質の分解速度を減じた値である前駆物質の正味の生成速度が算出され、この前駆物質の正味の生成速度と、Ｓｏｏｔの生成速度とに基づいて前駆物質発生速度が算出される。従って、前駆物質発生速度の算出に際し、上記３つの反応を全て考慮しつつ、計算負荷を少なくすることができる。
【００１８】
上記本発明に係るすす発生量推定装置において、前記すす発生速度算出手段が前記前駆物質の発生速度に基づいて前記すすの発生速度を算出するように構成されている場合、同すす発生速度算出手段は、前記前駆物質から前記すすが生成される速度を算出するすす生成速度算出手段と、前記生成されたすすが酸化される速度を算出するすす酸化速度算出手段とを備え、前記算出されたすすの生成速度、及びすすの酸化速度に基づいて前記すすの発生速度を算出するように構成されることが好適である。
【００１９】
前駆物質からＳｏｏｔが生成される反応は、Ｓｏｏｔ発生速度を増大させる方向に働く。一方、生成されたＳｏｏｔは酸化され得る。係るＳｏｏｔの酸化反応は、Ｓｏｏｔ発生速度を減少させる方向に働く。即ち、上記算出されたＳｏｏｔの生成速度、及びＳｏｏｔの酸化速度に基づいてＳｏｏｔ発生速度を算出するように構成すれば、これら２つの反応が考慮されてＳｏｏｔ発生速度が算出されるから、同Ｓｏｏｔ発生速度が精度良く算出され得る。
【００２０】
上記本発明に係るすす発生量推定装置において、前記すす生成速度算出手段が備えられている場合、前記すす生成速度算出手段は、前記前駆物質同士の衝突に起因して前記すすが生成される速度である第１すす生成速度を算出する第１すす生成速度算出手段と、前記前駆物質と前記すすの衝突に起因して新たな前記すすが生成される速度である第２すす生成速度を算出する第２すす生成速度算出手段とを備え、前記算出された第１すす生成速度、及び第２すす生成速度に基づいて前記すす生成速度を算出するように構成されることが好適である。
【００２１】
前駆物質からＳｏｏｔが生成される反応過程においては、前駆物質同士が衝突することで凝集によりＳｏｏｔが生成される第１反応と、Ｓｏｏｔが生成された後において前駆物質とＳｏｏｔとが衝突することで凝集により新たなＳｏｏｔが生成される第２反応とが考えられる。上記構成によれば、これら２つの反応にそれぞれ対応する第１Ｓｏｏｔ生成速度、及び第２Ｓｏｏｔ生成速度が別個独立に考慮されてＳｏｏｔ生成速度が算出されるから、同Ｓｏｏｔ生成速度が精度良く算出され得る。
【００２２】
この場合、前記第１すす生成速度算出手段は、前記前駆物質同士の衝突に起因して前記すすが生成される確率を考慮して前記第１すす生成速度を算出するように構成され、前記第２すす生成速度算出手段は、前記前駆物質と前記すすの衝突に起因して新たな前記すすが生成される確率を考慮して前記第２すす生成速度を算出するように構成されることが好適である。
【００２３】
より具体的には、この場合、前記第１すす生成速度算出手段は、前記前駆物質同士の衝突数（混合気単位質量あたりの衝突数）、及び／又は、衝突した前記前駆物質同士が凝集する確率を考慮して前記第１すす生成速度を算出するように構成され、前記第２すす生成速度算出手段は、前記前駆物質と前記すすの衝突数（混合気単位質量あたりの衝突数）、及び／又は、衝突した前記前駆物質と前記すすとが凝集する確率を考慮して前記第２すす生成速度を算出するように構成される。
【００２４】
上記第１Ｓｏｏｔ生成速度は、前駆物質同士の衝突に起因してＳｏｏｔが生成される確率に依存する。この確率は、前駆物質同士の衝突数（混合気単位質量あたりの衝突数）と、衝突した前駆物質同士が凝集する確率とにより表すことができる。同様に、上記第２Ｓｏｏｔ生成速度は、前駆物質とＳｏｏｔの衝突に起因して新たなＳｏｏｔが生成される確率に依存する。この確率は、前駆物質とＳｏｏｔの衝突数（混合気単位質量あたりの衝突数）と、衝突した前駆物質とＳｏｏｔとが凝集する確率とにより表すことができる。
【００２５】
従って、上記構成によれば、第１Ｓｏｏｔ生成速度、及び第２Ｓｏｏｔ生成速度がより一層精度良く算出され得るから、Ｓｏｏｔ発生速度がより一層精度良く算出され得、この結果、Ｓｏｏｔ発生量がより一層精度良く推定され得る。
【００２６】
ここで、一般に、分子間同士の凝集は、分子間力が大きいほど、従って、分子の質量が大きいほど発生し易くなる。重合体であるＳｏｏｔの分子の質量は、単量体である前駆物質の分子の質量よりも大きい。従って、衝突した前駆物質とＳｏｏｔとが凝集する確率は、衝突した前駆物質同士が凝集する確率よりも大きい値に設定されることが好ましい。
【００２７】
更には、前記第１すす生成速度算出手段は、前記燃料と前記燃焼室内に吸入されている筒内ガスとが混ざり合って形成されている混合気中の前記前駆物質の濃度と、前記混合気の温度とに基づいて前記前駆物質同士の衝突数を算出するように構成され、前記第２すす生成速度算出手段は、前記混合気中の前記前駆物質の濃度と、前記混合気中のすすの濃度と、前記混合気の温度とに基づいて前記前駆物質と前記すすの衝突数を算出するように構成されると好ましい。
【００２８】
一般に、ガス中における物質Ａと物質Ｂのガス単位質量あたりの衝突数は、ガス中の物質Ａの濃度と、ガス中の物質Ｂの濃度と、ガスの温度とに依存する。従って、上述した前駆物質同士の混合気単位質量あたりの衝突数は、混合気中の前駆物質の濃度と、混合気の温度とに依存する。同様に、上述した前駆物質とＳｏｏｔの混合気単位質量あたりの衝突数は、混合気中の前駆物質の濃度と、混合気中のＳｏｏｔの濃度と、混合気の温度とに依存する。即ち、上記構成によれば、前駆物質同士の（混合気単位質量あたりの）衝突数、及び前駆物質とＳｏｏｔの（混合気単位質量あたりの）衝突数を精度良く算出することができる。
【００２９】
ところで、上記本発明に係るすす発生量推定装置において前記燃料から前記すすの前駆物質を経由して前記すすが発生する反応過程についての反応モデルが用いられる場合、混合気の温度が所定温度以上のとき、前記反応モデルとして、前記燃料から前記すすの前駆物質を経由して前記すすが発生する反応過程についてのモデルに代えて、前記燃料から前記すすが直接発生する反応過程についてのモデルが用いられてもよい。
【００３０】
混合気の温度が高温になると、燃料からＳｏｏｔが発生する反応過程における上述した「Ｓｏｏｔ発生（反応）の遅れ」の程度が小さくなる傾向がある。換言すれば、混合気の温度が高温になると、燃料からＳｏｏｔが発生する反応過程においてＳｏｏｔの前駆物質を経由させた反応モデルに代えて燃料からＳｏｏｔが直接発生する反応過程についてのモデルを用いても、Ｓｏｏｔの発生量を十分に精度良く推定することができる。
【００３１】
上記構成は、係る知見に基づくものである。これによれば、混合気の温度が高温（所定温度以上）になると、Ｓｏｏｔ発生量の推定精度の低下を抑制しつつＳｏｏｔ発生量の推定に係る計算負荷を少なくすることができる。
【００３２】
この場合、前記すす発生量推定手段は、前記混合気の温度が前記所定温度以上のとき、前記混合気の温度に基づいて前記燃料から前記すすが生成される速度の値を求める関数であって前記混合気の温度が高くなるほど減少する値を採る関数を用いて前記すすの生成速度を算出する高温時すす生成速度算出手段を更に備え、前記すす発生速度算出手段は、前記混合気の温度が前記所定温度以上のとき、前記高温時すす生成速度算出手段により算出された前記すすの生成速度に基づいて前記すすの発生速度を算出するように構成されることが好適である。
【００３３】
上述したように、「混合気温度が或る温度を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」が実験的・経験的に知られている。従って、混合気温度が所定温度以上のとき、上記構成のように、混合気温度に基づいて燃料からＳｏｏｔが生成される速度の値を求める関数であって混合気温度が高くなるほど減少する値を採る関数を用いてＳｏｏｔ生成速度を算出するように構成すれば、混合気温度が所定温度以上のとき、上記現象を精度良く模擬することができる。換言すれば、混合気温度が所定温度以上のとき、燃料からＳｏｏｔが直接発生する反応過程についてのモデルを用いて計算負荷を少なくしつつ、Ｓｏｏｔ発生速度（従って、Ｓｏｏｔ発生量）を十分に精度良く推定することができる。
【００３４】
ところで、本発明に係るすす発生量推定装置においては、前記反応モデルとして、前記燃料から前記すすが直接発生する反応過程についてのモデルが常に用いられてもよい。この場合、前記すす発生量推定手段は、前記燃料と前記燃焼室内に吸入されている筒内ガスとが混ざり合って形成されている混合気の温度に基づいて前記燃料から前記すすが生成される速度の値を求める関数であって、前記混合気の温度が特定温度となるときに最大値を採り、同混合気の温度が同特定温度から離れるに従って同最大値から減少していく値を採る関数を用いて、前記すすの生成速度を算出するすす生成速度算出手段を更に備え、前記すす発生速度算出手段は、前記算出されたすすの生成速度に基づいて前記すすの発生速度を算出するように構成されることが好適である。
【００３５】
一般に、Ｓｏｏｔ発生速度は、混合気温度が或る温度に達するまでは同混合気温度の上昇に応じて増大していき、混合気温度が同或る温度を超えると、同混合気温度の上昇に応じて減少していく。従って、上記構成のように、混合気温度に基づいて燃料からＳｏｏｔが生成される速度の値を求める関数であって混合気温度が特定温度となるときに最大値を採り同混合気温度が同特定温度から離れるに従って同最大値から減少していく値を採る関数を用いてＳｏｏｔ生成速度を算出すると、Ｓｏｏｔ発生速度（従って、Ｓｏｏｔ発生量）を或る程度精度良く推定することができる。
【００３６】
この場合、前記すす生成速度算出手段は、前記燃料から前記すすが発生する反応過程における時間遅れを考慮して前記すすの生成速度を算出するように構成されることが好適である。これによれば、上述した実験的・経験的に知られている「Ｓｏｏｔ発生（反応）の遅れ」を模擬することができる。この結果、混合気温度にかかわらず、燃料からＳｏｏｔが直接発生する反応過程についてのモデルを用いて計算負荷を少なくしつつ、Ｓｏｏｔ発生速度（従って、Ｓｏｏｔ発生量）を十分に精度良く推定することができる。
【００３７】
上記本発明に係るすす発生量推定装置においては、前記すす発生量推定手段は、混合気の当量比が所定値以上の場合にのみ、前記すすの発生量の推定に係る計算を行うように構成されることが好適である。Ｓｏｏｔは、混合気の当量比が小さいと殆ど生成されないことが実験的・経験的に知られている。
【００３８】
従って、上記構成によれば、混合気の当量比が所定値未満であってＳｏｏｔが殆ど生成されない場合においてＳｏｏｔ発生量の推定に係る計算を省略することができ、この結果、Ｓｏｏｔ発生量の推定精度を低下することなく計算負荷を少なくすることができる。
【００３９】
更には、前記すす発生量推定手段は、前記混合気の当量比が前記所定値以上であって、且つ、前記混合気の当量比と同混合気の温度との関係が所定の範囲内にある場合にのみ、Ｓｏｏｔ発生量の推定に係る計算を行うように構成されると好適である。より正確には、混合気の当量比と同混合気の温度との関係が所定の範囲外にある場合、Ｓｏｏｔが殆ど生成されないことが実験的・経験的に判っている。
【００４０】
従って、上記構成によれば、混合気の当量比が所定値未満である場合のみならず、混合気の当量比が同所定値以上であってＳｏｏｔが殆ど生成されない場合においてもＳｏｏｔ発生量の推定に係る計算を省略することができ、この結果、Ｓｏｏｔ発生量の推定精度を低下することなく計算負荷をより一層少なくすることができる。
【００４１】
上記本発明に係るすす発生量推定装置においては、前記すす発生量推定手段は、前記燃料の性状を考慮して前記すすの発生量を推定するように構成されることが好適である。前記燃料の性状としては、例えば、前記燃料に含まれる芳香族炭化水素の成分の含有率等が挙げられる。
【００４２】
燃料からＳｏｏｔが生成される反応の程度は燃料の性状に大きく依存し、例えば、燃料に含まれる芳香族炭化水素の成分の含有率が大きいほど燃料からＳｏｏｔが生成され易くなる。従って、上記構成によれば、Ｓｏｏｔ発生量をより一層精度良く推定することができる。
【発明を実施するための最良の形態】
【００４３】
以下、本発明による内燃機関（ディーゼル機関）のすす発生量推定装置の実施形態について図面を参照しつつ説明する。
【００４４】
図１は、本発明の実施形態に係る内燃機関のすす発生量推定装置を４気筒内燃機関（ディーゼル機関）１０に適用したシステム全体の概略構成を示している。このシステムは、燃料供給系統を含むエンジン本体２０、エンジン本体２０の各気筒の燃焼室（筒内）にガスを導入するための吸気系統３０、エンジン本体２０からの排ガスを放出するための排気系統４０、排気還流を行うためのＥＧＲ装置５０、及び電気制御装置６０を含んでいる。
【００４５】
エンジン本体２０の各気筒の上部には燃料噴射弁（噴射弁、インジェクタ）２１が配設されている。各燃料噴射弁２１は、図示しない燃料タンクと接続された燃料噴射用ポンプ２２に燃料配管２３を介して接続されている。燃料噴射用ポンプ２２は、電気制御装置６０と電気的に接続されていて、同電気制御装置６０からの駆動信号（後述する指令最終燃料噴射圧力Pcrfinに応じた指令信号）により燃料の実際の噴射圧力（吐出圧力）が同指令最終燃料噴射圧力Pcrfinになるように同燃料を昇圧するようになっている。
【００４６】
これにより、燃料噴射弁２１には、燃料噴射用ポンプ２２から前記指令最終燃料噴射圧力Pcrfinまで昇圧された燃料が供給されるようになっている。また、燃料噴射弁２１は、電気制御装置６０と電気的に接続されていて、同電気制御装置６０からの駆動信号（指令燃料噴射量（質量）Qfinに応じた指令信号）により噴射期間TAUだけ開弁し、これにより各気筒の燃焼室内に前記指令最終燃料噴射圧力Pcrfinにまで昇圧された燃料を前記指令燃料噴射量Qfinだけ直接噴射するようになっている。
【００４７】
吸気系統３０は、エンジン本体２０の各気筒の燃焼室にそれぞれ接続された吸気マニホールド３１、吸気マニホールド３１の上流側集合部に接続され同吸気マニホールド３１とともに吸気通路を構成する吸気管３２、吸気管３２内に回動可能に保持されたスロットル弁３３、電気制御装置６０からの駆動信号に応答してスロットル弁３３を回転駆動するスロットル弁アクチュエータ３３ａ、スロットル弁３３の上流において吸気管３２に順に介装されたインタクーラー３４と過給機３５のコンプレッサ３５ａ、及び吸気管３２の先端部に配設されたエアクリーナ３６とを含んでいる。
【００４８】
排気系統４０は、エンジン本体２０の各気筒にそれぞれ接続された排気マニホールド４１、排気マニホールド４１の下流側集合部に接続された排気管４２、排気管４２に配設された過給機３５のタービン３５ｂ、及び排気管４２に介装されたディーゼルパティキュレートフィルタ（以下、「ＤＰＮＲ」と称呼する。）４３を含んでいる。排気マニホールド４１及び排気管４２は排気通路を構成している。
【００４９】
ＥＧＲ装置５０は、排気ガスを還流させる通路（ＥＧＲ通路）を構成する排気還流管５１と、排気還流管５１に介装されたＥＧＲ制御弁５２と、ＥＧＲクーラー５３とを備えている。排気還流管５１はタービン３５ｂの上流側排気通路（排気マニホールド４１）とスロットル弁３３の下流側吸気通路（吸気マニホールド３１）を連通している。ＥＧＲ制御弁５２は電気制御装置６０からの駆動信号に応答し、再循環される排気ガス量（排気還流量、ＥＧＲガス流量）を変更し得るようになっている。
【００５０】
電気制御装置６０は、互いにバスで接続されたＣＰＵ６１、ＣＰＵ６１が実行するプログラム、テーブル（ルックアップテーブル、マップ）、及び定数等を予め記憶したＲＯＭ６２、ＣＰＵ６１が必要に応じてデータを一時的に格納するＲＡＭ６３、電源が投入された状態でデータを格納するとともに同格納したデータを電源が遮断されている間も保持するバックアップＲＡＭ６４、並びにＡＤコンバータを含むインターフェース６５等からなるマイクロコンピュータである。
【００５１】
インターフェース６５は、吸気管３２に配置された熱線式エアフローメータ７１、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に設けられた吸気温センサ７２、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気管圧力センサ７３、クランクポジションセンサ７４、アクセル開度センサ７５、燃料噴射用ポンプ２２の吐出口の近傍の燃料配管２３に配設された燃料温度センサ７６、及び、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気酸素濃度センサ７７と接続されていて、これらのセンサからの信号をＣＰＵ６１に供給するようになっている。
【００５２】
また、インターフェース６５は、燃料噴射弁２１、燃料噴射用ポンプ２２、スロットル弁アクチュエータ３３ａ、及びＥＧＲ制御弁５２と接続されていて、ＣＰＵ６１の指示に応じてこれらに駆動信号を送出するようになっている。
【００５３】
熱線式エアフローメータ７１は、吸気通路内を通過する吸入空気の質量流量（単位時間当りの吸入空気量、単位時間あたりの新気量）を計測し、同質量流量Ga（空気流量Ga）を表す信号を発生するようになっている。吸気温センサ７２は、エンジン１０のシリンダ（即ち、燃焼室、筒内）に吸入されるガスの温度（即ち、吸気温度）を検出し、同吸気温度Tbを表す信号を発生するようになっている。吸気管圧力センサ７３は、エンジン１０のシリンダに吸入されるガスの圧力（即ち、吸気管圧力）を検出し、同吸気管圧力Pbを表す信号を発生するようになっている。
【００５４】
クランクポジションセンサ７４は、各気筒の絶対クランク角度を検出し、実クランク角度CAactを表すとともにエンジン１０の回転速度であるエンジン回転速度NEをも表す信号を発生するようになっている。アクセル開度センサ７５は、アクセルペダルＡＰの操作量を検出し、アクセル操作量Accpを表す信号を発生するようになっている。燃料温度センサ７６は、燃料配管２３を通過する燃料の温度を検出し、燃料温度Tcrを表す信号を発生するようになっている。吸気酸素濃度センサ７７は、吸気中の酸素濃度を検出し、吸気酸素濃度RO2inを表す信号を発生するようになっている。
【００５５】
（Ｓｏｏｔ発生量の推定方法の概要）
次に、上記のように構成されたすす発生量推定装置（以下、「本装置」と云う。）によるＳｏｏｔ発生量の推定方法について説明する。
【００５６】
図２は、或る一つの気筒のシリンダ内(筒内、燃焼室内)に吸気マニホールド３１からガスが吸入され、燃焼室内に吸入されたガスが排気マニホールド４１へ排出される様子を模式的に示した図である。図２に示したように、燃焼室内に吸入されるガス（従って、筒内ガス）には、吸気管３２の先端部からスロットル弁３３を介して吸入された新気と、排気還流管５１からＥＧＲ制御弁５２を介して吸入されたＥＧＲガスが含まれる。
【００５７】
吸入される新気量（質量）と吸入されるＥＧＲガス量（質量）の和に対するＥＧＲガス量の割合（即ち、ＥＧＲ率）は、運転状態に応じて電気制御装置６０（ＣＰＵ６１）により適宜制御されるスロットル弁３３の開度、及びＥＧＲ制御弁５２の開度に応じて変化する。
【００５８】
かかる新気、及びＥＧＲガスは、吸気行程において開弁している吸気弁Vinを介してピストンの下降に伴って燃焼室内に吸入されて筒内ガスとなる。筒内ガスは、ピストンが圧縮下死点に達する時点近傍で吸気弁Vinが閉弁することにより燃焼室内に密閉され、その後の圧縮行程においてピストンの上昇に伴って圧縮される。
【００５９】
そして、ピストンが圧縮上死点近傍に達すると（具体的には、後述する燃料噴射開始時期（クランク角度）CAinjが到来すると）、本装置は、前記指令燃料噴射量Qfinに応じた噴射期間TAUだけ燃料噴射弁２１を開弁することで燃料を燃焼室内に直接噴射する。この結果、燃料噴射弁２１の噴孔から噴射された（液体の）燃料は、圧縮により高温になっている筒内ガスから受ける熱により直ちに燃料蒸気になるとともに、時間の経過に伴って同筒内ガスを取り込みながら混合気となって燃焼室内において円錐状に拡散していく。
【００６０】
ここで、指令燃料噴射量Qfinの燃料は、実際には、上記燃料噴射開始時期CAinjから噴射期間TAUだけ連続して噴射されるが、以下においては便宜上、指令燃料噴射量Qfinの燃料が燃料噴射開始時期CAinjにて一時（瞬時）に噴射されるものとして説明を続ける。
【００６１】
図３（ａ）は、燃料噴射弁２１の噴孔から、指令燃料噴射量（質量）Qfinの燃料が一時に噴射された時点（即ち、噴射後経過時間ｔ＝０）での質量Qfinの燃料（燃料蒸気）の様子を模式的に示した図である。図３（ｂ）は、その後の或る時点（任意の噴射後経過時間ｔ）での図３（ａ）に示した質量Qfinの燃料蒸気の様子を模式的に示した図である。
【００６２】
図３（ｂ）に示すように、質量Qfinの燃料蒸気は、燃料噴射開始時期CAinj（即ち、噴射後経過時間ｔ＝０）において噴射された後、噴霧角θをもって円錐状に拡散しながら筒内ガスを順次取り込んでいく。そして、質量Qfinの燃料蒸気は、任意の噴射後経過時間ｔにおいて、噴射後経過時間ｔの関数である質量Gの筒内ガス（以下、「混合気形成筒内ガス」と云うこともある。）と混ざり合って質量（Qfin＋G）の混合気となっているものと仮定する。
【００６３】
本装置は、係る混合気内において発生するＳｏｏｔの発生速度（具体的には、混合気内のＳｏｏｔ濃度の変化速度）を噴射後経過時間ｔ＝０から微小時間Δｔ毎に求め、それぞれの値を時間で積分（積算）していくことで任意の噴射後経過時間ｔにおける混合気内のＳｏｏｔ濃度（Ｓｏｏｔ発生量）を推定する。以下、先ず、係る推定に用いられる反応モデルについて説明する。
【００６４】
<前駆物質を経由した反応モデルを用いたＳｏｏｔ発生量の推定>
近年、Ｓｏｏｔの生成メカニズムについての研究によれば、燃料からＳｏｏｔの前駆物質Ｐ（precursor）を経由してＳｏｏｔが発生する反応過程が明らかになりつつある。ここで、本例では、前駆物質Ｐは、ＰＡＨ、芳香族中間生成物、飽和重合体、及びＣ２Ｈ２等、凝集により後にＳｏｏｔに成長し得る全ての単量体と定義され、Ｓｏｏｔは、単量体である前駆物質Ｐの２分子以上が凝集してなる重合体と定義される。
【００６５】
図４は、前駆物質Ｐの例としてＰＡＨ（多環芳香族炭化水素類）を挙げた場合における、燃料から前駆物質Ｐを経由してＳｏｏｔが発生する反応過程を模式的に示した図である。図４に示すように、この反応過程では、燃料分子は、先ず、熱分解反応等によって単量体である前駆物質Ｐとなる。
【００６６】
その後、前駆物質Ｐの粒子は、多環化反応等により成長していく。そして、前駆物質Ｐが或る分子量に達すると、前駆物質Ｐ同士の分子間力が大きくなって凝集（具体的には、クラスター化、鎖状化等）が発生し、これにより、重合体であるＳｏｏｔが生成される。なお、一般には、上記成長していく前駆物質Ｐの分子に含まれる炭素原子数が１００個程度になると、上記凝集が発生する傾向がある。
【００６７】
そこで、本装置では、上述したＳｏｏｔ濃度の推定に際し、燃料（Ｆｕｅｌ）からＳｏｏｔが発生する反応過程において前駆物質Ｐが経由される図５に示す反応モデルが用いられる。以下、係る反応モデルについて説明する。なお、以下、「濃度」とは「質量濃度」を意味するものとする。
【００６８】
図５に示すように、この反応モデルでは、燃料からの前駆物質Ｐの生成反応、生成された前駆物質Ｐの熱分解反応、生成された前駆物質ＰからのＳｏｏｔの生成反応、及び生成されたＳｏｏｔの酸化反応が考慮される。ここで、燃料からの前駆物質Ｐの生成反応は、混合気内の前駆物質Ｐ濃度[Ｐ]mixの変化速度（以下、「前駆物質Ｐの発生速度d[Ｐ]mix/dt」と称呼する。）を増大させる方向に働く。一方、生成された前駆物質Ｐの熱分解反応、及び生成された前駆物質ＰからのＳｏｏｔ生成反応は、前駆物質Ｐの発生速度d[Ｐ]mix/dtを減少させる方向に働く。従って、この反応モデルでは、前駆物質Ｐの発生速度d[Ｐ]mix/dtは、下記(1)式に従って表すことができる。
【００６９】
【数１】
0004603951.tif 000002
【００７０】
上記(1)式において、dmpf/dtは、蒸気燃料からの前駆物質Ｐの生成反応に基づく混合気内の前駆物質Ｐ濃度[Ｐ]mixの増加速度（以下、「前駆物質Ｐの生成速度」と称呼する。）であり、dmpd/dtは、上記生成された前駆物質Ｐの熱分解反応に基づく混合気内の前駆物質Ｐ濃度[Ｐ]mixの減少速度（以下、「前駆物質Ｐの分解速度」と称呼する。）であり、dmsf/dtは、上記生成された前駆物質ＰからのＳｏｏｔの生成反応に基づく混合気内のＳｏｏｔ濃度[Soot]mixの増加速度（換言すれば、混合気内の前駆物質Ｐ濃度[Ｐ]mixの減少速度。以下、「Ｓｏｏｔの生成速度」と称呼する。）である。
【００７１】
他方、生成された前駆物質ＰからのＳｏｏｔの生成反応は、混合気内のＳｏｏｔ濃度[Soot]mixの変化速度（以下、「Ｓｏｏｔの発生速度d[Soot]mix/dt」と称呼する。）を増大させる方向に働く。一方、生成されたＳｏｏｔの酸化反応は、Ｓｏｏｔの発生速度d[Soot]mix/dtを減少させる方向に働く。従って、この反応モデルでは、Ｓｏｏｔの発生速度d[Soot]mix/dtは、下記(2)式に従って表すことができる。
【００７２】
【数２】
0004603951.tif 000003
【００７３】
上記(2)式において、dmsf/dtは、上述したＳｏｏｔの生成速度であり、dmso/dtは、上記生成されたＳｏｏｔの酸化反応に基づく混合気内のＳｏｏｔ濃度[Soot]mixの減少速度（以下、「Ｓｏｏｔの酸化速度」と称呼する。）である。以下、上記(1)式の右辺各項、及び(2)式の右辺各項について順に説明する。
【００７４】
<<前駆物質Ｐの生成速度の取得>>
上記(1)式の右辺第１項の前駆物質Ｐの生成速度dmpf/dtは、気相反応についての反応速度を算出する場合に使用され得る所謂アレニウス（Arrhenius）の式に基づく下記(3)式に従って求めることができる。燃料からの前駆物質Ｐの生成反応は、気相反応とみなせるからである。
【００７５】
【数３】
0004603951.tif 000004
【００７６】
上記(3)式において、Apf，αpfは定数、Epfは活性化エネルギ（定数）、Ｒは気体定数（定数。以下も同じ。）である。また、[Fuel]mixは混合気内の燃料濃度（以下も同じ）、Pgは筒内ガスの圧力（以下も同じ）、Tmixは混合気の温度（以下も同じ）であり、これらは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、前駆物質Ｐの生成速度dmpf/dtも、上記(3)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【００７７】
<<前駆物質Ｐの分解速度の取得>>
上記(1)式の右辺第２項の前駆物質Ｐの分解速度dmpd/dtも、上記(3)式と同様、アレニウスの式に基づく下記(4)式に従って求めることができる。上記生成された前駆物質Ｐの熱分解反応も、気相反応とみなせるからである。
【００７８】
【数４】
0004603951.tif 000005
【００７９】
上記(4)式において、Apd，αpdは定数、Epdは活性化エネルギ（定数）である。また、[Ｐ]mixは混合気内の前駆物質Ｐ濃度であり、上記(1)式に従って噴射後経過時間ｔ＝０から微小時間Δｔ毎に得られる前駆物質Ｐの発生速度d[Ｐ]mix/dtを時間で積分（積算）していくことで、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、前駆物質Ｐの分解速度dmpd/dtも、上記(4)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【００８０】
<<Ｓｏｏｔの生成速度の取得>>
上記(1)式の右辺第３項のＳｏｏｔの生成速度dmsf/dtは、下記(5)式に従って表すことができる。下記(5)式において、dmsfpp/dtは前駆物質Ｐ同士の衝突に起因するＳｏｏｔの生成反応に基づくＳｏｏｔの生成速度（以下、「第１Ｓｏｏｔ生成速度」と称呼する。）であり、dmsfps/dtは前駆物質ＰとＳｏｏｔの衝突に起因するＳｏｏｔの生成反応に基づくＳｏｏｔの生成速度（以下、「第２Ｓｏｏｔ生成速度」と称呼する。）である。
【００８１】
【数５】
0004603951.tif 000006
【００８２】
このように、この反応モデルでは、前駆物質からＳｏｏｔが生成される反応過程において、前駆物質Ｐ同士の衝突によるＳｏｏｔ生成反応と、Ｓｏｏｔが生成された後における前駆物質ＰとＳｏｏｔとの衝突によるＳｏｏｔ生成反応とが別個独立に考慮されてＳｏｏｔの生成速度dmsf/dtが算出される。
【００８３】
前駆物質ＰからＳｏｏｔが生成される反応は、分子量の比較的大きい粒子に関わる反応である。従って、前駆物質ＰからＳｏｏｔが生成される反応を気相反応とみなすことは困難である。この反応モデルでは、上記(5)式の右辺第１項の第１Ｓｏｏｔ生成速度dmsfpp/dt、及び上記(5)式の右辺第２項の第２Ｓｏｏｔ生成速度dmsfps/dtは、以下のように求める。
【００８４】
一般に、粒子ａと粒子ｂの衝突による凝集反応の速度は、ガス中（例えば、混合気中）における粒子ａと粒子ｂの間の「衝突数」（ガス単位質量あたりの衝突数）と、衝突した粒子ａと粒子ｂが凝集する確率（以下、「凝集確率」と称呼する。）の積で表すことができる。
【００８５】
また、一般に、ガス中における粒子ａと粒子ｂの間の「衝突数」は、粒子ａと粒子ｂの間の「相対速度」と、粒子ａと粒子ｂの衝突に関わる面積である「衝突断面積」と、粒子ａの「粒子数密度」（ガス単位質量あたりの粒子ａの粒子数）と、粒子ｂの「粒子数密度」（ガス単位質量あたりの粒子ｂの粒子数）の積で表すことができる。ここで、「衝突断面積」は、粒子ａの半径と粒子ｂの半径の和を半径とする円の面積であり、下記(6)式で表すことができる。ここで、Da，Dbはそれぞれ、粒子ａ、粒子ｂの直径（以下も同じ）である。
【００８６】
【数６】
0004603951.tif 000007
【００８７】
また、「粒子数密度」は、ガス中における粒子の濃度を１粒子の質量で除した値である。従って、粒子ａの「粒子数密度」は「Ｍa/ｍa」で、粒子ｂの「粒子数密度」は「Ｍb/ｍb」で表すことができる。ここで、Ｍa，Ｍbはそれぞれ、粒子ａ、粒子ｂのガス中における濃度であり、ｍa，ｍbはそれぞれ、粒子ａ、粒子ｂの１粒子の質量（以下も同じ）である。また、ρa，ρbをそれぞれ、粒子ａ、粒子ｂの１粒子の密度とすると、粒子ａの「粒子数密度」は「(6・Ｍa)/(ρa・π・Da3)」と、粒子ｂの「粒子数密度」は「(6・Ｍb)/(ρb・π・Db3)」と表すこともできる。
【００８８】
また、粒子ａと粒子ｂの間の「相対速度」は、マクスウエル−ボルツマン（Maxwell−Boltzmann）の分布関数を利用すれば、下記(7)式で表すことができる。ここで、μは粒子ａと粒子ｂの換算質量であって、μ＝(ｍa・ｍb)/(ｍa＋ｍb)で表される。ｋはボルツマン定数であり、Ｔはガスの温度である。
【００８９】
【数７】
0004603951.tif 000008
【００９０】
従って、ガス中における粒子ａと粒子ｂの間の「衝突数」は、下記(8)式で表すことができる。
【００９１】
【数８】
0004603951.tif 000009
【００９２】
以上のことから、衝突した粒子ａと粒子ｂが凝集する上記凝集確率をＡ１とし、下記(9)式に示すように係数Ａ２を定義すると、粒子ａと粒子ｂの衝突による凝集反応の速度は、一般に、下記(10)式に従って表すことができる。
【００９３】
【数９】
0004603951.tif 000010
【００９４】
【数１０】
0004603951.tif 000011
【００９５】
ここで、上記第１Ｓｏｏｔ生成速度dmsfpp/dtは、前駆物質Ｐ同士の衝突による凝集反応の速度であり、上記第２Ｓｏｏｔ生成速度dmsfps/dtは、前駆物質ＰとＳｏｏｔとの衝突による凝集反応の速度である。従って、この反応モデルでは、上記第１Ｓｏｏｔ生成速度dmsfpp/dt、及び第２Ｓｏｏｔ生成速度dmsfps/dtは、上記(10)式を用いて求められる。
【００９６】
先ず、上記第２Ｓｏｏｔ生成速度dmsfps/dtについて説明する。この場合、上記粒子ａが前駆物質Ｐの粒子に、上記粒子ｂがＳｏｏｔの粒子に対応すると考える。そうすると、上記第２Ｓｏｏｔ生成速度dmsfps/dtは、上記(10)式に対応する下記(11)式に従って求めることができる。ここで、Asfps1は、衝突した前駆物質ＰとＳｏｏｔとが凝集する凝集確率であり、Asfps2は、下記(12)式に示すように定義される、上記係数Ａ２に対応する係数である。
【００９７】
【数１１】
0004603951.tif 000012
【００９８】
【数１２】
0004603951.tif 000013
【００９９】
上記(12)式において、μpsは、前駆物質Ｐの１粒子とＳｏｏｔの１粒子の換算質量であって、μps＝(ｍp・ｍsoot)/(ｍp＋ｍsoot)で表される。ｍp，ｍsootはそれぞれ、前駆物質Ｐの１粒子の質量、Ｓｏｏｔの１粒子の質量（以下も同じ）である。Dp，Dsootはそれぞれ、前駆物質Ｐの１粒子の直径、Ｓｏｏｔの１粒子の直径（以下も同じ）である。ρp，ρsootはそれぞれ、前駆物質Ｐの１粒子の密度、Ｓｏｏｔの１粒子の密度（以下も同じ）である。
【０１００】
ここで、この反応モデルでは、ｍp，ｍsoot，Dp，Dsoot，ρp，ρsootは全て一定値として扱う。従って、上記(12)式から理解できるように、上記係数Asfps2は一定値となる。また、上記凝集確率Asfps1も一定であると仮定する。また、このモデルでは、前駆物質Ｐ及びＳｏｏｔ以外の他の粒子の衝突による上記第２Ｓｏｏｔ生成速度への影響も考慮するため、上記(11)式の右辺の因数として筒内ガス圧力Pgに依存する項「Pgαsfps」が加えられる（αsfpsは定数）。即ち、このモデルでは、上記第２Ｓｏｏｔ生成速度dmsfps/dtは、上記係数Asfps2と上記凝集確率Asfps1の積を係数Asfps（定数）と置くと、下記(13)式に従って求められる。
【０１０１】
【数１３】
0004603951.tif 000014
【０１０２】
上記(13)式において、上述したように、混合気内の前駆物質Ｐ濃度[Ｐ]mix、筒内ガス圧力Pg、混合気温度Tmixは、上述したように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。また、混合気内のＳｏｏｔ濃度[Soot]mixも、後述するように噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、第２Ｓｏｏｔ生成速度dmsfps/dtも、上記(13)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【０１０３】
このように、この反応モデルでは、第２Ｓｏｏｔ生成速度dmsfps/dtは、前駆物質Ｐの粒子とＳｏｏｔの粒子の衝突数、及び衝突した前駆物質Ｐの粒子とＳｏｏｔの粒子とが凝集する上記凝集穫率Asfps1を考慮して算出される。換言すれば、第２Ｓｏｏｔ生成速度dmsfps/dtは、前駆物質Ｐの粒子とＳｏｏｔの粒子との衝突に起因して新たなＳｏｏｔが生成される確率を考慮して算出される。加えて、上記「前駆物質Ｐの粒子とＳｏｏｔの粒子の衝突数」は、混合気中の前駆物質Ｐ濃度[Ｐ]mixと、混合気中のＳｏｏｔ濃度[Soot]mixと、混合気温度Tmixとに基づいて算出される。
【０１０４】
次に、上記第１Ｓｏｏｔ生成速度dmsfpp/dtについて説明する。この場合、上記粒子ａ、上記粒子ｂ共に、前駆物質Ｐの粒子に対応すると考える。そうすると、換算質量μ＝mp/2＝(ρp・π・Dp3)/6・(1/2)と表すことができる。また、この場合、同じ種類の粒子の衝突が扱われていることに起因して、一回の衝突現象を２回数えることになる。即ち、前駆物質Ｐ同士の「衝突数」は、上記(8)式から得られる値に(1/2)を乗じた値になる。
【０１０５】
以上のことを考慮すると、上記第１Ｓｏｏｔ生成速度dmsfpp/dtは、上記(10)式に対応する下記(14)式に従って求めることができる。ここで、Asfpp1は、衝突した前駆物質Ｐ同士が凝集する凝集確率であり、Asfpp2は、下記(15)式に示すように定義される、上記係数Ａ２に対応する係数である。
【０１０６】
【数１４】
0004603951.tif 000015
【０１０７】
【数１５】
0004603951.tif 000016
【０１０８】
ここで、この反応モデルでは、上述したように、Dp，ρpは一定値として扱う。従って、上記(15)式から理解できるように、上記係数Asfpp2は一定値となる。また、上記凝集確率Asfpp1も一定であると仮定する。また、このモデルでは、前駆物質Ｐ以外の他の粒子の衝突による上記第１Ｓｏｏｔ生成速度への影響も考慮するため、上記(14)式の右辺の因数として筒内ガス圧力Pgに依存する項「Pgαsfpp」が加えられる（αsfppは定数）。即ち、このモデルでは、上記第１Ｓｏｏｔ生成速度dmsfpp/dtは、上記係数Asfpp2と上記凝集確率Asfpp1の積を係数Asfpp（定数）と置くと、下記(16)式に従って求められる。
【０１０９】
【数１６】
0004603951.tif 000017
【０１１０】
上記(16)式において、上述したように、混合気内の前駆物質Ｐ濃度[Ｐ]mix、筒内ガス圧力Pg、混合気温度Tmixは、上述したように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、第１Ｓｏｏｔ生成速度dmsfpp/dtも、上記(16)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【０１１１】
このように、この反応モデルでは、第１Ｓｏｏｔ生成速度dmsfpp/dtは、前駆物質Ｐの粒子同士の衝突数、及び衝突した前駆物質Ｐの粒子同士が凝集する上記凝集穫率Asfpp1を考慮して算出される。換言すれば、第１Ｓｏｏｔ生成速度dmsfpp/dtは、前駆物質Ｐの粒子同士の衝突に起因してＳｏｏｔが生成される確率を考慮して算出される。加えて、上記「前駆物質Ｐの粒子同士の衝突数」は、混合気中の前駆物質Ｐ濃度[Ｐ]mixと、混合気温度Tmixとに基づいて算出される。
【０１１２】
ここで、上述した「衝突した前駆物質Ｐの粒子とＳｏｏｔの粒子とが凝集する凝集穫率Asfps1」は、上述した「衝突した前駆物質Ｐの粒子同士が凝集する凝集穫率Asfpp1」よりも大きい。これは、重合体であるＳｏｏｔの分子の質量が単量体である前駆物質の分子の質量よりも大きいことに起因して、前駆物質Ｐの分子とＳｏｏｔの分子の間の分子間力が前駆物質Ｐ同士の分子間力よりも大きいことに基づく。
【０１１３】
このように、第１Ｓｏｏｔ生成速度dmsfpp/dt、及び第２Ｓｏｏｔ生成速度dmsfps/dtが、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができるから、上記(1)式の右辺第３項のＳｏｏｔ生成速度dmsf/dtも、上記(5)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【０１１４】
<<前駆物質Ｐ濃度の取得>>
以上より、上記(1)式の右辺の項の全てを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができるから、上記(1)式に従って、前駆物質Ｐの発生速度d[Ｐ]mix/dtを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。従って、係る微小時間Δｔ毎のそれぞれの前駆物質Ｐの発生速度d[Ｐ]mix/dtの値を時間で積分（積算）していくことで、任意の噴射後経過時間ｔにおける混合気内の前駆物質Ｐ濃度[Ｐ]mixを推定することができる。
【０１１５】
次に、上記(2)式の右辺各項について順に説明する。上記(2)式の右辺第１項のＳｏｏｔ生成速度dmsf/dtは、上述した上記(1)式の右辺第３項のＳｏｏｔ生成速度dmsf/dtと同じであり、上記(5)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【０１１６】
<<Ｓｏｏｔの酸化速度の取得>>
上記(2)式の右辺第２項のＳｏｏｔの酸化速度dmso/dtは、この反応モデルでは、日本機械学会論文集(Ｂ編)  48巻432号「直接噴射式ディーゼル機関の燃焼モデルと性能予測」（以下、「非特許文献１」と称呼する。）にて紹介された実験式である下記(17)式に従って求める。
【０１１７】
【数１７】
0004603951.tif 000018
【０１１８】
上記(17)式において、Aso，αsoは定数、Esoは活性化エネルギ（定数）である。また、[O2]mixは混合気内の酸素濃度（以下も同じ）であり、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、Ｓｏｏｔの酸化速度dmso/dtも、上記(17)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。
【０１１９】
<<Ｓｏｏｔ濃度の取得>>
以上より、上記(2)式の右辺の項の全てを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができるから、上記(2)式に従って、Ｓｏｏｔの発生速度d[Soot]mix/dtを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。従って、係る微小時間Δｔ毎のそれぞれのＳｏｏｔの発生速度d[Soot]mix/dtの値を時間で積分（積算）していくことで、任意の噴射後経過時間ｔにおける混合気内のＳｏｏｔ濃度[Soot]mixを推定することができる。
【０１２０】
更には、このようにして、噴射後経過時間ｔにおける混合気内のＳｏｏｔ濃度[Soot]mixを推定することができれば、この混合気内のＳｏｏｔ濃度[Soot]mixの値に混合気質量Mmixの値を乗じることで噴射後経過時間ｔにおけるＳｏｏｔ発生質量（以下、「Ｓｏｏｔ発生質量Ｓｏｏｔnew」と称呼する。）を求めることができる。
【０１２１】
以上のように、本装置が用いる反応モデルでは、前駆物質Ｐの発生速度d[Ｐ]mix/dt（従って、混合気内の前駆物質Ｐ濃度[Ｐ]mix。前駆物質の発生量）に基づいてＳｏｏｔの発生速度d[Soot]mix/dtが算出され、同算出されたＳｏｏｔの発生速度d[Soot]mix/dtを時間で積分していくことで混合気内のＳｏｏｔ濃度[Soot]mix（Ｓｏｏｔの発生量）が推定される。
【０１２２】
このように、燃料からＳｏｏｔが発生する反応過程を、燃料から前駆物質Ｐが発生する反応過程と、前駆物質ＰからＳｏｏｔが発生する反応過程の２段階とすることで以下の効果を得ることができる。
【０１２３】
第１に、燃料からの前駆物質Ｐの発生により、先ず、前駆物質Ｐ濃度[Ｐ]mixが増大していく。そして、前駆物質Ｐ濃度[Ｐ]mixが増大していくことでＳｏｏｔ濃度[Soot]mixが増大していく。これにより、燃料からＳｏｏｔが発生する反応過程において実験的・経験的に知られている「Ｓｏｏｔ発生（反応）の遅れ」を精度良く模擬することができる。
【０１２４】
第２に、上記(1)式に従って、前駆物質Ｐの発生速度d[Ｐ]mix/dtが算出される際、上記(4)に従って求められる前駆物質Ｐの分解速度dmpd/dtが考慮されている。この前駆物質Ｐの分解速度dmpd/dtは混合気温度Tmixが高くなるほど大きくなる。従って、混合気温度Tmixが高くなるほど前駆物質Ｐの発生速度d[Ｐ]mix/dt（従って、前駆物質Ｐ濃度[Ｐ]mix）が減少していく。これにより、実験的・経験的に知られている「混合気温度が或る温度を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」をも精度良く模擬することができる。
【０１２５】
以上、前駆物質を経由した反応モデルを用いたＳｏｏｔ発生量（混合気内のＳｏｏｔ濃度[Soot]mix）の推定について説明した。次に、係るＳｏｏｔ濃度[Soot]mixの推定に必要となる、混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、混合気内の酸素濃度[O2]mix、及び筒内ガス圧力Pgの推定方法について説明する。以下、先ず、これらの物理量の推定に必要となる、任意の噴射後経過時間ｔにおける空気過剰率λの取得方法について説明する。
【０１２６】
<空気過剰率λの取得>
噴射後経過時間ｔにおける空気過剰率λは、下記(18)式に示すように定義される。下記(18)式において、stoichは、単位質量の燃料の燃焼に必要な筒内ガスの質量（以下、「筒内ガス理論空燃比stoich」と呼ぶ。）である。筒内ガス理論空燃比stoichの値は、吸気中の酸素濃度に応じて変化すると考えられるから、上記吸気酸素濃度RO2inを引数とする所定の関数に従って取得され得る。また、上述したように、Gは混合気形成筒内ガスの質量であり、Qfinは指令燃料噴射量（質量）である（図３を参照）。
【０１２７】
【数１８】
0004603951.tif 000019
【０１２８】
このように定義される空気過剰率λは、例えば、日本機械学会論文集  25-156(1959年)，820ページ  「ディーゼル機関の噴霧到達距離に関する研究」  和栗雄太郎，藤井勝，網谷竜夫，恒屋礼次郎  （以下、「非特許文献２」と称呼する。）にて紹介された実験式である下記(19)式、及び下記(20)式に基づいて噴射後経過時間ｔの関数として求めることができる。
【０１２９】
【数１９】
0004603951.tif 000020
【０１３０】
【数２０】
0004603951.tif 000021
【０１３１】
上記(20)式において、ｔは上記噴射後経過時間であり、dλ/dtは噴射後経過時間ｔの関数である燃料希釈率である。また、ｃは収縮係数、ｄは燃料噴射弁２１の噴孔径、ρfは（液体の）燃料密度、Ｌは論理希釈ガス量であって、これらの各値は全て定数である。
【０１３２】
上記(20)式において、ΔPは有効噴射圧力であって、上記最終燃料噴射圧力Pcrfinから噴射開始時点（即ち、噴射後経過時間ｔ＝０）での筒内ガス圧力Pg0を減じた値である。筒内ガス圧力Pg0は、圧縮行程（及び膨張行程）における筒内ガスの状態が吸気弁Vinの閉弁時（即ち、筒内ガスが密閉された時点。以下、「ＩＶＣ」と呼ぶ。）以降、断熱変化するとの仮定のもと、下記(21)式に従って求めることができる。
【０１３３】
【数２１】
0004603951.tif 000022
【０１３４】
上記(21)式において、PgivcはＩＶＣにおける筒内ガス圧力である。上述したように、ＩＶＣは圧縮下死点近傍であるから、ＩＶＣにおいて筒内ガス圧力は吸気管圧力Pbと略等しいと考えられる。従って、ＩＶＣにおいて吸気管圧力センサ７３により検出される吸気管圧力PbがPgivcとして使用され得る。Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積であり、Vg(CAinj)は噴射後経過時間ｔ＝０におけるクランク角度CAに対応する筒内容積である。筒内容積Vgは機関１０の設計諸元に基づいてクランク角度CAの関数Vg(CA)として取得することができるから、Vg(Caivc)，Vg(CAinj)も取得することができる。κは筒内ガスの比熱比（本例では、一定）である。
【０１３５】
また、上記(20)式において、θは図３（ｂ）に示した噴霧角である。噴霧角θは、噴射開始時点（即ち、噴射後経過時間ｔ＝０）における筒内ガスの密度ρg0、及び上記有効噴射圧力ΔPに応じて変化すると考えられるから、筒内ガスの密度ρg0、及び有効噴射圧力ΔPと噴霧角θとの関係を予め規定したテーブルMapθに基づいて取得することができる。筒内ガスの密度ρg0は、筒内ガスの全質量Mgを、噴射後経過時間ｔ＝０における上記筒内容積Vg(CAinj)で除することで取得することができる。筒内ガスの全質量Mgは、ＩＶＣにおける気体の状態方程式に基づく下記(22)式に従って取得され得る。下記(22)式において、TgivcはＩＶＣにおける筒内ガス温度である。ＩＶＣは圧縮下死点近傍であるから、ＩＶＣにおいて筒内ガス温度は吸気温度Tbと略等しいと考えられる。従って、ＩＶＣにおいて吸気温センサ７２により検出される吸気温度TbがTgivcとして使用され得る。Rは筒内ガスのガス定数（本例では、一定）である。
【０１３６】
【数２２】
0004603951.tif 000023
【０１３７】
また、上記(20)式において、ρgは噴射後経過時間ｔにおける筒内ガス密度であって、前記筒内ガスの全質量Mgを、噴射後経過時間ｔにおける上記筒内容積Vg(CA)で除することで、噴射後経過時間ｔの関数として取得することができる。
【０１３８】
以上のことから、有効噴射圧力ΔPと噴霧角θとを上述のようにして求めれば、噴射後経過時間ｔの値と同噴射後経過時間ｔの関数である筒内ガス密度ρgの値とを使用して、上記(20)式に従って燃料希釈率dλ/dtが噴射後経過時間ｔの関数として求められる。そして、噴射後経過時間ｔ＝０から微小時間Δｔ（例えば、０．１msec）毎に求めた燃料希釈率dλ/dtの値を上記(19)式に従って時間で積分（積算）していくことで噴射後経過時間ｔにおける空気過剰率λを噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。
【０１３９】
なお、上記(20)式から取得される燃料希釈率dλ/dtの値は常に正の値となることから上記(19)式から取得される空気過剰率λの値は噴射後経過時間ｔの増大に従って増加していく。そうすると、上記(18)式から理解できるように、混合気形成筒内ガスの質量Gが噴射後経過時間ｔの増大に従って増加していく。このことは、噴射後の燃料蒸気が円錐状に拡散していくことに伴って燃料蒸気と混ざり合う（燃料蒸気が取り込む）筒内ガス（従って、混合気形成筒内ガス）の量が増大していくことに対応している。
【０１４０】
<混合気温度Tmixの取得>
次に、上述のように取得される空気過剰率λの値を利用して任意の噴射後経過時間ｔにおける混合気温度Tmixを取得する方法について説明する。一般に、混合気の熱エネルギー（エンタルピ）Hmixは、混合気温度Tmixを用いて下記(23)式に従って表すことができる。
【０１４１】
【数２３】
0004603951.tif 000024
【０１４２】
上記(23)式において、Mmixは混合気の総質量（混合気質量）、Cmixは混合気の定圧比熱である。従って、混合気のエンタルピHmix、混合気質量Mmix、及び混合気の定圧比熱Cmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていく（更新していく）ことで、下記(24)式に従って混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。以下、先ず、混合気質量Mmixの求め方について説明する。
【０１４３】
【数２４】
0004603951.tif 000025
【０１４４】
<<混合気質量Mmix>>
上述したように、質量Qfinの燃料蒸気は、任意の噴射後経過時間ｔにおいて質量Gの混合気形成筒内ガスと混ざり合って質量（Qfin＋G）の混合気となっているから、任意の噴射後経過時間ｔにおける混合気質量Mmixは(Qfin＋G)である。ここで、上記(18)式より「G＝stoich・λ・Qfin」と表すことができるから、混合気質量Mmixは、空気過剰率λを用いて下記(25)式にて表すことができる。
【０１４５】
【数２５】
0004603951.tif 000026
【０１４６】
よって、上述したように噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得され得る空気過剰率λの値を上記(25)式に順次適用していくことで、混合気質量Mmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。
【０１４７】
<<混合気の定圧比熱Cmix>>
次に、混合気の定圧比熱Cmixの求め方について説明する。一般に、混合気の定圧比熱Cmixは、同混合気内の酸素濃度[O2]mix、及び混合気温度Tmixに大きく依存すると考えられる。ここで、混合気内の酸素濃度[O2]mixは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。従って、混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができれば、下記(26)式に従って、混合気の定圧比熱Cmixを微小時間Δｔ毎に求めることができる。
【０１４８】
【数２６】
0004603951.tif 000027
【０１４９】
上記(26)式において、funcCmixは、混合気の酸素濃度[O2]mix、及び混合気温度Tmixを引数とする混合気の定圧比熱Cmixを求めるための関数である。なお、上記(26)式を使用して混合気の定圧比熱Cmixを微小時間Δｔ毎に求めていく際における[O2]mix，Tmixの引数値としては、それぞれ現時点（即ち、噴射後経過時間ｔ）よりも微小時間Δｔ前の値が使用される。
【０１５０】
<<混合気のエンタルピHmix>>
次に、混合気のエンタルピHmixの求め方について説明する。いま、噴射後経過時間（ｔ−Δt）における混合気のエンタルピHmix(ｔ−Δｔ)が既知である場合において、同噴射後経過時間(ｔ−Δｔ)から噴射後経過時間ｔまでの微小時間Δｔの間における混合気のエンタルピの増加分ΔHmixについて考える。この混合気のエンタルピの増加分ΔHmixとしては、微小時間Δｔの間に混合気に新たに取り込まれる筒内ガスの熱エネルギーΔHgと、同微小時間Δｔの間において混合気内で発生する化学反応による反応熱Hrとが挙げられる。
【０１５１】
先ず、上記筒内ガスの熱エネルギーΔHgは、下記(27)式により表すことができる。ここにおいて、gは微小時間Δｔの間に混合気に新たに取り込まれる筒内ガスの質量である。この質量gは、噴射後経過時間ｔにおける混合気形成筒内ガスの質量から噴射後経過時間（ｔ−Δｔ）における混合気形成筒内ガスの質量を減じた値である。従って、上述した関係「G＝stoich・λ・Qfin」を利用して下記(28)式により求めることができる。(28)式において、λ(ｔ)，λ(ｔ−Δｔ)はそれぞれ、噴射後経過時間ｔ，(ｔ−Δｔ)における空気過剰率であり、上記(19)式、(20)式から求めることができる。
【０１５２】
【数２７】
0004603951.tif 000028
【０１５３】
【数２８】
0004603951.tif 000029
【０１５４】
また、上記(27)式において、Tgは、噴射後経過時間ｔにおける筒内ガスの温度であり、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(29)式に従って求めることができる。下記(29)式において、上述したように、TgivcはＩＶＣにおける筒内ガス温度であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。
【０１５５】
【数２９】
0004603951.tif 000030
【０１５６】
また、上記(27)式において、Cgは、噴射後経過時間ｔにおける筒内ガスの定圧比熱であり、混合気の定圧比熱Cmixを求める上記(26)式と同様、下記(30)式に従って求めることができる。下記(30)式において、funcCgは、吸気中の酸素濃度[O2]in、及び筒内ガス温度Tgを引数とする筒内ガスの定圧比熱Cgを求めるための関数である。
【０１５７】
【数３０】
0004603951.tif 000031
【０１５８】
なお、上記(30)式を使用して筒内ガスの定圧比熱Cgを微小時間Δｔ毎に求めていく際における[O2]inの引数値としては、吸気酸素濃度センサ７７により検出される上記吸気酸素濃度RO2inが使用される。また、上記筒内ガス温度Tgの引数値としては、現時点（即ち、噴射後経過時間ｔ）における値が使用される。以上により、上記(27)式の右辺の項の全てを求めることができるから、同(27)式に従って上記筒内ガスの熱エネルギーΔHgを求めることができる。
【０１５９】
次に、上記微小時間Δｔの間において混合気内で発生する化学反応による反応熱Hrは、下記(31)式で表すことができる。下記(31)式において、Hfは所定の定数であり、qrは上記微小時間Δｔの間において混合気内で発生する化学反応による燃料消費量である。
【０１６０】
【数３１】
0004603951.tif 000032
【０１６１】
上記燃料消費量qrの対象となる化学反応としては、前駆物質Ｐの生成反応、着火反応（熱炎反応）、低温酸化反応（冷炎反応）のみならず、その他の種々の化学反応が含まれる。この燃料消費量qrは、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、混合気温度Tmixに大きく依存すると考えられるから、下記(32)式に従って表すことができる。
【０１６２】
【数３２】
0004603951.tif 000033
【０１６３】
上記(32)式において、funcqrは、混合気内の酸素濃度[O2]mix、混合気内の燃料濃度[Fuel]mix、及び混合気温度Tmixを引数とする上記燃料消費量qrを求めるための関数である。混合気内の燃料濃度[Fuel]mixも、上記混合気内の酸素濃度[O2]mixと同様、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。上記(32)式を使用して上記燃料消費量qrを微小時間Δｔ毎に求めていく際における[O2]mix，[Fuel]mixの引数値としては、ぞれぞれ現時点（即ち、噴射後経過時間ｔ）よりも微小時間Δｔ前の値が使用される。
【０１６４】
また、上記(32)式の混合気温度Tmixの引数値としては化学反応前の混合気温度Tpreが使用される。この化学反応前の混合気温度Tpreは、上記(28)式にて算出される質量gの筒内ガスが混合気に新たに取り込まれた後であって、且つ、噴射後経過時間(ｔ−Δｔ)からの上記微小時間Δｔの間における化学反応が発生する前の段階における混合気温度であり、下記(33)式に従って求めることができる。
【０１６５】
【数３３】
0004603951.tif 000034
【０１６６】
上記(33)式において、Mmix(ｔ−Δｔ)，Cmix(ｔ−Δｔ)はそれぞれ、噴射後経過時間(ｔ−Δｔ)における混合気質量、及び混合気の定圧比熱であって、上記(25)式、及び上記(26)式によりそれぞれ取得され得る。更には、噴射後経過時間(ｔ−Δｔ)における混合気のエンタルピHmix(ｔ−Δｔ)は既知である。以上により、上記化学反応前の混合気温度Tpreを求めることができる。従って、上記(32)式の右辺の引数値を全て取得することができるから、(32)式，(31)式に従って上記化学反応による反応熱Hrを求めることができる。
【０１６７】
以上のことから、噴射後経過時間（ｔ−Δt）における混合気のエンタルピHmix(ｔ−Δｔ)が既知である場合において、同噴射後経過時間(ｔ−Δｔ)から噴射後経過時間ｔまでの微小時間Δｔの間における混合気のエンタルピの増加分ΔHmix(＝ΔHg＋Hr)が求められるから、噴射後経過時間ｔにおける混合気のエンタルピHmix(ｔ)（＝Hmix(ｔ−Δｔ)＋ΔHmix）を求めることができる。
【０１６８】
更には、噴射後経過時間ｔ＝０における混合気は、筒内ガスが取り込まれる前の状態（即ち、燃料蒸気のみ）であるから（図３（ａ）を参照）、この時点での混合気のエンタルピHmix(0)は、下記(34)式にて求めることができる。ここにおいて、Cfは燃料（蒸気）の定圧比熱（ここでは、定数）である。
【０１６９】
【数３４】
0004603951.tif 000035
【０１７０】
また、Tfは、燃料蒸気そのものの温度であり、液体の燃料が噴射直後に燃料蒸気に変化する際の単位質量当たりの潜熱Qvaporを考慮して下記(35)式に従って求めることができる。下記(35)式において、Tcrは噴射後経過時間ｔ＝０において燃料温度センサ７６により検出される液体の燃料温度である。αcrは燃料が燃料噴射用ポンプ２２の吐出口近傍から燃料噴射弁２１までの燃料配管２３を通過する際の熱損失分を考慮するための補正係数である。
【０１７１】
【数３５】
0004603951.tif 000036
【０１７２】
従って、噴射後経過時間ｔ＝０における混合気のエンタルピHmix(0)も求めることができる。以上より、混合気のエンタルピHmixを、噴射後経過時間ｔ＝０から微小時間Δｔ毎に取得することができる。
【０１７３】
このようにして、混合気のエンタルピHmix、混合気質量Mmix、及び混合気の定圧比熱Cmixがそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求められるから、上記(24)式に従って混合気温度Tmixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。
【０１７４】
<混合気内の燃料濃度[Fuel]mixの取得>
次に、混合気内の燃料濃度[Fuel]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の燃料濃度[Fuel]mixは、上記(25)式により取得される噴射後経過時間ｔにおける混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する燃料の質量」の割合である。
【０１７５】
ここで、「噴射後経過時間ｔにおいて混合気内に存在する燃料の質量」は、噴射後経過時間ｔ＝０にて噴射された燃料量（指令燃料噴射量Qfin）から、噴射から現時点（噴射後経過時間ｔ）までの間に化学反応で消費された燃料分を減じた値である。従って、噴射後経過時間ｔにおける混合気内の燃料濃度[Fuel]mixは、下記(36)式で表すことができる。
【０１７６】
【数３６】
0004603951.tif 000037
【０１７７】
上記(36)式において、「Σqr」は、噴射から現時点（噴射後経過時間ｔ）まで微小時間Δｔ毎に上記(32)式に従って逐次取得・更新されていくそれぞれの燃料消費量qrの和である。このように、燃料消費量qr、及び混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(36)式に従って混合気内の燃料濃度[Fuel]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。
【０１７８】
<混合気内の酸素濃度[O2]mixの取得>
次に、混合気内の酸素濃度（質量濃度）[O2]mixを取得する方法について説明する。噴射後経過時間ｔにおける混合気内の酸素濃度[O2]mixは、噴射後経過時間ｔにおける上記混合気質量Mmixに対する、「噴射後経過時間ｔにおいて混合気内に存在する酸素の質量」の割合である。
【０１７９】
ここで、「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」は、上述した噴射後経過時間ｔにおける混合気形成筒内ガスの質量Gから、噴射から現時点（噴射後経過時間ｔ）までの間に化学反応で消費された筒内ガス分を減じた値である。上記微小時間Δｔにおける燃料消費量qrの燃料と反応して同微小時間Δｔの間において混合気内で化学反応により消費される筒内ガスの消費量grは、下記(37)式にて表すことができる。
【０１８０】
【数３７】
0004603951.tif 000038
【０１８１】
従って、「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」は、「G−Σgr」と表すことができる。ここで、「Σgr」は、噴射から現時点（噴射後経過時間ｔ）まで微小時間Δｔ毎に上記(37)式に従って逐次取得・更新されていくそれぞれの筒内ガス消費量grの和である。
【０１８２】
「噴射後経過時間ｔにおいて混合気内に存在する酸素の質量」は、上記「噴射後経過時間ｔにおいて混合気内に存在する筒内ガスの質量」に、筒内ガス内の酸素濃度（従って、上記吸気中の酸素濃度[O2]in）を乗じることで求めることができる。以上のことから、噴射後経過時間ｔにおける混合気内の酸素濃度[O2]mixは、下記(38)式で表すことができる。
【０１８３】
【数３８】
0004603951.tif 000039
【０１８４】
このように、混合気形成筒内ガスの質量G、筒内ガス消費量gr、及び混合気質量Mmixをそれぞれ噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることで、上記(38)式に従って混合気内の酸素濃度[O2]mixを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めていくことができる。
【０１８５】
<筒内ガス圧力Pgの取得>
次に、筒内ガス圧力Pgを取得する方法について説明する。筒内ガス圧力Pgは、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(39)式に従って求めることができる。下記(39)式において、上述したように、PgivcはＩＶＣにおける筒内ガス圧力であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。以上、Ｓｏｏｔ濃度[Soot]mixの推定に必要となる、空気過剰率λ、混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、混合気内の酸素濃度[O2]mix、及び筒内ガス圧力Pgの推定方法について説明した。
【０１８６】
【数３９】
0004603951.tif 000040
【０１８７】
（燃料噴射制御の概要）
本装置は、上述した反応モデルを用いた混合気内のＳｏｏｔ濃度[Soot]mixの推定に係わる計算を筒内ガスの量が確定するＩＶＣの直後に開始し、燃料噴射開始時期CAinjの到来前に上記Ｓｏｏｔ発生質量Ｓｏｏｔnewの推定を完了する。
【０１８８】
そして、本装置は、目標Ｓｏｏｔ発生質量Ｓｏｏｔterを機関の運転状態から求め、上記推定されたＳｏｏｔ発生質量Ｓｏｏｔnewが目標Ｓｏｏｔ発生質量Ｓｏｏｔterに対して十分に大きい場合、Ｓｏｏｔ発生質量Ｓｏｏｔnewが小さくなるように燃料噴射圧力をフィードバック制御する。
【０１８９】
具体的には、本装置は、Ｓｏｏｔ発生質量Ｓｏｏｔnewから目標Ｓｏｏｔ発生質量Ｓｏｏｔterを減じた値が所定量よりも大きいとき、燃料噴射開始圧力を基本燃料噴射圧力Pcrbaseよりも所定量だけ高くする。これにより、Ｓｏｏｔ発生質量Ｓｏｏｔnewが小さくなる方向に制御される。以上が、燃料噴射制御の概要である。
【０１９０】
（実際の作動）
次に、上記のように構成された内燃機関のすす発生量推定装置の実際の作動について説明する。
【０１９１】
<混合気温度等及びエミッション量の算出>
ＣＰＵ６１は、図６〜図１０に一連のフローチャートにより示した混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンを所定時間の経過毎に、気筒毎に、繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１はステップ６００から処理を開始し、ステップ６０５に進んで吸気弁Vinが開状態から閉状態へと変化したか否か（即ち、ＩＶＣが到来したか否か）を判定し、「Ｎｏ」と判定する場合、ステップ６９５に直ちに進んで本ルーチンを一旦終了する。
【０１９２】
いま、或る気筒においてＩＶＣが到来したものとすると、ＣＰＵ６１はステップ６０５に進んだとき「Ｙｅｓ」と判定してステップ６１０に進み、ＩＶＣ時クランク角度CAivcをクランクポジションセンサ７４から取得される現時点での実クランク角度CAactの値に設定し、ＩＶＣ時筒内ガス圧力Pgivcを吸気管圧力センサ７３から得られる現時点での吸気管圧力Pbの値に設定し、ＩＶＣ時筒内ガス温度Tgivcを吸気温センサ７２から得られる現時点での吸気温度Tbの値に設定し、吸気酸素濃度[O2]inを吸気酸素濃度センサ７７から得られる現時点での吸気酸素濃度RO2inの値に設定する。
【０１９３】
続いて、ＣＰＵ６１はステップ６１５に進んで、上記設定されたＩＶＣ時筒内ガス圧力Pgivcと、上記設定されたＩＶＣ時筒内ガス温度Tgivcと、上記(22)式とに基づいて筒内ガスの全質量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との関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。
【０１９７】
次に、ＣＰＵ６１はステップ６３５に進み、燃料温度センサ７６から得られる現時点での燃料温度Tcrと、上記(35)式とに基づいて燃料蒸気温度Tfを求める。続いて、ＣＰＵ６１はステップ６４０に進んで、上記吸気酸素濃度[O2]inと、[O2]inを引数とする筒内ガス理論空燃比stoichを求めるための関数funcstoichとに基づいて筒内ガス理論空燃比stoichを求める。
【０１９８】
次に、ＣＰＵ６１はステップ６４５に進み、現時点でのエンジン回転速度NEと、微小時間Δｔ（例えば、０．１msec）と、NE，Δｔを引数とする微小クランク角度ΔCAを求めるための関数funcΔCAとに基づいて、同微小時間Δｔに相当するクランク角度である微小クランク角度ΔCAを求める。この微小クランク角度ΔCAは、エンジン回転速度NEが現時点（即ち、ＩＶＣ直後）での値である場合における、微小時間Δｔに相当するクランク角度である。
【０１９９】
続いて、ＣＰＵ６１はステップ６５０に進み、先のステップ６１５にて求めた筒内ガスの全質量Mgを、先のステップ６２５にて求めた燃料噴射開始時期CAinjから得られる燃料噴射開始時筒内容積Vg(CAinj)で除することで燃料噴射開始時筒内ガス密度ρg0を求める。
【０２００】
続いて、ＣＰＵ６１はステップ６５５に進み、先のステップ６１０にて求めたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記燃料噴射開始時筒内容積Vg(CAinj)と、上記(21)式とに基づいて燃料噴射開始時筒内ガス圧力Pg0を求める。
【０２０１】
次に、ＣＰＵ６１はステップ６６０に進んで、先のステップ６３０にて求めた基本燃料噴射圧力Pcrbaseから上記燃料噴射開始時筒内ガス圧力Pg0を減じることで有効噴射圧力ΔPを求め、続くステップ６６５にて、上記求めた有効噴射圧力ΔPと、筒内ガス密度ρg0と、テーブルMapθとに基づいて噴霧角θ（図３を参照）を求める。
【０２０２】
次いで、ＣＰＵ６１は図７のステップ７０５に進み、空気過剰率前回値λbを初期値「０」に設定し、続くステップ７１０にて混合気形成筒内ガス質量Gの値を初期値「０」に設定するとともに、続くステップ７１５にて燃料消費量積算値sumqr、及び筒内ガス消費量積算値sumgrを共に初期値「０」に設定する。
【０２０３】
続いて、ＣＰＵ６１はステップ７２０に進んで、混合気のエンタルピHmixを、上記(34)式に相当する式に従って、初期値（即ち、先のステップ６２０にて求めた指令燃料噴射量Qfinと、燃料の定圧比熱Cfと、先のステップ６３５にて求めた燃料蒸気温度Tfの積）に設定する。
【０２０４】
次に、ＣＰＵ６１はステップ７２５に進み、混合気の定圧比熱Cmixを初期値である上記燃料の定圧比熱Cfに設定し、続くステップ７３０にて混合気質量Mmixを初期値である指令燃料噴射量Qfinに設定する。
【０２０５】
次いで、ＣＰＵ６１はステップ７３５に進んで、混合気内の前駆物質Ｐ濃度[Ｐ]mix、混合気内のＳｏｏｔ濃度[Soot]mix、酸素濃度[O2]mixをそれぞれ初期値「０」に設定するとともに、混合気内の燃料濃度[Fuel]mixを初期値「１」に設定する。続いて、ＣＰＵ６１はステップ７４０に進み、噴射後経過時間ｔを初期値「０」に設定するとともに、クランク角度CAを初期値である燃料噴射開始時期CAinjに設定する。これにより、噴射後経過時間ｔが燃料噴射開始時期CAinjからカウントされることになる。このようにして、各種初期値が決定される。
【０２０６】
次に、ＣＰＵ６１は図８のルーチンに進み、混合気温度Tmixの算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ８０５に進み、噴射後経過時間ｔの値（現時点では、「０」）を微小時間Δｔだけ増大・更新させるとともに、クランク角度CAの値（現時点では、「CAinj」）を微小クランク角度ΔCAだけ増大・更新させる。このように、クランク角度CAの値が噴射後経過時間ｔに対応する値に維持されていく。これにより、以降、噴射後経過時間ｔ＝Δｔとなり、クランク角度CA＝CAinj＋ΔCAとなる。
【０２０７】
続いて、ＣＰＵ６１はステップ８１０に進んで、先のステップ６１５にて求めた筒内ガスの全質量Mgを、上記ステップ８０５にて更新されたクランク角度CAに対応する筒内容積Vg(CA)で除することで噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス密度ρgを求める。
【０２０８】
次に、ＣＰＵ６１はステップ８１５に進んで、上記ステップ６１０にて求められたＩＶＣ時筒内ガス圧力Pgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記クランク角度CAに対応する筒内容積Vg(CA)と、上記(39)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス圧力Pgを求める。
【０２０９】
次いで、ＣＰＵ６１はステップ８２０に進み、上記ステップ６１０にて求められたＩＶＣ時筒内ガス温度Tgivcと、上記ＩＶＣ時筒内容積Vg(CAivc)と、上記クランク角度CAに対応する筒内容積Vg(CA)と、上記(29)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス温度Tgを求める。
【０２１０】
次に、ＣＰＵ６１はステップ８２５に進んで、上記ステップ６１０にて求められた上記吸気酸素濃度[O2]inと、ステップ８２０にて求めた筒内ガス温度Tgと、[O2]in，Tgを引数とする筒内ガスの定圧比熱Cgを求めるための関数funcCgと、上記(30)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガスの定圧比熱Cgを求める。
【０２１１】
続いて、ＣＰＵ６１はステップ８３０に進み、先のステップ８１０にて求めた筒内ガス密度ρgと、先のステップ６６５にて求めた噴霧角θと、先のステップ６６０にて求めた有効噴射圧力ΔPと、先のステップ８０５にて更新した噴射後経過時間ｔと、上記(20)式とに基づいて燃料希釈率dλ/dtを求める。
【０２１２】
次いで、ＣＰＵ６１はステップ８３５に進んで、上記(19)式に従って、空気過剰率λを、その時点での空気過剰率前回値λb（現時点では、ステップ７０５の処理により「０」）に、上記求めた燃料希釈率dλ/dtに微小時間Δｔを乗じた値「dλ/dt・Δｔ」を加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気の空気過剰率λが求められる。
【０２１３】
次に、ＣＰＵ６１はステップ８４０に進み、先のステップ６４０にて求めた筒内ガス理論空燃比stoichと、ステップ８３５にて求めた空気過剰率λと、空気過剰率前回値λb（現時点では、ステップ７０５の処理により「０」。次回からは、後述するステップ８８５にて設定されている値）と、ステップ６２０にて設定された指令燃料噴射量Qfinと、上記(28)式に相当する式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気に新たに取り込まれた筒内ガス質量gを求める。
【０２１４】
続いて、ＣＰＵ６１はステップ８４５に進んで、混合気形成筒内ガス質量Gを、その時点での値（現時点では、ステップ７１０の処理により「０」）に上記新たに取り込まれた筒内ガス質量gを加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気形成筒内ガス質量Gが求められる。
【０２１５】
次に、ＣＰＵ６１はステップ８５０に進み、混合気質量Mmixを、その時点での値（現時点では、ステップ７３０の処理により指令燃料噴射量Qfin）に上記新たに取り込まれた筒内ガス質量gを加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気質量Mmixが求められる。
【０２１６】
続いて、ＣＰＵ６１はステップ８５５に進んで、化学反応前の混合気のエンタルピHpreを、その時点での混合気のエンタルピHmix（現時点では、ステップ７２０の処理により値「Qfin・Cf・Tf」）に、上記(27)式に従って求められる「上記新たに取り込まれた筒内ガスの熱エネルギーΔHg＝g・Cg・Tg」を加えた値に設定する。
【０２１７】
次に、ＣＰＵ６１はステップ８６０に進み、上記求めた化学反応前の混合気のエンタルピHpreを、上記ステップ８５０にて求めた混合気質量Mmixにその時点での混合気の定圧比熱Cmix（現時点では、ステップ７２５の処理により燃料の定圧比熱Cf。次回からは、後述するステップ９２５にて設定されている値）を乗じた値で除することにより上記(33)式に相当する式に従って化学反応前の混合気の温度Tpreを求める。
【０２１８】
次いで、ＣＰＵ６１はステップ８６５に進んで、その時点での混合気内の酸素濃度[O2]mix（現時点では、ステップ７３５の処理により「０」。次回からは、後述するステップ９２０にて設定されている値）と、燃料濃度[Fuel]mix（現時点では、ステップ７３５の処理により「１」。次回からは、後述するステップ９１５にて設定されている値）と、上記求めた化学反応前の混合気の温度Tpreと、上記(32)式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrを求める。
【０２１９】
続いて、ＣＰＵ６１はステップ８７０に進み、上記求めた消費燃料量qrと、上記(31)式とに基づいて微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による反応熱Hrを求め、続くステップ８７５にて混合気のエンタルピHmixを、上記求めた化学反応前の混合気のエンタルピHpreに上記求めた反応熱Hrを加えた値に設定・更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気のエンタルピHmixが求められる。
【０２２０】
そして、ＣＰＵ６１はステップ８８０に進んで、上記ステップ８７５にて求めた混合気のエンタルピHmixと、上記ステップ８５０にて求めた混合気質量Mmixと、その時点での混合気の定圧比熱Cmix（現時点では、ステップ７２５の処理により燃料の定圧比熱Cf。次回からは、後述するステップ９２５にて設定されている値）と、上記(24)式とに基づいて混合気温度Tmixを算出する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj＋ΔCA）における混合気温度Tmixが求められる。
【０２２１】
次に、ＣＰＵ６１はステップ８８５に進んで、次回の計算の準備のため、空気過剰率前回値λbを上記ステップ８３５にて求めた空気過剰率λの値に設定する。この値は、以降、上述したステップ８３５にて使用される。このようにして、噴射後経過時間ｔにおける混合気温度Tmixが算出される。
【０２２２】
次に、ＣＰＵ６１は図９のルーチンに進み、各種濃度の算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ９０５に進み、上記ステップ８６５にて求めた微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrと、ステップ６４０にて求めた筒内ガス理論空燃比stoichと、上記(37)式とに基づいて、微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による筒内ガス消費量grを求める。
【０２２３】
続いて、ＣＰＵ６１はステップ９１０に進んで、燃料消費量積算値sumqrを、その時点での値（現時点では、ステップ７１５の処理により「０」）にステップ８６５にて求めた上記燃料消費量qrを加えた値に設定・更新し、筒内ガス消費量積算値sumgrを、その時点での値（現時点では、ステップ７１５の処理により「０」）にステップ９０５にて求めた上記筒内ガス消費量grを加えた値に設定・更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における燃料消費量積算値sumqr、及び筒内ガス消費量積算値sumgrが求められる。
【０２２４】
次いで、ＣＰＵ６１はステップ９１５に進み、ステップ６２０にて求めた指令燃料噴射量Qfinと、上記求めた燃料消費量積算値sumqrと、ステップ８５０にて求めた混合気質量Mmixと、上記(36)式に相当する式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における混合気内の燃料濃度[Fuel]mixを求める。
【０２２５】
次に、ＣＰＵ６１はステップ９２０に進んで、ステップ８４５にて求めた混合気形成筒内ガス質量Gと、上記求めた筒内ガス消費量積算値sumgrと、ステップ６１０にて設定された吸気酸素濃度[O2]inと、ステップ８５０にて求めた混合気質量Mmixと、上記(38)式とに基づいて、噴射後経過時間ｔ（従って、クランク角度CA）における混合気内の酸素濃度[O2]mixを求める。
【０２２６】
続いて、ＣＰＵ６１はステップ９２５に進み、上記ステップ９２０にて求めた混合気内の酸素濃度[O2]mixと、ステップ８８０にて求めた混合気温度Tmixと、上記(26)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における混合気の定圧比熱Cmixを求める。この値は、以降、ステップ８６０、８８０にて使用される。
【０２２７】
次いで、ＣＰＵ６１はステップ９３０に進み、上記ステップ９１５にて求めた混合気内の燃料濃度[Fuel]mixと、上記ステップ８１５にて求めた筒内ガス圧力Pgと、上記ステップ８８０にて求めた混合気温度Tmixと、上記(3)式とに基づいて前駆物質Ｐの生成速度dmpf/dtを求める。ここで、ステップ９３０は前駆物質生成速度算出手段に対応する。
【０２２８】
続いて、ＣＰＵ６１はステップ９３５に進んで、後述するステップ９６５にて更新されている混合気内の前駆物質Ｐ濃度[Ｐ]mixと、上記筒内ガス圧力Pgと、上記混合気温度Tmixと、上記(4)式とに基づいて前駆物質Ｐの分解速度dmpd/dtを求める。ここで、ステップ９３５は前駆物質分解速度算出手段に対応する。
【０２２９】
次に、ＣＰＵ６１はステップ９４０に進み、上記混合気内の前駆物質Ｐ濃度[Ｐ]mixと、上記筒内ガス圧力Pgと、上記混合気温度Tmixと、上記(16)式とに基づいて上記第１Ｓｏｏｔ生成速度dmsfpp/dtを求め、続くステップ９４５にて、上記混合気内の前駆物質Ｐ濃度[Ｐ]mixと、後述するステップ９７５にて更新されている混合気内のＳｏｏｔ濃度[Soot]mixと、上記筒内ガス圧力Pgと、上記混合気温度Tmixと、上記(13)式とに基づいて上記第２Ｓｏｏｔ生成速度dmsfps/dtを求める。ここで、ステップ９４０は第１すす生成速度算出手段に対応し、ステップ９４５は第２すす生成速度算出手段に対応する。
【０２３０】
次に、ＣＰＵ６１はステップ９５０に進んで、上記第１Ｓｏｏｔ生成速度dmsfpp/dtと、上記第２Ｓｏｏｔ生成速度dmsfps/dtと、上記(5)式とに基づいてＳｏｏｔの生成速度dmsf/dtを求める。ここで、ステップ９５０はすす生成速度算出手段に対応する。
【０２３１】
続いて、ＣＰＵ６１はステップ９５５に進み、上記混合気内のＳｏｏｔ濃度[Soot]mixと、上記混合気内の酸素濃度［０２］mixと、上記筒内ガス圧力Pgと、上記混合気温度Tmixと、上記(17)式とに基づいてＳｏｏｔの酸化速度dmso/dtを求める。ここで、ステップ９５５はすす酸化速度算出手段に対応する。
【０２３２】
次いで、ＣＰＵ６１はステップ９６０に進んで、上記ステップ９３０にて求めた前駆物質Ｐの生成速度dmpf/dtと、上記ステップ９３５にて求めた前駆物質Ｐの分解速度dmpd/dtと、上記ステップ９５０にて求めたＳｏｏｔの生成速度dmsf/dtと、上記(1)式とに基づいて前駆物質Ｐの発生速度d[Ｐ]mix/dtを求める。ここで、ステップ９６０は前駆物質発生速度算出手段に対応する。
【０２３３】
続いて、ＣＰＵ６１はステップステップ９６５に進み、混合気内の前駆物質Ｐ濃度[Ｐ]mixを、その時点での値（現時点では、ステップ７３５の処理により「０」）に、上記求めた前駆物質Ｐの発生速度d「Ｐ」mix/dtに微小時間Δｔを乗じた値「d[Ｐ]mix/dt・Δｔ」を加えた値に設定・更新する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj＋ΔCA）における混合気内の前駆物質Ｐ濃度[Ｐ]mixが求められる。
【０２３４】
次に、ＣＰＵ６１はステップ９７０に進み、上記求めたＳｏｏｔの生成速度dmsf/dtと、上記求めたＳｏｏｔの酸化速度dmso/dtと、上記(2)式とに基づいてＳｏｏｔの発生速度d[Soot]mix/dtを求め、続くステップ９７５にて、混合気内のＳｏｏｔ濃度[Soot]mixを、その時点での値（現時点では、ステップ７３５の処理により「０」）に、上記求めたＳｏｏｔの発生速度d[Soot]mix/dtに微小時間Δｔを乗じた値「d[Soot]mix/dt・Δｔ」を加えた値に設定・更新する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj＋ΔCA）における混合気内のＳｏｏｔ濃度[Soot]mixが求められる。ここで、ステップ９７０はすす発生速度算出手段に対応し、ステップ９７５はすす発生量推定手段に対応する。
【０２３５】
続いて、ＣＰＵ６１は図１０のステップ１００５に進んで、クランク角度CAが圧縮上死点（以下、「ＴＤＣ」と称呼する。）以降であってステップ８８０で求められた噴射後経過時間ｔ（従って、クランク角度CA）における混合気温度TmixがＳｏｏｔの生成反応限界温度Tminより低いか、或いは、クランク角度CAがＴＤＣ以降の所定の終了判定クランク角度CAendに一致したか否かを判定する。
【０２３６】
現時点では、上述のごとく、クランク角度CAは、燃料噴射開始時期クランク角度CAinjに微小クランク角度ΔCAを加えた値（従って、ＴＤＣ前）であるから上記終了判定クランク角度CAendに達していない。
【０２３７】
従って、現時点では、ＣＰＵ６１はステップ１００５にて「Ｎｏ」と判定して図８のステップ８０５に戻り、噴射後経過時間ｔ（現時点では、「1・Δｔ」）を微小時間Δｔだけ増大・更新させるとともに、クランク角度CA（現時点では、「CAinj＋ΔCA」）を微小クランク角度ΔCAだけ増大・更新させた後、上述した図８のステップ８１０〜図１０のステップ１００５の処理を再び実行する。
【０２３８】
これにより、噴射後経過時間ｔ＝２・Δｔ（クランク角度CA＝CAinj＋２・ΔCA）における、混合気の空気過剰率λ（ステップ８３５を参照。）、混合気温度Tmix（ステップ８８０を参照。）、及び各種濃度（[Fuel]mix，[O2]mix，[Ｐ]mix，[Soot]mix。ステップ９１５、９２０、９６５、９７５を参照。）が算出される。
【０２３９】
そして、図１０のステップ１００５の判定において「Ｎｏ」と判定される毎に、図８のステップ８０５〜図１０のステップ１００５の処理が繰り返し実行されていく。これにより、図１０のステップ１００５の判定において「Ｎｏ」と判定される限りにおいて、混合気の空気過剰率λ、混合気温度Tmix、及び各種濃度（[Fuel]mix，[O2]mix，[Ｐ]mix，[Soot]mix）が微小時間Δｔ毎に（即ち、CAinjから微小クランク角度ΔCA毎に）更新されていく。
【０２４０】
そして、ＴＤＣ後の膨張行程における筒内容積の増大等に伴って混合気温度Tmixが低下して上記Ｓｏｏｔの生成反応限界温度Tminより低くなるか、或いは、クランク角度CAが上記終了判定クランク角度CAendに達した場合、ＣＰＵ６１は図１０のステップ１００５に進んだとき「Ｙｅｓ」と判定してステップ１０１０以降に進む。
【０２４１】
ＣＰＵ６１はステップ１０１０に進むと、先のステップ９７５の処理にて更新されている現時点での混合気内のＳｏｏｔ濃度[Soot]mixの値に、先のステップ８５０にて更新されている現時点での混合気質量Mmixを乗じることでＳｏｏｔ発生質量Ｓｏｏｔnewを求める。
【０２４２】
次いで、ＣＰＵ６１はステップ１０１５に進み、Ｓｏｏｔ総質量Ｓｏｏｔを、上記求めたＳｏｏｔ発生質量Ｓｏｏｔnewに所定の残留値Ｓｏｏｔ0を加えた値に設定する。ここで、残留値Ｓｏｏｔ0は、燃料噴射開始時において既に筒内ガス中に含まれている（具体的には、ＥＧＲガス中に含まれている）Ｓｏｏｔ質量に相当する。
【０２４３】
次に、ＣＰＵ６１はステップ１０２０に進み、現時点（即ち、ＩＶＣの直後の時点）におけるエンジン回転速度NEと、上記ステップ６２０にて決定されている指令燃料噴射量Qfinと、テーブルMapＳｏｏｔterとに基づいて目標Ｓｏｏｔ発生質量Ｓｏｏｔterを求める。
【０２４４】
次いで、ＣＰＵ６１はステップ１０２５に進んで、Ｓｏｏｔ発生量偏差ΔＳｏｏｔを、ステップ１０１５にて求められているＳｏｏｔ総質量Ｓｏｏｔから上記求めた目標Ｓｏｏｔ発生質量Ｓｏｏｔterを減じた値に設定する。
【０２４５】
続いて、ＣＰＵ６１はステップ１０３０に進み、上記求めたＳｏｏｔ発生量偏差ΔＳｏｏｔが基準値Ｓｏｏｔrefより大きいか否かを判定し、「Ｙｅｓ」と判定する場合、ステップ１０３５に進んで最終燃料噴射圧力Pcrfinを、上記ステップ６３０にて決定されている基本燃料噴射圧力Pcrbaseに所定値ΔPcrを加えた値に設定する。これにより、Ｓｏｏｔ発生質量Ｓｏｏｔnewが小さくなる方向に燃料噴射圧力が補正されることになる。
【０２４６】
一方、ステップ１０３０にて「Ｎｏ」と判定する場合、ＣＰＵ６１はステップ１０４０に進み、最終燃料噴射圧力Pcrfinを上記基本燃料噴射圧力Pcrbaseと等しい値に設定する。即ち、この場合、燃料噴射圧力が補正されない。
【０２４７】
そして、ＣＰＵ６１はステップ１０４５に進むと、燃料噴射圧力が上記設定された最終燃料噴射圧力Pcrfinになるように燃料噴射用ポンプ２２（の駆動回路）に対して制御指示を行い、ステップ６９５に進んで図６〜図１０の一連の本ルーチンを一旦終了する。以降、ＣＰＵ６１は、次のＩＶＣが到来するまでの間、ステップ６０５に進む毎に「Ｎｏ」と判定し続ける。
【０２４８】
この結果、本ルーチンの実行により、ＩＶＣが到来する毎に、燃料噴射形態（噴射量、噴射圧力、噴射時期）が決定されるとともにＳｏｏｔ発生質量Ｓｏｏｔnew（従って、Ｓｏｏｔ総質量Ｓｏｏｔ）が直ちに推定され、係る推定結果に基づいて噴射圧力が補正されていく。
【０２４９】
また、ＣＰＵ６１は、図１４にフローチャートにより示した燃料噴射制御を行うためのルーチンを所定時間の経過毎に、気筒毎に、繰り返し実行するようになっている。従って、所定のタイミングになると、ＣＰＵ６１はステップ１４００から処理を開始し、ステップ１４０５に進んで実クランク角度CAactが先のステップ６２５にて決定されている燃料噴射開始時期CAinjに一致したか否かを判定し、「Ｎｏ」と判定する場合、ステップ１４９５に直ちに進んで本ルーチンを一旦終了する。
【０２５０】
いま、実クランク角度CAactが上記燃料噴射開始時期CAinjに一致したものとすると、ＣＰＵ６１はステップ１４１０に進んで、対応する燃料噴射弁２１に対してステップ６２０にて決定されている指令燃料噴射量Qfinの燃料の噴射指示（具体的には、燃料噴射期間TAUに亘る開弁指示）を行い、ステップ１４９５に進んで本ルーチンを一旦終了する。これにより、指令燃料噴射量Qfinの燃料が先のステップ１０３５、１０４０の何れかにて設定されている最終燃料噴射圧力Pcrfinをもって噴射される。
【０２５１】
以上、説明したように、本発明によるすす発生量推定装置の実施形態によれば、燃料からの前駆物質Ｐの生成反応、生成された前駆物質Ｐの熱分解反応、生成された前駆物質ＰからのＳｏｏｔの生成反応を考慮して前駆物質Ｐの発生速度d[Ｐ]mix/dt（従って、混合気内の前駆物質Ｐ濃度[Ｐ]mix）を求め、この混合気内の前駆物質Ｐ濃度[Ｐ]mixに依存する前駆物質ＰからのＳｏｏｔの生成反応と、生成されたＳｏｏｔの酸化反応を考慮してＳｏｏｔの発生速度d[Soot]mix/dt（従って、混合気内のＳｏｏｔ濃度[Soot]mix（Ｓｏｏｔ発生量））を推定する。
【０２５２】
このように、燃料からＳｏｏｔが発生する反応過程を、燃料から前駆物質Ｐが発生する反応過程と、前駆物質ＰからＳｏｏｔが発生する反応過程の２段階とする反応モデルを用いることで、燃料からＳｏｏｔが発生する反応過程において実験的・経験的に知られている「Ｓｏｏｔ発生（反応）の遅れ」、並びに、「混合気温度が或る温度を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」を精度良く模擬することができる。この結果、Ｓｏｏｔの発生量が精度良く推定され得る。
【０２５３】
本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記実施形態では、図５に示した反応モデルのように、前駆物質Ｐの発生速度d[Ｐ]mix/dtの算出に際し（図９のステップ９６０を参照）、前駆物質Ｐの生成速度dmpf/dt、及び前駆物質Ｐの分解速度dmpd/dtをそれぞれ、アレニウス型の式である上記(3)式、及び上記(4)式を用いて別個独立に算出しているが（図９のステップ９３０、９３５を参照）、図１５に示した反応モデルのように、前駆物質Ｐの生成速度dmpf/dtから前駆物質Ｐの分解速度dmpd/dtを減じた値である前駆物質Ｐの正味生成速度dmpfd/dtを、ガウス（Gauss）型の式である下記(40)式を用いて一括して算出し、この前駆物質Ｐの正味生成速度dmpfd/dtから、上記(5)式（(16)式＋(13)式）を用いて算出されるＳｏｏｔの生成速度dmsf/dtを減じることで、前駆物質Ｐの発生速度d[Ｐ]mix/dtを算出するように構成してもよい。
【０２５４】
この場合、図９のルーチンにおいて、ステップ９３０、９３５に代えて上記前駆物質Ｐの正味生成速度dmpfd/dtを下記(40)式を用いて算出するステップが設けられ、且つ、ステップ９６０の処理は「d[Ｐ]mix/dt ← dmpfd/dt−dmsf/dt」に変更される。
【０２５５】
【数４０】
0004603951.tif 000041
【０２５６】
この上記(40)式において、Apfd，αpfd，T1（ピーク温度），σ1（標準偏差）は定数である。また、[Fuel]mixは混合気内の燃料濃度、Pgは筒内ガス圧力、Tmixは混合気温度である。上記(40)式によれば、前駆物質Ｐの正味生成速度dmpfd/dtは、混合気温度Tmixが値T1となるときに最大値（ピーク値）を採り、混合気温度Tmixが同値T1から離れるに従って同最大値から減少するように算出される。上記(40)式は、前駆物質正味生成速度算出手段に対応する。
【０２５７】
上記(40)式を用いて前駆物質Ｐの正味生成速度dmpfd/dtを算出すれば、上述した「混合気温度が或る温度（上記値T1）を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」を模擬することができるから、前駆物質Ｐの発生速度d[Ｐ]mix/dtの算出精度を大きく下げることなくＣＰＵ６１の計算負荷を少なくすることができる。
【０２５８】
また、上記実施形態においては、図５に示した反応モデルのように、燃料からＳｏｏｔが発生する反応過程においてＳｏｏｔの前駆物質Ｐを経由させた反応モデルを常に使用しているが、図１６に示した反応モデルのように、燃料からＳｏｏｔが直接発生する反応過程についてのモデルを常に用いてもよい。
【０２５９】
図１６に示した反応モデルでは、Ｓｏｏｔの酸化速度dmso/dtは、上記実施形態と同様、上記(17)式を用いて算出される。一方、燃料からＳｏｏｔが生成される速度（即ち、Ｓｏｏｔの生成速度dmsf/dt）は、ガウス（Gauss）型の式である下記(41)式にて算出される値に所定の遅れ処理（例えば、一次遅れ処理）を施すことで求められる。
【０２６０】
【数４１】
0004603951.tif 000042
【０２６１】
この上記(41)式において、A2，α21，α22，T2（ピーク温度），σ2（標準偏差）は定数である。また、[Fuel]mixは混合気内の燃料濃度、Pgは筒内ガス圧力、Tmixは混合気温度である。上記(41)式によれば、Ｓｏｏｔの生成速度dmsf/dtは、混合気温度Tmixが値T2となるときに最大値（ピーク値）を採り、混合気温度Tmixが同値T2から離れるに従って同最大値から減少するように算出される。
【０２６２】
これにより、上記遅れ処理の実行により上述した「Ｓｏｏｔ発生（反応）の遅れ」を模擬でき、且つ、上述した「混合気温度が或る温度（上記値T2）を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」をも模擬することができる。従って、これによっても、Ｓｏｏｔの発生量が精度良く推定され得る。
【０２６３】
なお、この場合、図９のルーチンにおいて、ステップ９３０〜９６５（ステップ９５５を除く。）が省略され、且つ、ステップ９５０に代えてＳｏｏｔの生成速度dmsf/dtを上記(41)式を用いて算出するステップが設けられる。
【０２６４】
ところで、図１７は、混合気温度Tmix及び混合気の当量比φ（＝1／空気過剰率λ）と、Ｓｏｏｔ発生の有無との関係を示したグラフである。図１７に示した斜線領域がＳｏｏｔが発生する領域に対応している。この図１７については、例えば、日本燃料学会誌  第４６巻  １３６号  （２００４）  「数値シミュレーションを用いたディーゼル機関無煙低音燃焼法のメカニズム解析」  稲垣和久、秋浜一弘等に詳しく紹介されている。
【０２６５】
この図１７から理解できるように、混合気の当量比φが所定値（例えば、「２」）未満では、Ｓｏｏｔが発生し得ない。換言すれば、混合気の当量比φが上記所定値未満では、Ｓｏｏｔの発生量の推定に係る計算を行う必要がない。
【０２６６】
以上のことから、上記実施形態においては、混合気の当量比φにかかわらずＳｏｏｔの発生量の推定に係る計算を行うように構成されているが、混合気の当量比φが所定値（例えば、「２」）以上の場合にのみ、Ｓｏｏｔの発生量の推定に係る計算を行うように構成してもよい。これによると、Ｓｏｏｔ発生量の推定精度を低下することなく計算負荷を少なくすることができる。
【０２６７】
更には、図１７から理解できるように、図１７の斜線領域に対応する「混合気の当量比φと混合気温度Tmixとの関係が所定の範囲内にある場合」にのみ、Ｓｏｏｔ発生量の推定に係る計算を行うように構成してもよい。これによると、混合気の当量比φが上記所定値未満である場合のみならず、混合気の当量比φが同所定値以上であってＳｏｏｔが生成されない場合においてもＳｏｏｔ発生量の推定に係る計算を省略することができ、この結果、Ｓｏｏｔ発生量の推定精度を低下することなく計算負荷をより一層少なくすることができる。
【０２６８】
ところで、上記実施形態においては、図５に示した反応モデルのように、燃料からＳｏｏｔが発生する反応過程においてＳｏｏｔの前駆物質Ｐを経由させた反応モデルを常に使用することで、上述した「Ｓｏｏｔ発生（反応）の遅れ」の模擬を達成している。しかしながら、混合気温度Tmixが或る温度（高温）に達すると、上述した「Ｓｏｏｔ発生（反応）の遅れ」の程度が小さくなる。
【０２６９】
換言すれば、混合気温度Tmixが上記或る温度になると、上記「Ｓｏｏｔの前駆物質を経由させた反応モデル」に代えて図１８に示した反応モデルのように、燃料からＳｏｏｔが直接発生する反応過程についてのモデルを用いても、Ｓｏｏｔの発生量を十分に精度良く推定することができる。
【０２７０】
図１８に示した反応モデルでは、Ｓｏｏｔの酸化速度dmso/dtは、上記実施形態と同様、上記(17)式を用いて算出される。一方、燃料からＳｏｏｔが生成される速度（即ち、Ｓｏｏｔの生成速度dmsf/dt）は、下記(42)式にて算出される。
【０２７１】
【数４２】
0004603951.tif 000043
【０２７２】
この上記(42)式において、A3，α31，α32，zは定数である。また、[Fuel]mixは混合気内の燃料濃度、Pgは筒内ガス圧力、Tmixは混合気温度である。上記(42)式によれば、Ｓｏｏｔの生成速度dmsf/dtは、混合気温度Tmixが（上記或る温度以上において）上昇するにつれて減少していくように算出される。この上記(42)式は、高温時すす生成速度算出手段に対応する。
【０２７３】
これにより、上述した「混合気温度が或る温度を超えた場合において同混合気温度の上昇に応じてＳｏｏｔ発生速度が減少していく現象」を模擬することができる。従って、これによっても、混合気温度Tmixが上記或る温度以上の場合において、Ｓｏｏｔ発生量の推定精度を低下することなく計算負荷を少なくすることができる。
【０２７４】
なお、この場合、図９のルーチンにおいて、ステップ９３０〜９６５（ステップ９５５を除く。）が省略され、且つ、ステップ９５０に代えてＳｏｏｔの生成速度dmsf/dtを上記(42)式を用いて算出するステップが設けられる。
【０２７５】
ところで、上記実施形態においては、燃料の性状を考慮してＳｏｏｔの発生量を推定していない。しかしながら、一般に、例えば、燃料に含まれる芳香族炭化水素の成分の含有率（燃料内におけるアロマ成分の質量濃度。以下、「アロマ成分含有率ｘ」と称呼する。）が大きいほど燃料からＳｏｏｔが生成され易くなることが経験的に知られている。これは、燃料におけるアロマ成分含有率が大きいほど炭素間の結合力が大きくなり、且つ、結合した炭素同士が分離し難くなることに起因すると考えられる。
【０２７６】
従って、Ｓｏｏｔの発生量の推定に際し、アロマ成分含有率ｘを考慮することがＳｏｏｔ発生量を精度良く推定する上で好適である。この場合、例えば、上記実施形態における上記(3)式（図９のステップ９３０）において、燃料濃度[Fuel]mixの値を、「[Fuel]mix・(1＋a・ｘb)」に置き換えることが好適である。ここで、a,bは正の定数である。或いは、上記(3)式において、燃料濃度[Fuel]mixの項を、「([Fuel]mix)c」に置き換えることが好適である。ここで、c＝1＋d・x  （dは正の定数）である。
【０２７７】
これらによれば、アロマ成分含有率ｘの増加に応じて、燃料濃度[Fuel]mixの値を見かけ上大きくして前駆物質Ｐの生成速度dmpf/dt（従って、Ｓｏｏｔの発生量）を大きくすることができる。
【０２７８】
更には、上記実施形態における上記(3)式（図９のステップ９３０）において、活性化エネルギーEpfの値を、「Epf−ｘe・E1」に置き換えてもよい。ここで、e,E1は正の定数である。これによれば、アロマ成分含有率ｘの増加に応じて、活性化エネルギーEpfの値を見かけ上小さくしてＳｏｏｔの生成速度dmsf/dt（従って、Ｓｏｏｔの発生量）を大きくすることができる。
【０２７９】
また、上記実施形態においては、Ｓｏｏｔの酸化速度dmso/dtを上記(17)式に従って算出しているが、所謂ＮＳＣモデルを利用してＳｏｏｔの酸化速度dmso/dtを算出してもよい。なお、ＮＳＣモデルについては、Nagle,
J.Strickland-Constable, R.D., Proc. of the Fifth Carbon Conf., Vol.1 (1962)に詳しく紹介されている。
【０２８０】
また、上記実施形態においては、上記(32)式（図８のステップ８６５）にて算出される燃料消費量qrの対象となる化学反応として、前駆物質Ｐの生成反応、着火反応（熱炎反応）、及び低温酸化反応（冷炎反応）のみならず、その他の種々の化学反応が含まれているが、前駆物質Ｐの生成反応、着火反応、及び低温酸化反応に比してその他の種々の化学反応に係わる燃料消費量は十分に小さいことを考慮して、上記燃料消費量qrの対象となる化学反応として、前駆物質Ｐの生成反応、着火反応、及び低温酸化反応のみを扱うように構成してもよい。
【０２８１】
また、上記実施形態においては、クランク角度CAがＴＤＣ以降であって、且つ、混合気温度TmixがＳｏｏｔの生成反応限界温度Tminより低くなった時点以降、ＣＰＵ６１の計算負荷低減のために前駆物質Ｐの発生速度、及びＳｏｏｔの発生速度の時間積分処理が終了するように構成されているが、これに加え、クランク角度CAがＴＤＣ以前であって、且つ、混合気温度TmixがＳｏｏｔの生成反応限界温度Tminを超えるまでの間も、前駆物質Ｐの発生速度、及びＳｏｏｔの発生速度の時間積分処理を行わないように構成してもよい。これにより、Ｓｏｏｔの発生量の算出に係わる不必要な計算が省略でき、より一層ＣＰＵ６１の計算負荷を低減することができる。
【０２８２】
加えて、上記実施形態においては、種々の反応速度を算出する式（具体的には、上記(3)式、(4)式、(13)式、(16)式、(17)式）の何れかにおいて、筒内ガス圧力Pgに係る項を省略してもよい。
【図面の簡単な説明】
【０２８３】
【図１】本発明の実施形態に係る内燃機関のすす発生量推定装置を４気筒内燃機関（ディーゼル機関）に適用したシステム全体の概略構成図である。
【図２】或る一つの気筒のシリンダ内(筒内)に吸気マニホールドからガスが吸入され、筒内に吸入されたガスが排気マニホールドへ排出される様子を模式的に示した図である。
【図３】（ａ）は、燃料が一時に噴射された時点での同燃料の様子を模式的に示した図であり、（ｂ）は、筒内ガスと混ざり合いながら混合気となって円錐状に拡散していく燃料の様子を模式的に示した図である。
【図４】前駆物質の例としてＰＡＨ（多環芳香族炭化水素類）を挙げた場合における、燃料から前駆物質を経由してＳｏｏｔが発生する反応過程を模式的に示した図である。
【図５】図１に示したすす発生量推定装置が使用する、燃料からＳｏｏｔが発生する反応過程において前駆物質が経由される反応モデルを表した図である。
【図６】図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。
【図７】図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。
【図８】図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。
【図９】図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。
【図１０】図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第５番目の部分を示したフローチャートである。
【図１１】図１に示したＣＰＵが図６に示したルーチンを実行する際に参照する指令燃料噴射量を決定するためのテーブルである。
【図１２】図１に示したＣＰＵが図６に示したルーチンを実行する際に参照する燃料燃料噴射時期を決定するためのテーブルである。
【図１３】図１に示したＣＰＵが図６に示したルーチンを実行する際に参照する基本燃料噴射圧力を決定するためのテーブルである。
【図１４】図１に示したＣＰＵが実行する燃料噴射制御を行うためのルーチンを示したフローチャートである。
【図１５】本発明の実施形態の変形例に係るすす発生量推定装置が使用する、燃料からＳｏｏｔが直接発生する反応過程についての反応モデルを表した図である。
【図１６】本発明の実施形態の他の変形例に係るすす発生量推定装置が使用する、燃料からＳｏｏｔが直接発生する反応過程についての反応モデルを表した図である。
【図１７】混合気温度及び混合気の当量比と、Ｓｏｏｔ発生の有無との関係を示したグラフである。
【図１８】本発明の実施形態の他の変形例に係るすす発生量推定装置が使用する、燃料からＳｏｏｔが直接発生する反応過程についての反応モデルを表した図である。
【符号の説明】
【０２８４】
２１…燃料噴射弁、６０…電気制御装置、６１…ＣＰＵ、７２…吸気温センサ、７３…吸気管圧力センサ、７４…クランクポジションセンサ、７６…燃料温度センサ、７７…吸気酸素濃度センサ
【図１】
0004603951.tif 000044
【図２】
0004603951.tif 000045
【図３】
0004603951.tif 000046
【図４】
0004603951.tif 000047
【図５】
0004603951.tif 000048
【図６】
0004603951.tif 000049
【図７】
0004603951.tif 000050
【図８】
0004603951.tif 000051
【図９】
0004603951.tif 000052
【図１０】
0004603951.tif 000053
【図１１】
0004603951.tif 000054
【図１２】
0004603951.tif 000055
【図１３】
0004603951.tif 000056
【図１４】
0004603951.tif 000057
【図１５】
0004603951.tif 000058
【図１６】
0004603951.tif 000059
【図１７】
0004603951.tif 000060
【図１８】
0004603951.tif 000061