JP2009138688A - 内燃機関のすす発生量推定装置 - Google Patents

Abstract

【課題】比較的低温雰囲気下でのＳｏｏｔの酸化反応を考慮してＳｏｏｔの発生量を精度良く推定できるすす発生量推定装置を提供すること。
【解決手段】この装置は、Ｓｏｏｔの酸化反応として、Ｓｏｏｔを構成する多数の炭素原子のうちで安定状態にあるものと酸化剤（Ｏ_２、Ｏ）との衝突によりＳｏｏｔが酸化される反応（Ｃ＋Ｏ_２）、Ｓｏｏｔを構成する多数の炭素原子のうちでラジカル部位を有するものと酸化剤（Ｏ_２、Ｏ）との衝突によりＳｏｏｔが酸化される反応（Ｃラジカル＋Ｏ_２）、及び、Ｓｏｏｔを構成する多数の炭素原子のうちで安定状態にあるものと水酸化ラジカル（ＯＨ）との衝突によりＳｏｏｔが酸化される反応（Ｃ＋ＯＨ）が考慮されて、生成されたＳｏｏｔの酸化速度が算出され、このＳｏｏｔの酸化速度に基づいてＳｏｏｔの発生量が推定される。
【選択図】図５

Description

本発明は、内燃機関の燃焼室内において燃料の反応に起因して発生するすす（カーボン微粒子。以下、「Ｓｏｏｔ」とも称呼する。）の発生量を推定するすす発生量推定装置に関する。

内燃機関（特に、ディーゼル機関）の燃焼室内で発生する粒子状物質（パティキュレート・マター（ＰＭ））を構成する主たる成分の一つはＳｏｏｔである。このＳｏｏｔの発生量を精度良く制御してＳｏｏｔの発生量を少なくするためには、Ｓｏｏｔの発生量を精度良く推定する必要がある。

例えば、下記特許文献１に記載の内燃機関のすす発生量推定装置では、燃料噴射後において燃料の反応に起因してすすが生成される速度（Ｓｏｏｔ生成速度）が逐次算出される。燃料噴射後において生成されたすすが酸化される速度（Ｓｏｏｔ酸化速度）も逐次算出される。このようにして求められる燃料噴射後において時々刻々と変化していくＳｏｏｔ生成速度及びＳｏｏｔ酸化速度に基づいて燃料噴射後においてすすが発生する速度（Ｓｏｏｔ発生速度）が逐次算出される（Ｓｏｏｔ発生速度＝Ｓｏｏｔ生成速度−Ｓｏｏｔ酸化速度）。そして、燃料噴射後において時々刻々と変化していくＳｏｏｔ発生速度を逐次積算していくことでＳｏｏｔの発生量が推定される。
特開２００７−４６４７７号公報

上記文献に記載の装置では、Ｓｏｏｔ酸化速度は、日本機械学会論文集(Ｂ編) 48巻432号「直接噴射式ディーゼル機関の燃焼モデルと性能予測」（以下、「非特許文献１」と称呼する。）にて紹介された実験式（所謂廣安（酸化）モデル）に従って計算されている。また、Ｓｏｏｔ酸化速度を算出するモデルとして、Nagle J et al., ‘Oxidation of carbon between 1000-2000℃’, Proc.of the Fifth Carbon
Conf., Vol.1(1962) （以下、「非特許文献２」と称呼する。）にて紹介された実験式（所謂ＮＳＣモデル）も広く知られている。

上述の非特許文献１、２に紹介されたモデルでは、Ｓｏｏｔの酸化反応として、比較的高温雰囲気下（例えば、１３００Ｋ以上）において、Ｓｏｏｔを構成する多数の炭素（原子）のうちで安定状態にあるもの（ラジカル部位を有さないもの）と酸化剤（酸素分子Ｏ_２、酸素原子Ｏ等）との衝突により炭素（原子）間の結合（Ｃ=Ｃ結合、共有結合）が崩壊してＳｏｏｔが酸化される反応のみが想定されている。ここで、「ラジカル部位」は、活性基、反応基、不対電子とも称呼される。

以下、この反応を「高温域での酸化剤による炭素の酸化反応」とも称呼することもある。この「高温域での酸化剤による炭素の酸化反応」は、比較的低温雰囲気下（例えば、１３００Ｋ未満）では殆ど発生し得ない。従って、上述のモデルを使用すると、比較的低温雰囲気下においてＳｏｏｔ酸化速度が無視できる程度に十分小さい値に算出される。

ところで、本発明者の実験・研究によれば、比較的低温雰囲気下（例えば、１３００Ｋ未満）でも、実際には、Ｓｏｏｔの酸化反応が無視できない程度で発生し得ることが分かってきている。そして、この反応は、Ｓｏｏｔを構成する多数の炭素（原子）のうちでラジカル部位を有するもの（以下、「炭素ラジカル」とも称呼される。）と酸化剤（酸素分子Ｏ_２、酸素原子Ｏ等）との衝突によりＳｏｏｔが酸化される反応（ラジカル反応）であることも分かってきている。以下、この反応を「低温域での酸化剤による炭素ラジカルの酸化反応」とも称呼することもある。

従って、比較的低温雰囲気下にて支配的となる「低温域での酸化剤による炭素ラジカルの酸化反応」を考慮していない上述のモデルでＳｏｏｔ酸化速度を逐次計算していくと、比較的低温雰囲気下にてＳｏｏｔ酸素速度が真値よりも小さめに計算され続ける。よって、比較的低温雰囲気下にてＳｏｏｔ発生速度（＝Ｓｏｏｔ生成速度−Ｓｏｏｔ酸化速度）が真値よりも大きめに計算され続け、この結果、Ｓｏｏｔ発生速度の積算値であるＳｏｏｔの発生量が真値よりも大きめに推定されるという問題が発生し得る。

本発明は、上述の問題を解決するためになされたものであり、その目的は、比較的低温雰囲気下でのＳｏｏｔの酸化反応を考慮してＳｏｏｔの発生量を精度良く推定できるすす発生量推定装置を提供することにある。

本発明に係るすす発生量推定装置は、内燃機関の燃焼室内において燃料の反応に起因して発生するすすの発生量を推定するすす発生量推定手段を備えている。そして、本発明に係るすす発生量推定装置の特徴は、前記すす発生量推定手段が、前記燃焼室内の温度が所定温度未満の場合において、生成されたすすを構成する炭素であってラジカル部位を有するもの（上記炭素ラジカル）と酸化剤（Ｏ_２，Ｏ）との衝突により前記生成されたすすが酸化される反応を考慮して前記すすの発生量を推定するように構成されたことにある。

より具体的には、例えば、前記すす発生量推定手段は、前記燃料の反応に起因して前記すすが生成される速度であるすす生成速度（＝Ｓｏｏｔ生成速度）を算出するすす生成速度算出手段と、前記生成されたすすが酸化される速度であるすす酸化速度（＝Ｓｏｏｔ酸化速度）を算出するすす酸化速度算出手段と、前記すす生成速度と前記すす酸化速度とに基づいて前記すすが発生する速度であるすす発生速度（＝Ｓｏｏｔ発生速度）を算出するすす発生速度算出手段とを備え、前記すす発生速度に基づいて前記すすの発生量を推定するように構成され得る。

このように前記すす発生量推定手段が構成される場合において、前記すす酸化速度算出手段は、前記ラジカル部位を有する炭素と前記酸化剤との衝突に起因して前記生成されたすすが酸化される速度である第１すす酸化速度を算出する第１すす酸化速度算出手段を備え、前記燃焼室内の温度が前記所定温度未満の場合において、前記第１すす酸化速度に基づいて前記すす酸化速度を算出するように構成され得る。この第１すす酸化速度は、上述の「低温域での酸化剤による炭素ラジカルの酸化反応」に基づく酸化速度に対応する。

ここにおいて、前記「燃焼室内の温度」とは、例えば、燃焼室内において燃料と筒内ガスとが混ざり合って形成される混合気の温度等を指す。また、前記「ラジカル部位」とは活性基、反応基、不対電子に対応している。また、前記「所定温度」とは、例えば、１３００Ｋ（或いは、１３００Ｋを含む所定の範囲、例えば、１５００Ｋ）である。

上記構成によれば、比較的低温雰囲気下において支配的となる「低温域での酸化剤による炭素ラジカルの酸化反応」（前記第１すす酸化速度に対応）が考慮されてＳｏｏｔ酸化速度が（逐次）計算される。従って、比較的低温雰囲気下においてＳｏｏｔ酸化速度が（逐次）精度良く計算され得る。加えて、このように精度良く計算され得るＳｏｏｔ酸化速度に基づいてＳｏｏｔ発生速度が（逐次）精度良く計算され得る。Ｓｏｏｔ発生速度とは、例えば、混合気内のＳｏｏｔ量の変化速度、混合気内のＳｏｏｔ濃度の変化速度等であり、具体的には、Ｓｏｏｔ生成速度からＳｏｏｔ酸化速度を減じることで算出され得る。

そして、このように精度良く計算され得るＳｏｏｔ発生速度に基づいて（具体的には、Ｓｏｏｔ発生速度を時間積分（積算）することで）Ｓｏｏｔの発生量が推定される。Ｓｏｏｔの発生量とは、例えば、混合気内のＳｏｏｔ量、混合気内のＳｏｏｔ濃度等である。このように、Ｓｏｏｔの発生量が比較的低温雰囲気下において支配的となる「低温域での酸化剤による炭素ラジカルの酸化反応」が考慮されて推定されるから、Ｓｏｏｔの発生量が精度良く推定され得る。なお、上述の「Ｓｏｏｔ発生速度の時間積分」が行われる期間は、例えば、燃料噴射時点から、膨張行程において燃料の反応が終了すると判定される程度まで燃焼室内の温度が低下した時点まで、の期間等である。

ここで、一般に、低圧縮比の内燃機関では燃焼温度が低い。従って、燃料噴射から膨張行程において燃料の反応が終了するまでの期間内において比較的低温雰囲気（例えば、１３００Ｋ未満）となる期間が占める割合が大きくなる。係る観点からみると、上記本発明に係るすす発生量推定装置は、例えば、圧縮比が１６以下の低圧縮比の内燃機関（特に、ディーゼル機関）に適用される場合において、Ｓｏｏｔ酸化速度が（逐次）精度良く計算され得る期間が長くなって、Ｓｏｏｔの発生量が特に精度良く推定され得る。

上記本発明に係るすす発生量推定装置においては、前記すす酸化速度算出手段は、前記生成されたすすを構成する炭素であって安定状態にあるものと酸化剤との衝突に起因して前記生成されたすすが酸化される速度である第２すす酸化速度を算出する第２すす酸化速度算出手段を備え、前記燃焼室内の温度が前記所定温度未満の場合では、前記第１すす酸化速度のみに基づいて前記すす酸化速度を算出するとともに、前記燃焼室内の温度が前記所定温度以上の場合では、前記第１すす酸化速度と前記第２すす酸化速度とに基づいて前記すす酸化速度を算出するように構成されることが好適である。この第２すす酸化速度は、上述の「高温域での酸化剤による炭素の酸化反応」に基づく酸化速度に対応する。なお、燃焼室内の温度が所定温度未満の場合（＝比較的低温雰囲気に対応）とは、具体的には、燃料噴射から、燃焼反応により増大していく燃焼室内温度が所定温度に達するまでの期間に対応し、燃焼室内の温度が所定温度以上の場合（＝比較的高温雰囲気に対応）とは、具体的には、燃焼室内温度が上昇しながら所定温度に達した時点以降（且つ、ピーク温度を経て所定温度に復帰するまで）の期間に対応する。

比較的高温雰囲気下では、上述の「低温域での酸化剤による炭素ラジカルの酸化反応」に加えて上述の「高温域での酸化剤による炭素の酸化反応」が複合的に発生する。この「高温域での酸化剤による炭素の酸化反応」は、温度が上昇するにつれてより活発となる。上記構成によれば、比較的高温雰囲気下において、「低温域での酸化剤による炭素ラジカルの酸化反応」（前記第１すす酸化速度に対応）に加えて「高温域での酸化剤による炭素の酸化反応」（前記第２すす酸化速度に対応）が考慮されてＳｏｏｔ酸化速度が計算される。従って、比較的高温雰囲気下においても、Ｓｏｏｔ酸化速度（従って、Ｓｏｏｔの発生量）が精度良く計算され得る。

加えて、上述のように、比較的低温雰囲気下では、「高温域での酸化剤による炭素の酸化反応」が殆ど発生しない（即ち、前記第２すす酸化速度が無視できる程度に十分小さい値に算出される）。換言すれば、「高温域での酸化剤による炭素の酸化反応」を考慮することなく「低温域での酸化剤による炭素ラジカルの酸化反応」のみを考慮してＳｏｏｔ酸化速度を計算してもＳｏｏｔ酸化速度の計算精度が低下しない。上記構成によれば、比較的低温雰囲気下において、「低温域での酸化剤による炭素ラジカルの酸化反応」（前記第１すす酸化速度に対応）のみが考慮されてＳｏｏｔ酸化速度が計算される。従って、比較的低温雰囲気下において、Ｓｏｏｔ酸化速度（従って、Ｓｏｏｔの発生量）の計算精度を下げることなく、前記第２すす酸化速度の計算そのものを省略してＳｏｏｔ酸化速度の計算に要する負荷を軽減できる。

ところで、燃焼室内では、燃料と酸素分子との反応等に起因して、水酸化ラジカル（ＯＨ）が発生し得る。本発明者の実験・研究によれば、比較的高温雰囲気下（例えば、１３００Ｋ以上）にて、Ｓｏｏｔを構成する多数の炭素（原子）のうちで安定状態にあるもの（ラジカル部位を有さないもの）と水酸化ラジカル（ＯＨ）との衝突によりＳｏｏｔが酸化される反応（ラジカル反応）が発生することも分かってきている。以下、この反応を「高温域での水酸化ラジカルによる炭素の酸化反応」と称呼することもある。

即ち、比較的高温雰囲気下では、上述の「低温域での酸化剤による炭素ラジカルの酸化反応」に加えて「高温域での水酸化ラジカルによる炭素の酸化反応」も複合的に発生する。この「高温域での水酸化ラジカルによる炭素の酸化反応」も、温度が上昇するにつれてより活発となる。

係る観点に基づき、上記本発明に係るすす発生量推定装置においては、前記すす酸化速度算出手段は、前記生成されたすすを構成する炭素であって安定状態にあるものと燃料の反応に起因して発生する水酸化ラジカルとの衝突に起因して前記生成されたすすが酸化される速度である第３すす酸化速度を算出する第３すす酸化速度算出手段を備え、前記燃焼室内の温度が前記所定温度未満の場合では、前記第１すす酸化速度のみに基づいて前記すす酸化速度を算出するとともに、前記燃焼室内の温度が前記所定温度以上の場合では、前記第１すす酸化速度と前記第３すす酸化速度とに基づいて前記すす酸化速度を算出するように構成されることが好適である。

上記構成によれば、比較的高温雰囲気下において、「低温域での酸化剤による炭素ラジカルの酸化反応」（前記第１すす酸化速度に対応）に加えて「高温域での水酸化ラジカルによる炭素の酸化反応」（前記第３すす酸化速度に対応）が考慮されてＳｏｏｔ酸化速度が計算される。従って、上記構成によっても、比較的高温雰囲気下において、Ｓｏｏｔ酸化速度（従って、Ｓｏｏｔの発生量）が精度良く計算され得る。

加えて、比較的低温雰囲気下では、上述の「高温域での水酸化ラジカルによる炭素の酸化反応」が殆ど発生しない（即ち、前記第３すす酸化速度が無視できる程度に十分小さい値に算出される）。換言すれば、「高温域での水酸化ラジカルによる炭素の酸化反応」を考慮することなく「低温域での酸化剤による炭素ラジカルの酸化反応」のみを考慮してＳｏｏｔ酸化速度を計算してもＳｏｏｔ酸化速度の計算精度が低下しない。従って、上記構成によっても、比較的低温雰囲気下において、Ｓｏｏｔ酸化速度（従って、Ｓｏｏｔの発生量）の計算精度を下げることなく、前記第３すす酸化速度の計算そのものを省略してＳｏｏｔ酸化速度の計算に要する負荷を軽減できる。

更には、前記すす酸化速度算出手段は、前記燃焼室内の温度が前記所定温度未満の場合では、前記第１すす酸化速度のみに基づいて前記すす酸化速度を算出するとともに、前記燃焼室内の温度が前記所定温度以上の場合では、前記第１すす酸化速度と前記第２すす酸化速度と前記第３すす酸化速度とに基づいて前記すす酸化速度を算出するように構成されることがより好ましい。

これによれば、比較的低温雰囲気下において前記第２、第３すす酸化速度の計算そのものを省略してＳｏｏｔ酸化速度の計算に要する負荷を大幅に軽減できる。加えて、比較的高温雰囲気下において「低温域での酸化剤による炭素ラジカルの酸化反応」に加えて「高温域での酸化剤による炭素の酸化反応」と「高温域での水酸化ラジカルによる炭素の酸化反応」とが共に考慮されてＳｏｏｔ酸化速度が計算されるから、比較的高温雰囲気下において、Ｓｏｏｔ酸化速度（従って、Ｓｏｏｔの発生量）がより一層精度良く計算され得る。

また、上記本発明に係るすす発生量推定装置においては、前記第１すす酸化速度算出手段は、燃料の反応に起因して発生する水酸化ラジカルの濃度を用いることなく記述される関係に基づいて前記第１すす酸化速度を算出するように構成されることが好適である。この場合、例えば、前記第１すす酸化速度算出手段は、前記燃焼室内のすすの濃度、前記燃焼室内の酸素の濃度、前記燃焼室内の温度のべき乗、及び前記燃焼室内の温度（並びに、活性化エネルギー）に基づいて決定される値を指数とする指数関数、の積により表される関係に基づいて前記第１すす酸化速度を算出することができる。

厳密には、「低温域での酸化剤による炭素ラジカルの酸化反応」の反応過程において、水酸化ラジカルが関係する。しかしながら、一般に、この反応過程に関して、水酸化ラジカルの濃度（混合気内の濃度）を用いて記述される関係（理論式、関数等）は複雑となり、この関係に基づく計算に要する負荷は大きくなる傾向がある。上記構成によれば、水酸化ラジカルの濃度を用いることなく記述される関係（簡略式、実験式等）に基づいて前記第１すす酸化速度を算出できるから、前記第１すす酸化速度の計算に要する負荷を大幅に軽減することができる。

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

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

エンジン本体２０の各気筒の上部には燃料噴射弁（噴射弁、インジェクタ）２１が配設されている。各燃料噴射弁２１は、図示しない燃料タンクと接続された燃料噴射用ポンプ２２に燃料配管２３を介して接続されている。燃料噴射用ポンプ２２は、電気制御装置６０と電気的に接続されていて、同電気制御装置６０からの駆動信号（後述する指令最終燃料噴射圧力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)式に従って表すことができる。

上記(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)式に従って表すことができる。

上記(2)式において、dmsf/dtは、上述したＳｏｏｔの生成速度であり、dmso/dtは、上記生成されたＳｏｏｔの酸化反応に基づく混合気内のＳｏｏｔ濃度[Soot]mixの減少速度（以下、「Ｓｏｏｔの酸化速度」と称呼する。）である。

（前駆物質Ｐ濃度の取得）
以下、先ず、上記(1)式に基づいて前駆物質Ｐ濃度を取得する方法について説明する。前駆物質Ｐ濃度を取得するためには、上記(1)式の右辺の各項の値を求める必要がある。

<前駆物質Ｐの生成速度の取得>
上記(1)式の右辺第１項の前駆物質Ｐの生成速度dmpf/dtは、気相反応についての反応速度を算出する場合に使用され得る所謂アレニウス（Arrhenius）の式に基づく下記(3)式に従って求めることができる。燃料からの前駆物質Ｐの生成反応は、気相反応とみなせるからである。

上記(3)式において、Apf，αpfは定数、Epfは活性化エネルギー（定数）、Ｒは気体定数（定数。以下も同じ。）である。また、[Fuel]mixは混合気内の燃料濃度（以下も同じ）、Pgは筒内ガスの圧力（以下も同じ）、Tmixは混合気の温度（以下も同じ）であり、これらは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、前駆物質Ｐの生成速度dmpf/dtも、上記(3)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

<前駆物質Ｐの分解速度の取得>
上記(1)式の右辺第２項の前駆物質Ｐの分解速度dmpd/dtも、上記(3)式と同様、アレニウスの式に基づく下記(4)式に従って求めることができる。上記生成された前駆物質Ｐの熱分解反応も、気相反応とみなせるからである。

上記(4)式において、Apd，αpdは定数、Epdは活性化エネルギ（定数）である。また、[Ｐ]mixは混合気内の前駆物質Ｐ濃度であり、上記(1)式に従って噴射後経過時間ｔ＝０から微小時間Δｔ毎に得られる前駆物質Ｐの発生速度d[Ｐ]mix/dtを時間で積分（積算）していくことで、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、前駆物質Ｐの分解速度dmpd/dtも、上記(4)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

<Ｓｏｏｔの生成速度の取得>
上記(1)式の右辺第３項のＳｏｏｔの生成速度dmsf/dtは、下記(5)式に従って表すことができる。下記(5)式において、dmsfpp/dtは前駆物質Ｐ同士の衝突に起因するＳｏｏｔの生成反応に基づくＳｏｏｔの生成速度（以下、「第１Ｓｏｏｔ生成速度」と称呼する。）であり、dmsfps/dtは前駆物質ＰとＳｏｏｔの衝突に起因するＳｏｏｔの生成反応に基づくＳｏｏｔの生成速度（以下、「第２Ｓｏｏｔ生成速度」と称呼する。）である。

このように、この反応モデルでは、前駆物質からＳｏｏｔが生成される反応過程において、前駆物質Ｐ同士の衝突によるＳｏｏｔ生成反応と、Ｓｏｏｔが生成された後における前駆物質ＰとＳｏｏｔとの衝突によるＳｏｏｔ生成反応とが別個独立に考慮されてＳｏｏｔの生成速度dmsf/dtが算出される。

前駆物質ＰからＳｏｏｔが生成される反応は、分子量の比較的大きい粒子に関わる反応である。従って、前駆物質ＰからＳｏｏｔが生成される反応を気相反応とみなすことは困難である。この反応モデルでは、上記(5)式の右辺第１項の第１Ｓｏｏｔ生成速度dmsfpp/dt、及び上記(5)式の右辺第２項の第２Ｓｏｏｔ生成速度dmsfps/dtは、以下のように求める。

一般に、粒子ａと粒子ｂの衝突による凝集反応の速度は、ガス中（例えば、混合気中）における粒子ａと粒子ｂの間の「衝突数」（ガス単位質量あたりの衝突数）と、衝突した粒子ａと粒子ｂが凝集する確率（以下、「凝集確率」と称呼する。）の積で表すことができる。

また、一般に、ガス中における粒子ａと粒子ｂの間の「衝突数」は、粒子ａと粒子ｂの間の「相対速度」と、粒子ａと粒子ｂの衝突に関わる面積である「衝突断面積」と、粒子ａの「粒子数密度」（ガス単位質量あたりの粒子ａの粒子数）と、粒子ｂの「粒子数密度」（ガス単位質量あたりの粒子ｂの粒子数）の積で表すことができる。ここで、「衝突断面積」は、粒子ａの半径と粒子ｂの半径の和を半径とする円の面積であり、下記(6)式で表すことができる。ここで、Da，Dbはそれぞれ、粒子ａ、粒子ｂの直径（以下も同じ）である。

また、「粒子数密度」は、ガス中における粒子の濃度を１粒子の質量で除した値である。従って、粒子ａの「粒子数密度」は「Ｍa/ｍa」で、粒子ｂの「粒子数密度」は「Ｍb/ｍb」で表すことができる。ここで、Ｍa，Ｍbはそれぞれ、粒子ａ、粒子ｂのガス中における濃度であり、ｍa，ｍbはそれぞれ、粒子ａ、粒子ｂの１粒子の質量（以下も同じ）である。また、ρa，ρbをそれぞれ、粒子ａ、粒子ｂの１粒子の密度とすると、粒子ａの「粒子数密度」は「(6・Ｍa)/(ρa・π・Da³)」と、粒子ｂの「粒子数密度」は「(6・Ｍb)/(ρb・π・Db³)」と表すこともできる。

また、粒子ａと粒子ｂの間の「相対速度」は、マクスウエル−ボルツマン（Maxwell−Boltzmann）の分布関数を利用すれば、下記(7)式で表すことができる。ここで、μは粒子ａと粒子ｂの換算質量であって、μ＝(ｍa・ｍb)/(ｍa＋ｍb)で表される。ｋはボルツマン定数であり、Ｔはガスの温度である。

従って、ガス中における粒子ａと粒子ｂの間の「衝突数」は、下記(8)式で表すことができる。

以上のことから、衝突した粒子ａと粒子ｂが凝集する上記凝集確率をＡ１とし、下記(9)式に示すように係数Ａ２を定義すると、粒子ａと粒子ｂの衝突による凝集反応の速度は、一般に、下記(10)式に従って表すことができる。

ここで、上記第１Ｓｏｏｔ生成速度dmsfpp/dtは、前駆物質Ｐ同士の衝突による凝集反応の速度であり、上記第２Ｓｏｏｔ生成速度dmsfps/dtは、前駆物質ＰとＳｏｏｔとの衝突による凝集反応の速度である。従って、この反応モデルでは、上記第１Ｓｏｏｔ生成速度dmsfpp/dt、及び第２Ｓｏｏｔ生成速度dmsfps/dtは、上記(10)式を用いて求められる。

先ず、上記第２Ｓｏｏｔ生成速度dmsfps/dtについて説明する。この場合、上記粒子ａが前駆物質Ｐの粒子に、上記粒子ｂがＳｏｏｔの粒子に対応すると考える。そうすると、上記第２Ｓｏｏｔ生成速度dmsfps/dtは、上記(10)式に対応する下記(11)式に従って求めることができる。ここで、Asfps1は、衝突した前駆物質ＰとＳｏｏｔとが凝集する凝集確率であり、Asfps2は、下記(12)式に示すように定義される、上記係数Ａ２に対応する係数である。

上記(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)式に従って求められる。

上記(13)式において、上述したように、混合気内の前駆物質Ｐ濃度[Ｐ]mix、筒内ガス圧力Pg、混合気温度Tmixは、上述したように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。また、混合気内のＳｏｏｔ濃度[Soot]mixも、後述するように噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上より、第２Ｓｏｏｔ生成速度dmsfps/dtも、上記(13)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

このように、この反応モデルでは、第２Ｓｏｏｔ生成速度dmsfps/dtは、前駆物質Ｐの粒子とＳｏｏｔの粒子の衝突数、及び衝突した前駆物質Ｐの粒子とＳｏｏｔの粒子とが凝集する上記凝集穫率Asfps1を考慮して算出される。換言すれば、第２Ｓｏｏｔ生成速度dmsfps/dtは、前駆物質Ｐの粒子とＳｏｏｔの粒子との衝突に起因して新たなＳｏｏｔが生成される確率を考慮して算出される。加えて、上記「前駆物質Ｐの粒子とＳｏｏｔの粒子の衝突数」は、混合気中の前駆物質Ｐ濃度[Ｐ]mixと、混合気中のＳｏｏｔ濃度[Soot]mixと、混合気温度Tmixとに基づいて算出される。

次に、上記第１Ｓｏｏｔ生成速度dmsfpp/dtについて説明する。この場合、上記粒子ａ、上記粒子ｂ共に、前駆物質Ｐの粒子に対応すると考える。そうすると、換算質量μ＝mp/2＝(ρp・π・Dp³)/6・(1/2)と表すことができる。また、この場合、同じ種類の粒子の衝突が扱われていることに起因して、一回の衝突現象を２回数えることになる。即ち、前駆物質Ｐ同士の「衝突数」は、上記(8)式から得られる値に(1/2)を乗じた値になる。

以上のことを考慮すると、上記第１Ｓｏｏｔ生成速度dmsfpp/dtは、上記(10)式に対応する下記(14)式に従って求めることができる。ここで、Asfpp1は、衝突した前駆物質Ｐ同士が凝集する凝集確率であり、Asfpp2は、下記(15)式に示すように定義される、上記係数Ａ２に対応する係数である。

ここで、この反応モデルでは、上述したように、Dp，ρpは一定値として扱う。従って、上記(15)式から理解できるように、上記係数Asfpp2は一定値となる。また、上記凝集確率Asfpp1も一定であると仮定する。また、このモデルでは、前駆物質Ｐ以外の他の粒子の衝突による上記第１Ｓｏｏｔ生成速度への影響も考慮するため、上記(14)式の右辺の因数として筒内ガス圧力Pgに依存する項「Pg^αsfpp」が加えられる（αsfppは定数）。即ち、このモデルでは、上記第１Ｓｏｏｔ生成速度dmsfpp/dtは、上記係数Asfpp2と上記凝集確率Asfpp1の積を係数Asfpp（定数）と置くと、下記(16)式に従って求められる。

上記(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)式の右辺の各項の値を求める必要がある。

<Ｓｏｏｔの生成速度の取得>
上記(2)式の右辺第１項のＳｏｏｔ生成速度dmsf/dtは、上述した上記(1)式の右辺第３項のＳｏｏｔ生成速度dmsf/dtと同じであり、上記(5)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

<Ｓｏｏｔの酸化速度の取得>
この反応モデルでは、上記(2)式の右辺第２項のＳｏｏｔ酸化速度dmso/dtの算出にあたり、Ｓｏｏｔの酸化反応として、以下の３つの酸化反応が考慮される。

<<高温域での酸化剤による炭素の酸化反応>>
この酸化反応（以下、「第１酸化反応」とも称呼する。）は、図６に示すように、Ｓｏｏｔを構成する多数の炭素（原子）Ｃのうちで安定状態にあるもの（ラジカル部位を有さないもの）と酸化剤（酸素分子Ｏ_２、酸素原子Ｏ等）との衝突により炭素（原子）Ｃ間の結合（Ｃ=Ｃ結合、共有結合）が崩壊してＳｏｏｔが酸化される反応である。この第１酸化反応に基づく混合気内のＳｏｏｔ濃度[Soot]mixの減少速度（以下、「第１Ｓｏｏｔ酸化速度dmso1/dt」と称呼する。）は、上述の発明の開示の欄で述べた非特許文献２に紹介された「ＮＳＣモデル」を表す下記(17)式に従って表すことができる。

上記(17)式において、Mwcは、炭素の分子量（炭素原子１モルあたりの質量）である。Rtotalは、下記(18)式に従って得られる値である。下記(18)において値ｘは下記(19)式に従って得られる値である。

上記(18)式、(19)式において、Po2は燃焼室内における酸素分圧であり、筒内ガス圧力Pgに後述する混合気内の酸素濃度[O2]mixを乗じることで得ることができる。Ka＝20・exp(-30,000/(R・Tmix))であり、Kb＝4.46×10^-3・exp(-15,200/(R・Tmix))であり、Kt＝1.51×10⁵・exp(-97,000/(R・Tmix))であり、Kz＝21.3・exp(4,100/(R・Tmix))である。ここで、上記値「30,000」，「15,200」，「97,000」はそれぞれ活性化エネルギーである。

上記(17)式において、Ｓｏｏｔ濃度[Soot]mixも、混合気温度Tmixも、上述のように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。従って、第１Ｓｏｏｔ酸化速度dmso1/dtも、上記(17)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

<<低温域での酸化剤による炭素ラジカルの酸化反応>>
この酸化反応（以下、「第２酸化反応」とも称呼する。）は、図７に示すように、Ｓｏｏｔを構成する多数の炭素（原子）のうちでラジカル部位を有するもの（以下、「炭素ラジカル」とも称呼される。）と酸化剤（酸素分子Ｏ_２、酸素原子Ｏ等）との衝突によりＳｏｏｔが酸化される反応（ラジカル反応）である。この第２酸化反応に基づく混合気内のＳｏｏｔ濃度[Soot]mixの減少速度（以下、「第２Ｓｏｏｔ酸化速度dmso2/dt」と称呼する。）は、水酸化ラジカルＯＨの濃度を用いることなく記述された実験式である下記(20)式に従って表すことができる。

上記(20)式において、Aso，αsoは定数、Eso2は活性化エネルギ（定数）である。ここで、この活性化エネルギーEso2は、上記(17)式（より具体的には、(18)式、(19)式）にて使用される上述の活性化エネルギーの値よりも小さい。なお、活性化エネルギーEso2は、発明の開示の欄で述べた非特許文献１に紹介された「廣安モデル」で使用される活性化エネルギーの値よりも小さい。活性化エネルギーEso2は、例えば、「０」であってもよい。

混合気内の酸素濃度[O2]mixは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。従って、第２Ｓｏｏｔ酸化速度dmso2/dtも、上記(20)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

<<高温域での水酸化ラジカルによる炭素の酸化反応>>
この酸化反応（以下、「第３酸化反応」とも称呼する。）は、図８に示すように、Ｓｏｏｔを構成する多数の炭素（原子）Ｃのうちで安定状態にあるもの（ラジカル部位を有さないもの）と水酸化ラジカル（ＯＨ）との衝突により炭素（原子）Ｃ間の結合（Ｃ=Ｃ結合、共有結合）が崩壊してＳｏｏｔが酸化される反応である。この水酸化ラジカル（ＯＨ）は、燃料と酸素分子との反応等に起因して燃焼室内にて発生し得るものである。この第３酸化反応に基づく混合気内のＳｏｏｔ濃度[Soot]mixの減少速度（以下、「第３Ｓｏｏｔ酸化速度dmso3/dt」と称呼する。）は、下記(21)式に従って表すことができる。

上記(21)式において、γohはNeoh’sモデルによる反応確率（定数）であり、Mwohは、水酸化ラジカルＯＨの分子量（水酸化ラジカル１モルあたりの質量）である。Naはアボガドロ数であり、Eso3は活性化エネルギ（定数）である。noh、nsootは、ＯＨの粒子数密度、Ｓｏｏｔの粒子数密度であり、それぞれ、下記(22)式、(23)式に従って表すことができる。(22)式において、[OH]mixは混合気内の水酸化ラジカル濃度である。

上記(21)式は、Tao, F. et al., Combust Flame
136(2004)270-282 にて紹介されたＯＨによるＳｏｏｔの酸化モデルを改良したものである。(21)式の右辺において、「γoh」は反応確率を表し、「Mwc/Na」は１反応あたりのＳｏｏｔ減少量を表し、「√((8・R・Tmix)/(Pg・Mwoh))」はＯＨの速度を表し、「π・(Dsoot/2)²」はＳｏｏｔ１粒子の断面積を表し、「noh・nsoot」は粒子数密度を表す。即ち、「√((8・R・Tmix)/(Pg・Mwoh))・π・(Dsoot/2)²・noh・nsoot」は単位時間、単位体積あたりの衝突数を表している。

上記(21)式において、混合気内の水酸化ラジカル濃度[OH]mixは、後述するように、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。従って、第３Ｓｏｏｔ酸化速度dmso3/dtも、上記(21)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。以上、この反応モデルでは、上記(2)式の右辺第２項のＳｏｏｔ酸化速度dmso/dtの算出にあたり、Ｓｏｏｔの酸化反応として上述の３つの酸化反応が考慮される。

<<温度に対する第１〜第３Ｓｏｏｔ酸化速度の特性>>
ところで、第１〜第３Ｓｏｏｔ酸化速度は、混合気温度Tmixに対して、図９において破線、２点鎖線、１点鎖線で示した特性をそれぞれ有している。図９において、実線は、第１〜第３Ｓｏｏｔ酸化速度の和の混合気温度Tmixに対する特性を示している。

図９に示すように、第２酸化反応は、比較的高温雰囲気下（例えば、１３００Ｋ以上）のみならず、比較的低温雰囲気下（例えば、１３００Ｋ未満）でも無視できない程度に発生し、第２Ｓｏｏｔ酸化速度dmso2/dtは、比較的低温雰囲気下から無視できない程度に十分大きい値に算出されて、混合気温度Tmixに増加に応じて緩やかに増大していく。

一方、第１、第３酸化反応は、比較的低温雰囲気下（例えば、１３００Ｋ未満）では殆ど発生せず、第１、第３Ｓｏｏｔ酸化速度dmso1/dt，dmso3/dtは、１３００Ｋ未満では無視できる程度に十分小さい値に算出され、１３００Ｋ以上では、混合気温度Tmixに増加に応じて急激に増大していく。

以上より、比較的低温雰囲気下（例えば、１３００Ｋ未満）では、図１０に示すように、第２酸化反応が支配的となり、第１、第３酸化反応が殆ど発生しない。一方、比較的高温雰囲気下（例えば、１３００Ｋ以上）では、図１１に示すように、第１〜第３酸化反応の全てが複合的に発生する。なお、これは、第２酸化反応についての活性化エネルギーEso2（(20)式を参照）が第１、第３酸化反応についての活性化エネルギー（(18)式、(19)式、並びに、(21)式のEso3を参照）よりも小さいことに基づく。

以上のことに鑑み、本装置では、混合気温度Tmixが１３００Ｋ未満の場合（具体的には、燃料噴射から、燃焼反応により増大していく混合気温度Tmixが１３００Ｋに達するまでの期間）、第２酸化反応のみが考慮された下記(24)式に従って上記(2)式の右辺第２項のＳｏｏｔ酸化速度dmso/dtが逐次算出される。これにより、Ｓｏｏｔ酸化速度dmso/dtの計算精度を下げることなく、上記(17)式、(21)式に表される第１、第３Ｓｏｏｔ酸化速度dmso1/dt，dmso3/dtの計算そのものを省略してＳｏｏｔ酸化速度dmso/dtの計算に要する負荷を軽減できる。

一方、混合気温度Tmixが１３００Ｋ以上の場合（具体的には、混合気温度Tmixが上昇しながら１３００Ｋに達した時点以降、ピーク温度を経て１３００Ｋに復帰するまでの期間）では、第１〜第３酸化反応の全てが考慮された下記(25)式に従って上記(2)式の右辺第２項のＳｏｏｔ酸化速度dmso/dtが逐次算出される。これにより、複合的に発生する全てのＳｏｏｔ酸化反応が考慮されてＳｏｏｔ酸化速度dmso/dtが精度良く計算できる。以上、上記(2)式のＳｏｏｔ酸化速度dmso/dtは、上記(24)式、(25)式に従って、噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求めることができる。

以上より、上記(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（Ｓｏｏｔの発生量）が推定される。ここで、燃料からＳｏｏｔが発生する反応過程を、燃料から前駆物質Ｐが発生する反応過程と、前駆物質ＰからＳｏｏｔが発生する反応過程の２段階とすることによる効果については、特開２００７−４６４７７号公報に詳細に記載されている。

以上、図５に示した反応モデルを用いたＳｏｏｔ発生量（混合気内のＳｏｏｔ濃度[Soot]mix）の推定について説明した。次に、係るＳｏｏｔ濃度[Soot]mixの推定に必要となる、混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、混合気内の酸素濃度[O2]mix、混合気内の水酸化ラジカル濃度[OH]mix、及び筒内ガス圧力Pgの推定方法について説明する。以下、先ず、これらの物理量の推定に必要となる、任意の噴射後経過時間ｔにおける空気過剰率λの取得方法について説明する。

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

このように定義される空気過剰率λは、例えば、日本機械学会論文集 25-156(1959年)，820ページ 「ディーゼル機関の噴霧到達距離に関する研究」 和栗雄太郎，藤井勝，網谷竜夫，恒屋礼次郎 （以下、「非特許文献３」と称呼する。）にて紹介された実験式である下記(27)式、及び下記(28)式に基づいて噴射後経過時間ｔの関数として求めることができる。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<混合気内の水酸化ラジカル濃度[OH]mixの取得>
次に、混合気内の水酸化ラジカル濃度（質量濃度）[OH]mixを取得する方法について説明する。水酸化ラジカルＯＨは、燃料と酸素分子との反応に起因して発生し得、ＯＨの発生速度（具体的には、混合気内のＯＨ量の変化速度、混合気内のＯＨ濃度の変化速度等）は、酸素濃度に大きく依存する。従って、ＯＨの発生速度d[OH]mix/dtは、混合気内の酸素濃度[O2]mixを引数とするＯＨの発生速度d[OH]mix/dtを求めるための関数funcを用いて下記(47)式に従って求めることができる。

従って、上記(47)式にて得られるＯＨの発生速度d[OH]mix/dtを噴射後経過時間ｔ＝０から微小時間Δｔ毎に逐次求め、これらの値を時間で積分（積算）していくことで、任意の噴射後経過時間ｔにおける混合気内のＯＨ濃度[OH]mixを求めることができる。

<筒内ガス圧力Pgの取得>
次に、筒内ガス圧力Pgを取得する方法について説明する。筒内ガス圧力Pgは、筒内ガスの状態がＩＶＣ以降断熱変化するとの仮定のもと、下記(48)式に従って求めることができる。下記(48)式において、上述したように、PgivcはＩＶＣにおける筒内ガス圧力であり、Vg(CAivc)はＩＶＣにおけるクランク角度CAに対応する筒内容積である。また、Vg(CA)は現時点（即ち、噴射後経過時間ｔ）における上記筒内容積Vg(CA)である。以上、Ｓｏｏｔ濃度[Soot]mixの推定に必要となる、空気過剰率λ、混合気温度Tmix、混合気内の燃料濃度[Fuel]mix、混合気内の酸素濃度[O2]mix、混合気内の水酸化ラジカル濃度[OH]mix、及び筒内ガス圧力Pgの推定方法について説明した。

（燃料噴射制御の概要）
本装置は、上述した反応モデルを用いた混合気内のＳｏｏｔ濃度[Soot]mixの推定に係わる計算を筒内ガスの量が確定するＩＶＣの直後に開始し、燃料噴射開始時期CAinjの到来前に上記Ｓｏｏｔ発生質量Ｓｏｏｔnewの推定を完了する。

そして、本装置は、目標Ｓｏｏｔ発生質量Ｓｏｏｔterを機関の運転状態から求め、上記推定されたＳｏｏｔ発生質量Ｓｏｏｔnewが目標Ｓｏｏｔ発生質量Ｓｏｏｔterに対して十分に大きい場合、Ｓｏｏｔ発生質量Ｓｏｏｔnewが小さくなるように燃料噴射圧力をフィードバック制御する。

具体的には、本装置は、Ｓｏｏｔ発生質量Ｓｏｏｔnewから目標Ｓｏｏｔ発生質量Ｓｏｏｔterを減じた値が所定量よりも大きいとき、燃料噴射開始圧力を基本燃料噴射圧力Pcrbaseよりも所定量だけ高くする。これにより、Ｓｏｏｔ発生質量Ｓｏｏｔnewが小さくなる方向に制御される。以上が、燃料噴射制御の概要である。

（実際の作動）
次に、上記のように構成された内燃機関のすす発生量推定装置の実際の作動について説明する。

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

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

続いて、ＣＰＵ６１はステップ１２１５に進んで、上記設定されたＩＶＣ時筒内ガス圧力Pgivcと、上記設定されたＩＶＣ時筒内ガス温度Tgivcと、上記(30)式とに基づいて筒内ガスの全質量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と、上記(43)式とに基づいて燃料蒸気温度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)と、上記(29)式とに基づいて燃料噴射開始時筒内ガス圧力Pg0を求める。

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

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

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

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

次いで、ＣＰＵ６１はステップ１３３５に進んで、混合気内の前駆物質Ｐ濃度[Ｐ]mix、混合気内のＳｏｏｔ濃度[Soot]mix、混合気内の酸素濃度[O2]mix、混合気内の水酸化ラジカル濃度[OH]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)と、上記(48)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における筒内ガス圧力Pgを求める。

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

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

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

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

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

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

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

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

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

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

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

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

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

次に、ＣＰＵ６１は図１５のルーチンに進み、各種濃度の算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ１５０５に進み、上記ステップ１４６５にて求めた微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による燃料消費量qrと、ステップ１２４０にて求めた筒内ガス理論空燃比stoichと、上記(45)式とに基づいて、微小時間Δｔ（噴射後経過時間(ｔ−Δｔ)〜ｔの間）において混合気内で発生する化学反応による筒内ガス消費量grを求める。

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

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

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

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

次いで、ＣＰＵ６１はステップ１５３０に進み、上記ステップ１５２０にて求めた混合気内の酸素濃度[O2]mixと、上記(47)式とに基づいて噴射後経過時間ｔ（従って、クランク角度CA）における水酸化ラジカルＯＨの発生速度d[OH]mix/dtを求め、続くステップ１５３５にて、水酸化ラジカル濃度[OH]mixを、その時点での水酸化ラジカル濃度[OH]mix（現時点では、ステップ１３３５の処理により「０」）に、上記求めたＯＨの発生速度d[OH]mix/dtに微小時間Δｔを乗じた値「d[OH]mix/dt・Δｔ」を加えた値に更新する。これにより、噴射後経過時間ｔ（従って、クランク角度CA）における混合気内の水酸化ラジカル濃度[OH]mixが求められる。

次に、ＣＰＵ６１はステップ１５４０に進んで、上記ステップ１５１５にて求めた混合気内の燃料濃度[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を求める。ここで、ステップ１５６０は前記「すす生成速度算出手段」に対応する。

続いて、ＣＰＵ６１はステップ１５６５に進み、上記ステップ１５４０にて求めた前駆物質Ｐの生成速度dmpf/dtと、上記ステップ１５４５にて求めた前駆物質Ｐの分解速度dmpd/dtと、上記ステップ１５６０にて求めたＳｏｏｔの生成速度dmsf/dtと、上記(1)式とに基づいて前駆物質Ｐの発生速度d[Ｐ]mix/dtを求める。

続いて、ＣＰＵ６１はステップステップ１５７０に進み、混合気内の前駆物質Ｐ濃度[Ｐ]mixを、その時点での値（現時点では、ステップ１３３５の処理により「０」）に、上記求めた前駆物質Ｐの発生速度d「Ｐ」mix/dtに微小時間Δｔを乗じた値「d[Ｐ]mix/dt・Δｔ」を加えた値に設定・更新する。これにより、噴射後経過時間ｔ＝Δｔ（従って、クランク角度CA＝CAinj＋ΔCA）における混合気内の前駆物質Ｐ濃度[Ｐ]mixが求められる。

次に、ＣＰＵ６１は図１６のルーチンに進み、Ｓｏｏｔ濃度の算出のための処理を開始する。具体的には、ＣＰＵ６１は先ずステップ１６０５に進み、上記混合気内のＳｏｏｔ濃度[Soot]mixと、上記混合気内の酸素濃度［０２］mixと、上記混合気温度Tmixと、上記(20)式とに基づいて第２Ｓｏｏｔ酸化速度dmso2/dtを求める。ここで、ステップ１６０５は前記「第１すす酸化速度算出手段」に対応する。

次いで、ＣＰＵ６１はステップ１６１０に進んで、上記混合気温度Tmix（ステップ１４８０にて更新されている）が１３００Ｋ未満であるか否かを判定し、「Ｙｅｓ」と判定される場合、ステップ１６１５に進んで、Ｓｏｏｔ酸化速度dmso/dtを上記求めた第２Ｓｏｏｔ酸化速度dmso2/dtと等しい値に設定する。ここで、ステップ１６１５は前記「すす酸化速度算出手段」の一部に対応する。

一方、ステップ１６１０にて「Ｎｏ」と判定される場合、ＣＰＵ６１はステップ１６２０に進んで、値ｘを上記(19)式に従って算出し、続くステップ１６２５にて値Rtotalを上記(18)式に従って算出し、続くステップ１６３０にて、上記混合気内のＳｏｏｔ濃度[Soot]mixと、上記求めた値Rtotalと、上記(17)式とに基づいて第１Ｓｏｏｔ酸化速度dmso1/dtを求める。ここで、ステップ１６３０は前記「第２すす酸化速度算出手段」に対応する。

次いで、ＣＰＵ６１はステップ１６３５に進み、ステップ１５３５にて更新されている水酸化ラジカル濃度[OH]mixと、上記(22)式とに基づいてＯＨの粒子数密度nohを求め、続くステップ１６４０にて、上記Ｓｏｏｔ濃度[Soot]mixと、上記(23)式とに基づいてＳｏｏｔの粒子数密度nsootを求め、続くステップ１６４５にて、上記混合気温度Tmixと、上記筒内ガス圧力Pgと、上記求めた粒子数密度noh，nsootと、上記(21)式とに基づいて第３Ｓｏｏｔ酸化速度dmso3/dtを求める。ここで、ステップ１６４５は前記「第３すす酸化速度算出手段」に対応する。

そして、Ｓｏｏｔ酸化速度dmso/dtを、ステップ１６３０にて求めた第１Ｓｏｏｔ酸化速度dmso1/dtと、ステップ１６０５にて求めた第２Ｓｏｏｔ酸化速度dmso2/dtと、ステップ１６４５にて求めた第３Ｓｏｏｔ酸化速度dmso3/dtとの和の値に設定する。ここで、ステップ１６５０は前記「すす酸化速度算出手段」の一部に対応する。

次に、ＣＰＵ６１はステップ１６５５に進み、ステップ１５６０にて求めたＳｏｏｔ生成速度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，[OH]mix，[Ｐ]mix，[Soot]mix。ステップ１５１５、１５２０、１５３５、１５７０、１６６０を参照。）が算出される。

そして、図１７のステップ１７０５の判定において「Ｎｏ」と判定される毎に、図１４のステップ１４０５〜図１７のステップ１７０５の処理が繰り返し実行されていく。これにより、図１７のステップ１７０５の判定において「Ｎｏ」と判定される限りにおいて、混合気の空気過剰率λ、混合気温度Tmix、及び各種濃度（[Fuel]mix，[O2]mix，[OH]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をもって噴射される。

以上、説明したように、本発明によるすす発生量推定装置の実施形態によれば、生成されたＳｏｏｔの酸化速度dmso/dtを算出するにあたり、「高温域での酸化剤による炭素の酸化反応（第１酸化反応、第１Ｓｏｏｔ酸化速度dmso1/dt）」、「低温域での酸化剤による炭素ラジカルの酸化反応（第２酸化反応、第２Ｓｏｏｔ酸化速度dmso2/dt）」、及び「高温域での水酸化ラジカルによる炭素の酸化反応（第３酸化反応、第３Ｓｏｏｔ酸化速度dmso3/dt）」が考慮される。従って、複合的に発生する全てのＳｏｏｔ酸化反応が考慮されてＳｏｏｔ酸化速度dmso/dtが精度良く計算できる。この結果、精度良く計算され得るＳｏｏｔ酸化速度dmso/dtに基づいてＳｏｏｔの発生量も精度良く計算できる。

加えて、混合気温度Tmixが１３００Ｋ未満の場合、殆ど反応が発生しない第１、第３酸化反応が考慮されずに第２酸化反応のみが考慮されてＳｏｏｔの酸化速度dmso/dtが計算される。これにより、混合気温度Tmixが１３００Ｋ未満の場合において、Ｓｏｏｔ酸化速度dmso/dtの計算精度を下げることなく、第１、第３Ｓｏｏｔ酸化速度dmso1/dt，dmso3/dtの計算そのものを省略してＳｏｏｔ酸化速度dmso/dtの計算に要する負荷を軽減できる。

本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記実施形態では、混合気温度Tmixが１３００Ｋ未満の場合、第１、第３酸化反応が考慮されずに第２酸化反応のみが考慮されてＳｏｏｔの酸化速度dmso/dtが計算されているが、混合気温度Tmixが１３００Ｋ未満の場合も、混合気温度Tmixが１３００Ｋ以上の場合と同様、第１〜第３酸化反応の全てを考慮してＳｏｏｔの酸化速度dmso/dtを計算してもよい。

また、上記実施形態においては、混合気温度Tmixが１３００Ｋ以上の場合、第１〜第３酸化反応の全てを考慮してＳｏｏｔの酸化速度dmso/dtが計算されているが、混合気温度Tmixが１３００Ｋ以上の場合において、第１、第３酸化反応の何れか一方を考慮せずにＳｏｏｔの酸化速度dmso/dtを計算してもよい。

また、上記実施形態においては、上記第１酸化反応を表すモデルとして非特許文献２に紹介されている「ＮＳＣモデル」が使用されているが、上記第１酸化反応を表すモデルとして、例えば、非特許文献１に紹介されている「廣安（酸化）モデル」が使用されてもよい。

また、上記実施形態においては、上記第２酸化反応を表すモデルとして実験式である上記(20)式が使用されているが、上記第２酸化反応を表すモデルとして、例えば、Frenklach.,第23回国際燃焼シンポジウム(1990)P.1559-1566に紹介されている、水酸化ラジカルＯＨの濃度を用いたモデルを使用してもよい。

また、上記実施形態においては、Ｓｏｏｔの生成速度dmsf/dtを表すモデルとして前駆物質Ｐを考慮した図５に示したモデルが使用されているが、Ｓｏｏｔの生成速度dmsf/dtを表すモデルとして、例えば、非特許文献１に紹介されている「廣安（生成）モデル」が使用されてもよい。

また、上記実施形態においては、クランク角度CAがＴＤＣ以降であって、且つ、混合気温度TmixがＳｏｏｔの生成反応限界温度Tminより低くなった時点以降、ＣＰＵ６１の計算負荷低減のために前駆物質Ｐの発生速度、及びＳｏｏｔの発生速度の時間積分処理が終了するように構成されているが、これに加え、クランク角度CAがＴＤＣ以前であって、且つ、混合気温度TmixがＳｏｏｔの生成反応限界温度Tminを超えるまでの間も、前駆物質Ｐの発生速度、及びＳｏｏｔの発生速度の時間積分処理を行わないように構成してもよい。これにより、Ｓｏｏｔの発生量の算出に係わる不必要な計算が省略でき、より一層ＣＰＵ６１の計算負荷を低減することができる。

加えて、上記実施形態においては、種々の反応速度を算出する式（具体的には、上記(3)式、(4)式、(13)式、(16)式の何れかにおいて、筒内ガス圧力Pgに係る項を省略してもよい。

本発明の実施形態に係る内燃機関のすす発生量推定装置を４気筒内燃機関（ディーゼル機関）に適用したシステム全体の概略構成図である。 或る一つの気筒のシリンダ内(筒内)に吸気マニホールドからガスが吸入され、筒内に吸入されたガスが排気マニホールドへ排出される様子を模式的に示した図である。 （ａ）は、燃料が一時に噴射された時点での同燃料の様子を模式的に示した図であり、（ｂ）は、筒内ガスと混ざり合いながら混合気となって円錐状に拡散していく燃料の様子を模式的に示した図である。 前駆物質の例としてＰＡＨ（多環芳香族炭化水素類）を挙げた場合における、燃料から前駆物質を経由してＳｏｏｔが発生する反応過程を模式的に示した図である。 図１に示したすす発生量推定装置が使用する、燃料からＳｏｏｔが発生する反応過程において前駆物質が経由される反応モデルを表した図である。 Ｓｏｏｔの酸化に関して「高温域での酸化剤による炭素の酸化反応（第１酸化反応）」を説明するための図である。 Ｓｏｏｔの酸化に関して「低温域での酸化剤による炭素ラジカルの酸化反応（第２酸化反応）」を説明するための図である。 Ｓｏｏｔの酸化に関して「高温域での水酸化ラジカルによる炭素の酸化反応（第３酸化反応）」を説明するための図である。 第１〜第３Ｓｏｏｔ酸化速度と、混合気温度との関係を示したグラフである。 混合気温度が１３００Ｋ未満の場合において、第２酸化反応のみを考慮する様子を示した図である。 混合気温度が１３００Ｋ以上の場合において、第１〜第３酸化反応の全てを考慮する様子を示した図である。 図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第１番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第２番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第３番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第４番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第５番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する混合気温度等、及びＳｏｏｔ発生量の算出を行うためのルーチンの第６番目の部分を示したフローチャートである。 図１に示したＣＰＵが図１２に示したルーチンを実行する際に参照する指令燃料噴射量を決定するためのテーブルである。 図１に示したＣＰＵが図１２に示したルーチンを実行する際に参照する燃料燃料噴射時期を決定するためのテーブルである。 図１に示したＣＰＵが図１２に示したルーチンを実行する際に参照する基本燃料噴射圧力を決定するためのテーブルである。 図１に示したＣＰＵが実行する燃料噴射制御を行うためのルーチンを示したフローチャートである。

符号の説明

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

Claims (10)

1. 内燃機関の燃焼室内において燃料の反応に起因して発生するすすの発生量を推定するすす発生量推定手段を備えた内燃機関のすす発生量推定装置において、
   前記すす発生量推定手段は、
   前記燃焼室内の温度が所定温度未満の場合において、生成されたすすを構成する炭素であってラジカル部位を有するものと酸化剤との衝突により前記生成されたすすが酸化される反応を考慮して前記すすの発生量を推定するように構成された内燃機関のすす発生量推定装置。
2. 請求項１に記載の内燃機関のすす発生量推定装置において、
   前記すす発生量推定手段は、
   前記燃料の反応に起因して前記すすが生成される速度であるすす生成速度を算出するすす生成速度算出手段と、
   前記生成されたすすが酸化される速度であるすす酸化速度を算出するすす酸化速度算出手段と、
   前記すす生成速度と前記すす酸化速度とに基づいて前記すすが発生する速度であるすす発生速度を算出するすす発生速度算出手段と、
   を備え、
   前記すす発生速度に基づいて前記すすの発生量を推定するように構成されていて、
   前記すす酸化速度算出手段は、
   前記ラジカル部位を有する炭素と前記酸化剤との衝突に起因して前記生成されたすすが酸化される速度である第１すす酸化速度を算出する第１すす酸化速度算出手段を備え、
   前記燃焼室内の温度が前記所定温度未満の場合において、前記第１すす酸化速度に基づいて前記すす酸化速度を算出するように構成された内燃機関のすす発生量推定装置。
3. 請求項２に記載の内燃機関のすす発生量推定装置において、
   前記すす酸化速度算出手段は、
   前記生成されたすすを構成する炭素であって安定状態にあるものと酸化剤との衝突に起因して前記生成されたすすが酸化される速度である第２すす酸化速度を算出する第２すす酸化速度算出手段を備え、
   前記燃焼室内の温度が前記所定温度未満の場合では、前記第１すす酸化速度のみに基づいて前記すす酸化速度を算出するとともに、前記燃焼室内の温度が前記所定温度以上の場合では、前記第１すす酸化速度と前記第２すす酸化速度とに基づいて前記すす酸化速度を算出するように構成された内燃機関のすす発生量推定装置。
4. 請求項２に記載の内燃機関のすす発生量推定装置において、
   前記すす酸化速度算出手段は、
   前記生成されたすすを構成する炭素であって安定状態にあるものと燃料の反応に起因して発生する水酸化ラジカルとの衝突に起因して前記生成されたすすが酸化される速度である第３すす酸化速度を算出する第３すす酸化速度算出手段を備え、
   前記燃焼室内の温度が前記所定温度未満の場合では、前記第１すす酸化速度のみに基づいて前記すす酸化速度を算出するとともに、前記燃焼室内の温度が前記所定温度以上の場合では、前記第１すす酸化速度と前記第３すす酸化速度とに基づいて前記すす酸化速度を算出するように構成された内燃機関のすす発生量推定装置。
5. 請求項２に記載の内燃機関のすす発生量推定装置において、
   前記すす酸化速度算出手段は、
   前記生成されたすすを構成する炭素であって安定状態にあるものと酸化剤との衝突に起因して前記生成されたすすが酸化される速度である第２すす酸化速度を算出する第２すす酸化速度算出手段と、
   前記生成されたすすを構成する炭素であって安定状態にあるものと燃料の反応に起因して発生する水酸化ラジカルとの衝突に起因して前記生成されたすすが酸化される速度である第３すす酸化速度を算出する第３すす酸化速度算出手段と、
   を備え、
   前記燃焼室内の温度が前記所定温度未満の場合では、前記第１すす酸化速度のみに基づいて前記すす酸化速度を算出するとともに、前記燃焼室内の温度が前記所定温度以上の場合では、前記第１すす酸化速度と前記第２すす酸化速度と前記第３すす酸化速度とに基づいて前記すす酸化速度を算出するように構成された内燃機関のすす発生量推定装置。
6. 請求項２乃至請求項５の何れか一項に記載の内燃機関のすす発生量推定装置において、
   前記第１すす酸化速度算出手段は、
   燃料の反応に起因して発生する水酸化ラジカルの濃度を用いることなく記述される関係に基づいて前記第１すす酸化速度を算出するように構成された内燃機関のすす発生量推定装置。
7. 請求項６に記載の内燃機関のすす発生量推定装置において、
   前記第１すす酸化速度算出手段は、
   前記燃焼室内のすすの濃度、前記燃焼室内の酸素の濃度、前記燃焼室内の温度のべき乗、及び前記燃焼室内の温度に基づいて決定される値を指数とする指数関数、の積により表される関係に基づいて前記第１すす酸化速度を算出するように構成された内燃機関のすす発生量推定装置。
8. 請求項２乃至請求項７の何れか一項に記載の内燃機関のすす発生量推定装置において、
   前記すす生成速度算出手段は、前記燃焼室への燃料噴射後において時々刻々と変化していく前記すす生成速度を逐次算出し、
   前記すす酸化速度算出手段は、前記燃料噴射後において時々刻々と変化していく前記すす酸化速度を逐次算出し、
   前記すす発生速度算出手段は、前記燃料噴射後において時々刻々と変化していく前記すす生成速度と前記すす酸化速度とに基づいて前記燃料噴射後において時々刻々と変化していく前記すす発生速度を逐次算出するように構成されていて、
   前記すす発生量推定手段は、
   前記燃料噴射後において時々刻々と変化していく前記すす発生速度を逐次積算していくことで前記すすの発生量を推定するように構成された内燃機関のすす発生量推定装置。
9. 請求項１乃至請求項８の何れか一項に記載の内燃機関のすす発生量推定装置において、
   前記所定温度は、１３００Ｋである内燃機関のすす発生量推定装置。
10. 請求項１乃至請求項９の何れか一項に記載の内燃機関のすす発生量推定装置において、
    前記すす発生量推定手段は、圧縮比が１６以下の内燃機関の燃焼室内において発生するすすの発生量を推定するように構成された内燃機関のすす発生量推定装置。

Images