JP2007247487A - 内燃機関の混合気の不均一度取得装置、及び混合気状態取得装置 - Google Patents

Abstract

【課題】 比較的簡易な計算により混合気内における燃料濃度の不均一度を精度良く取得することができる内燃機関の混合気の不均一度取得装置を提供すること。
【解決手段】 この装置では、混合気は、燃料質量分率が異なる球形のガス塊の集合体であると仮定される。燃料噴射開始時点以降、微小時間Δtの経過毎に、その時点で存在しているガス塊のうちの任意の２つの組み合わせの全てについて「衝突反応」が順に行われる。これにより、混合気内のガス塊の質量分布が微小時間Δt毎に取得・更新され、混合気の不均一度が取得・更新されていく。「衝突反応」とは、「燃料質量分率が異なる２つのガス塊の衝突により衝突した２つのガス塊のそれぞれの一部分が混合し、同混合した前記各一部分が対応するガス塊からそれぞれ離脱して前記２つのガス塊の燃料質量分率とは異なる燃料質量分率を有する他のガス塊の一部又は全部となる反応」をいう。
【選択図】 図９

Description

本発明は、内燃機関の燃焼室内にて同燃焼室内に吸入されている筒内ガスと噴射された燃料とが混ざり合って形成されていく混合気の不均一度を取得する内燃機関の混合気の不均一度取得装置、及びその不均一度を利用して混合気の状態を取得する混合気状態取得装置に関する。

内燃機関の燃焼制御において、燃焼騒音や発生トルク等を精度良く予測することは重要である。燃焼騒音や発生トルク等を精度良く予測するためには、燃焼室内におけるクランク角度に対する熱発生率（熱発生速度）を精度良く予測する必要がある。この熱発生率は燃焼室内にて形成される混合気内での燃料（燃料蒸気）の反応速度に依存し、燃料の反応速度は混合気の状態（燃料濃度、温度等）に大きく依存する。従って、熱発生率を精度良く予測するためには、混合気の状態を精度良く予測する必要がある。混合気の状態を予測する手法としては種々のものが知られている（例えば、下記特許文献１を参照）。
特開２００１−２５４６４５号公報

ところで、燃焼室内にて形成されていく混合気を微視的に見た場合、混合気内における燃料濃度（燃料質量濃度）は、混合気が占める領域内の位置に応じて異なる。従って、混合気内の燃料濃度に大きく影響される熱発生率も混合気が占める領域内の位置に応じて異なる。換言すれば、燃焼室内にて形成されていく混合気内の燃料濃度が不均一となるから、熱発生率も不均一になる。この結果、燃焼室内全体における熱発生率は混合気内の燃料濃度の不均一度に大きく影響され得ると考えられる。

以上のことから、熱発生率（従って、燃焼騒音、発生トルク等）を精度良く予測するためには、混合気内における燃料濃度の不均一度を精度良く取得することが重要であると考えられる。しかしながら、上記文献では、混合気内における燃料濃度の不均一度については一切記載されていない。他方、混合気内における燃料濃度の不均一度を精度良く取得する手法は、従来より種々の論文等にて紹介されてきている。

しかしながら、何れの手法も、例えば、確率密度関数等を用いた複雑な微分方程式等を解く必要があり膨大な計算負荷を伴うものであった。従って、複雑な微分方程式等を解く必要のない比較的簡易な計算により、混合気内における燃料濃度の不均一度を精度良く取得する手法の到来が望まれていたところである。

以上のことから、本発明の目的は、比較的簡易な計算により混合気内における燃料濃度の不均一度を精度良く取得することができる内燃機関の混合気の不均一度取得装置を提供することにある。

本発明に係る混合気の不均一度取得装置では、燃焼室内にて同燃焼室内に吸入されている筒内ガスと噴射された燃料とが混ざり合って混合気が形成されていく過程において、「衝突反応」が混合気内で繰り返し行われると仮定される。「衝突反応」とは、燃料質量割合が異なる２つのガス塊の衝突により前記２つのガス塊のそれぞれの一部分が混合し同混合した前記各一部分が対応する前記ガス塊からそれぞれ離脱して前記２つのガス塊の燃料質量割合とは異なる燃料質量割合を有する他のガス塊の一部又は全部となる反応である。

そして、係る仮定に基づいて混合気内における燃料濃度（燃料質量濃度）の不均一度を取得する不均一度取得手段が備えられている。このように、本発明では、混合気は、燃料質量割合が異なるガス塊の集合体であると考える。そして、２つ以上のガス塊のうちの任意の２つの組み合わせについて上記衝突反応が順に繰り返されることで混合気内におけるそれぞれのガス塊の量（質量、体積等）の分布が逐次変化していく。本発明では、例えば、係るガス塊の量の分布に基づいて混合気内における燃料濃度の不均一度（以下、単に「不均一度」と云うこともある。）が取得される。

なお、本発明では、「燃料質量割合」における「燃料」とは、燃料噴射により新たに燃焼室内に供給された燃料（燃料蒸気）を意味していて、筒内ガスに予め含まれる燃料（例えばＥＧＲ等により筒内ガスに含まれる未燃燃料）を含まないと扱うことができる。この場合、筒内ガスそのもの（燃料蒸気を含んでいない）の燃料質量割合は「０」となる。一方、燃料噴射による燃料蒸気そのもの（筒内ガスを含んでいない）の燃料質量割合は「１」となる。

具体的には、前記不均一度取得手段は、燃料質量割合が「０」の前記筒内ガスの塊、及び燃料質量割合が「１」の燃料蒸気の塊を前記衝突反応の対象となるガス塊として扱うように構成される。即ち、本発明では、混合気の一部として、筒内ガスそのもの（の塊）、及び燃料蒸気そのもの（の塊）が含まれる。なお、筒内ガスに予め含まれる燃料（例えばＥＧＲ等により筒内ガスに含まれる未燃燃料）を考慮して筒内ガスの塊の燃料質量割合を決定してもよい。

燃焼室内にて混合気が形成されていく過程では、時間の経過とともに筒内ガスが混合気の一部として混合気内に継続的に取り込まれていくとともに、燃料が噴射されている期間内では噴射された燃料（燃料蒸気）が混合気の一部として扱われていく。従って、混合気内での筒内ガスの塊（即ち、燃料質量割合が「０」のガス塊）の量における増加分は、時間の経過とともに新たに取り込まれていく筒内ガスの量に依存する。また、混合気内の燃料蒸気の塊（即ち、燃料質量割合が「１」のガス塊）の量における増加分は、時間の経過とともに新たに噴射されていく燃料の量に依存する。

前記不均一度取得手段は、前記混合気内における燃料質量割合が異なる２つ以上の前記ガス塊のうちの任意の２つの組み合わせの全てについて前記衝突反応が順に行われるものと仮定して前記不均一度を取得するように構成されることが好適である。これにより、例えば、混合気内におけるそれぞれのガス塊の量の分布（即ち、不均一度）を精度良く求めることができる。

また、上記本発明に係る不均一度取得装置においては、前記衝突した２つのガス塊の前記一部分に含まれるそれぞれの燃料の質量の和が前記２つのガス塊の前記一部分のそれぞれの質量の和に前記他のガス塊の燃料質量割合を乗じた値に等しいとの燃料についての質量保存則を利用して、前記他のガス塊の燃料質量割合を特定する燃料質量割合特定手段を更に備えることが好適である。ここで、「ガス塊の一部分に含まれる燃料の質量」は、ガス塊の一部分の質量にそのガス塊の燃料質量割合を乗じることで求めることができる。

このように、上記衝突反応の前後に亘って混合・離脱に係わった燃料についての質量保存則が成立すると仮定すれば、衝突反応により量が増加するガス塊である上記「他のガス塊」の燃料質量割合を簡易、且つ容易に特定することができる。

上記本発明に係る不均一度取得装置は、前記混合気内における燃料質量割合が異なる２つ以上の前記ガス塊の量をそれぞれ取得するガス塊量取得手段と、前記衝突した２つのガス塊の前記それぞれの一部分が混合する体積である混合体積を取得する混合体積取得手段とを更に備え、前記ガス塊量取得手段は、前記取得された混合体積を用いて前記衝突した２つのガス塊の量における前記一部分に相当するそれぞれの減少分と前記他のガス塊の量における増加分とを算出することで、前記２つのガス塊のそれぞれの量、及び前記他のガス塊の量を更新・取得するように構成されていて、前記不均一度取得手段は、前記更新・取得されたそれぞれのガス塊の量の分布に基づいて前記不均一度を取得するように構成される。ここにおいて、ガス塊の「量」とは、例えば、ガス塊の質量、体積等である。

これによれば、衝突反応により２つのガス塊が混合する体積である混合体積が求められる。この混合体積は、衝突反応の対象となる２つのガス塊からそれぞれ離脱して上記他のガス塊に移行する上記一部分のそれぞれの体積と等しい。上記一部分のそれぞれの量（質量）は、混合体積を用いて（例えば、上記混合体積にガス塊の密度を乗じることで）求めることができる。

加えて、上記一部分のそれぞれの量は上記２つのガス塊の量におけるそれぞれの減少分に等しく、上記一部分のそれぞれの量の和は上記他のガス塊の量における増加分に等しい。以上のことから、上記構成によれば、上記混合体積を用いることで、上記衝突反応による上記２つのガス塊の量におけるそれぞれの減少分及び上記他のガス塊の量における増加分を簡易、且つ容易に求めることができる。

従って、上記２つのガス塊のそれぞれの量、及び上記他のガス塊の量が簡易、且つ容易に更新され得るから、混合気内におけるそれぞれのガス塊の量の分布が簡易、且つ容易に更新され得る。この結果、係るガス塊の量の分布に基づいて不均一度を簡易、且つ容易に求めることができる。

なお、前記不均一度としては、例えば、前記更新・取得されたそれぞれのガス塊の量の分布に基づいて得られる燃料質量割合についての標準偏差、分散等を使用することが好適である。

上記混合体積取得手段が備えられている場合、混合体積取得手段は、前記衝突した２つのガス塊のそれぞれの量（体積）のうち小さい方に基づいて前記混合体積を取得するように構成されることが好適である。

衝突した２つのガス塊が混合する体積である上記混合体積は、衝突した２つのガス塊が係わりあった体積に等しいと考えることができる。２つのガス塊が係わりあった体積は、２つのガス塊のうち体積が小さい方のガス塊（以下、「小ガス塊」と称呼する。）全体が体積が大きい方のガス塊（以下、「大ガス塊」と称呼する。）の内部を通過する体積（以下、「通過体積」と称呼する。）に関連すると考えられる。

この通過体積は、小ガス塊の量（体積、具体的には、例えば、ガス塊を球形と考える場合には小ガス塊の球の断面積）に大きく依存すると考えられる。従って、上記のように構成すれば、通過体積を求めることができるから、通過体積に関連する混合体積を求めることができる。

この場合、前記混合体積取得手段は、前記大ガス塊の体積の前記混合気全体の体積に対する割合に基づいて前記混合体積を取得するように構成されると好ましい。小ガス塊が大ガス塊の内部を通過する割合（確率）は、大ガス塊の体積の混合気全体の体積に対する割合に等しいと仮定することができる。係る仮定によれば、混合体積は、通過体積に上記割合を乗じた値に等しく（或いは、応じた値に）なる。上記構成は係る知見に基づくものである。これによれば、混合体積をより一層精度良く求めることができる。

また、上記混合体積取得手段が備えられている場合、本発明に係る不均一度取得装置は、前記混合気内における燃料質量割合が異なる２つ以上の前記ガス塊の速度をそれぞれ取得するガス塊速度取得手段を更に備え、前記混合体積取得手段は、前記衝突した２つのガス塊の速度の差に基づいて前記混合体積を取得するように構成されることが好適である。

上記通過体積は、小ガス塊の量（例えば、小ガス塊の球の断面積）に加えて、小ガス塊と大ガス塊の速度差にも大きく依存する。従って、上記のように構成することにより、通過体積を求めることができるから、通過体積に関連する混合体積を求めることができる。

このように、前記ガス塊速度取得手段が備えられている場合、ガス塊速度取得手段は、前記衝突反応前における前記他のガス塊の運動量と前記衝突反応前における前記２つのガス塊の前記一部分のそれぞれの運動量との和が前記衝突反応後における前記他のガス塊の運動量に等しいとの運動量保存則を利用して、前記他のガス塊の速度を更新・取得するように構成されることが好適である。

このように、上記衝突反応の前後に亘って他のガス塊と混合・分離に係わった燃料との間で運動量保存則が成立すると仮定すれば、衝突反応後における他のガス塊の速度を簡易、且つ容易に更新することができる。

また、前記ガス塊速度取得手段が備えられている場合、ガス塊速度取得手段は、前記混合気内に新たに取り込んだ前記筒内ガスの速度を「０」として前記筒内ガスの塊の速度を取得するように構成されることが好適である。

時間の経過とともに筒内ガスが混合気内に逐次取り込まれていくのは、燃焼室内において静止している筒内ガス中において燃料噴霧の領域（従って、混合気の領域）が時間の経過とともに広がっていくからであると考えることができる。即ち、混合気内に新たに取り込まれる筒内ガスの速度は「０」であると考えることができる。上記構成は係る観点に基づくものである。

これにより、混合気内に新たに取り込まれた筒内ガスの運動量が「０」となる。従って、例えば、筒内ガスが混合気内に新たに取り込まれたことで変化する筒内ガスの塊（燃料質量割合が「０」のガス塊）の速度を運動量保存則等を用いて計算する場合、その計算を簡易とすることができる。なお、筒内のガスの流れ（スワール・スキッシュ等）を考慮して前記混合気内に新たに取り込んだ筒内ガスの速度を決定してもよい。

また、前記ガス塊速度取得手段が備えられている場合、ガス塊速度取得手段は、前記噴射された燃料の速度を利用して前記燃料蒸気の塊の速度を取得するように構成されることが好適である。混合気の一部として新たに加えられていく燃料蒸気の速度は、燃料噴射弁から噴射される燃料の速度と等しいと考えることができる。燃料噴射弁から噴射される燃料の速度は、種々の論文等にて紹介されている周知の式により簡易、且つ精度良く求めることができる。上記構成は係る観点に基づくものである。

これにより、混合気の一部として新たに加えられていく燃料蒸気の運動量を簡易、且つ精度良く求めることができる。従って、例えば、燃料蒸気が混合気の一部として新たに加えられたことで変化する燃料蒸気の塊（燃料質量割合が「１」のガス塊）の速度を運動量保存則等を用いて計算する場合、簡易、且つ精度良く燃料蒸気の塊の速度を計算することができる。

また、本発明に係る混合気状態取得装置は、上述した本発明に係る内燃機関の混合気の不均一度取得装置により取得された前記混合気内における燃料濃度の不均一度を利用して前記混合気の状態を取得する混合気状態取得手段を備えている。ここにおいて、混合気の「状態」とは、例えば、混合気の温度、燃料濃度等である。

これによれば、混合気内における不均一度が考慮されて混合気の状態が取得される。従って、不均一度が考慮されずに混合気の状態が取得される場合に比して混合気の状態が精度良く取得され得る。

この場合、前記混合気状態取得手段は、前記混合気の温度を取得する混合気温度取得手段と、前記混合気の温度が所定の着火温度を超えたか否かによって前記混合気内での燃料の反応速度を取得する手法を変更して前記燃料の反応速度（燃焼速度）を取得する燃料反応速度取得手段とを備えることが好適である。

より具体的には、前記燃料反応速度取得手段は、例えば、前記混合気の温度が前記所定の着火温度を超えていない場合には低温酸化反応に基づく予混合燃焼を考慮して前記燃料の反応速度を取得するとともに、前記混合気の温度が前記所定の着火温度を超えた場合には高温酸化反応に基づく予混合燃焼と拡散燃焼とを考慮して前記燃料の反応速度を取得するように構成される。ここにおいて、前記着火温度とは、燃料の高温酸化反応（所謂、熱炎反応）が開始される温度であり、これは、実験・シミュレーション等を通して取得することができる。

混合気の温度が着火温度を超えていない場合、混合気内の燃料の燃焼態様として低温酸化反応（所謂、冷炎反応）に基づく予混合燃焼が支配的であることが知られている。一方、混合気の温度が着火温度を超えた場合、混合気内の燃料の燃焼態様として高温酸化反応に基づく予混合燃焼に加えて拡散燃焼を考慮して燃料の反応速度を計算すると、その計算結果に基づく熱発生率が実測値に精度良く一致することが判明した（詳細は、後述する）。

以上のことから、上記のように構成すれば、混合気の温度が着火温度を超える前後に亘って燃料の反応速度（燃焼速度）を精度良く計算することができ、この結果、燃焼騒音、発生トルク等を精度良く予測することができる。

更には、前記混合気の温度が前記所定の着火温度を超えた場合であって燃料の噴射がなお継続する場合、前記燃料反応速度取得手段は、前記混合気を、前記混合気の温度が前記所定の着火温度を超える前に噴射された燃料に基づく着火前噴射部分と前記混合気の温度が前記所定の着火温度を超えた後に噴射された燃料に基づく着火後噴射部分とに区分して個別に扱うとともに、前記着火前噴射部分についての燃料の反応速度を前記高温酸化反応に基づく予混合燃焼と前記拡散燃焼とを考慮して取得し、前記着火後噴射部分についての燃料の反応速度を前記拡散燃焼のみを考慮して取得するように構成されることが好適である。

混合気の温度が前記所定の着火温度を超えた場合であって燃料の噴射がなお継続する場合については、上記着火前噴射部分では上記と同様に高温酸化反応に基づく予混合燃焼に加えて拡散燃焼が発生すると考える一方で上記着火後噴射部分では拡散燃焼のみが発生すると考えて混合気全体についての燃料の反応速度を計算すると、その計算結果に基づく熱発生率が実測値により一層精度良く一致することが判明した（詳細は、後述する）。

以上のことから、上記のように構成すれば、混合気の温度が前記所定の着火温度を超えた場合であって燃料の噴射がなお継続する場合であっても燃料の反応速度（燃焼速度）を精度良く計算することができ、この結果、燃焼騒音、発生トルク等を精度良く予測することができる。

以下、本発明による内燃機関（ディーゼル機関）の混合気の不均一度取得装置を含んだ混合気状態取得装置の実施形態について図面を参照しつつ説明する。この混合気状態取得装置は、混合気の状態に加えてクランク角度に対する燃料の反応に基づく熱発生率（熱発生速度）を推定する。更には、この装置は、混合気内における燃料濃度の不均一度も推定する。

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

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

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

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

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

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

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

インターフェース６５は、吸気管３２に配置された熱線式エアフローメータ７１、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に設けられた吸気温センサ７２、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気管圧力センサ７３、クランクポジションセンサ７４、アクセル開度センサ７５、スロットル弁３３の下流であって排気還流管５１が接続された部位よりも下流の吸気通路に配設された吸気酸素濃度センサ７６、及び排気マニホールド４１の下流側集合部に設けられた排気酸素濃度センサ７７と接続されていて、これらのセンサからの信号をＣＰＵ６１に供給するようになっている。また、インターフェース６５は、燃料噴射弁２１、燃料噴射用ポンプ２２、スロットル弁アクチュエータ３３ａ、及びＥＧＲ制御弁５２と接続されていて、ＣＰＵ６１の指示に応じてこれらに駆動信号を送出するようになっている。

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

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

（混合気状態の推定方法の概要）
次に、上記のように構成された混合気の不均一度取得装置を含んだ混合気状態取得装置（以下、「本装置」と云う。）による混合気状態の推定方法について説明する。

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

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

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

そして、ピストンが圧縮上死点近傍に達すると（具体的には、後述する燃料噴射（開始）時期（クランク角度）CAinjが到来すると）、本装置は、前記燃料噴射量Qfinに応じた噴射期間TAUだけ燃料噴射弁２１を開弁する。この結果、燃料噴射弁２１から液体の燃料が上記燃料噴射時期CAinjから噴射期間TAUだけ連続して燃焼室内にて（従って、筒内ガスに向けて）直接噴射される。この噴射された燃料は、圧縮により高温になっている筒内ガスから受ける熱により直ちに燃料蒸気になるとともに、時間の経過に伴って同筒内ガスを取り込みながら混合気（＝燃料蒸気＋筒内ガス）となって燃焼室内において円錐状に拡散していく。

図３は、ＩＶＣ後、燃料噴射開始前において上記筒内ガスのみが存在する燃焼室内の様子を示した図であり、図４は、燃料噴射開始後の或る時点において燃焼室内にて燃料蒸気（従って、混合気）が広がっている様子を模式的に示した図である。

図３、及び図４に示すように、本装置では、燃焼室内を、筒内ガスが占める領域と混合気が占める領域とに分けて取り扱い、これらをそれぞれ「ゾーン１」、「ゾーン２」と称呼する。従って、燃焼室は、燃料噴射開始前はゾーン１のみから構成され、燃料噴射開始後はゾーン１、及びゾーン２から構成される。

なお、後述するように、本装置では、混合気着火後もなお燃料噴射が継続する場合、混合気を、着火時点以前において噴射された燃料に基づく部分（前記着火前噴射部分）と着火時点以降において噴射された燃料に基づく部分（前記着火後噴射部分）とに更に分けて取り扱う。この場合、着火前噴射部分が占める領域をそのまま「ゾーン２」と称呼するとともに着火後噴射部分が占める領域を「ゾーン３」と称呼する。

以下、ゾーン１、ゾーン２、及びゾーン３に係わる変数、記号等にはそれぞれその末尾に「z1」「z2」「z3」を付すことにする。また、説明の便宜上、先ずは、混合気着火後において燃料噴射が継続する場合であっても、前記着火前噴射部分と前記着火後噴射部分とを分けることなく上記混合気が占める領域はゾーン２のみからなるものとして説明を続ける。ゾーン３については後に詳述する。

図４に示すように、噴射された燃料（従って、燃料蒸気）は、噴霧角θをもって円錐状に広がりながら筒内ガスを順次取り込んでいく（即ち、筒内ガスと順次混ざり合っていく）ものとする。そして、燃料噴射開始時期CAinjからの任意の経過時間（以下、「噴射後経過時間tz2」と称呼する。）において、質量Qz2の燃料（従って、燃料蒸気）が既に噴射されていてその燃料蒸気が体積Vz2の円錐状に広がっているものとする。

更には、この結果、噴射後経過時間tz2において、上記燃料蒸気が質量Gz2の筒内ガス（以下、「混合気形成筒内ガス」と云うこともある。）と混ざり合って質量Mz2＝(Qz2＋Gz2）、体積Vz2の混合気となっているものとする。換言すれば、噴射後経過時間tz2において、ゾーン２の体積がVz2に、ゾーン２内のガスの総質量がMz2になっているものとする。

本装置は、このように形成されていく混合気の状態（即ち、ゾーン２内のガス状態）を推定する。具体的には、燃料噴射開始前ではゾーン１内のガス状態が推定され、燃料噴射開始後ではゾーン１，２内のガス状態がそれぞれ個別に推定される。ガス状態としては、上記熱発生率の推定に必要となる各種物理量（ガス温度、ガス体積、ガス質量等）が推定される。以下、ゾーン１内のガス状態とゾーン２内のガス状態のそれぞれの推定方法について順に説明する。

<ゾーン１内のガス状態>
先ず、ＩＶＣ以降におけるゾーン１内のガス質量Gz1について説明する。ＩＶＣ後においてゾーン１は、燃料噴射の前後にかかわらず筒内ガスのみから構成される。従って、ゾーン１内のガス質量Gz1（＝ゾーン１内の筒内ガスの質量）は下記(1)式にて表すことができる。

(1)式において、Gz2は上述したようにゾーン２内に取り込まれている筒内ガス（即ち、上記混合気形成筒内ガス。以下、「ゾーン２内筒内ガス」とも称呼する。）の質量である。ゾーン２内筒内ガス質量Gz2は燃料噴射開始前では「０」に維持される。燃料噴射開始後におけるゾーン２内筒内ガス質量Gz2の取得方法については後述する。

(1)式において、GcはＩＶＣ以降燃焼室内に密閉されている筒内ガスの全質量（全筒内ガス質量）である。全筒内ガス質量はGcは、ＩＶＣにおける気体の状態方程式に基づく下記(2)式に従って取得され得る。

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

(2)式において、TgivcはＩＶＣにおける筒内ガス温度である。ＩＶＣは圧縮下死点近傍であるから、ＩＶＣにおいて筒内ガス温度は吸気温度Tbと略等しいと考えられる。従って、ＩＶＣにおいて吸気温センサ７２により検出される吸気温度TbがTgivcとして使用され得る。(2)式において、Rは筒内ガスのガス定数である。ここでのガス定数Rは、実際には一般ガス定数を筒内ガスの平均分子量で除した値であるが、本例では一定として扱う。

以上より、(2)式の右辺の各因数の値が全て取得され得るから(2)式に従って全筒内ガス質量Gcを取得することができる。従って、(1)式の右辺の各項の値が全て取得され得るから(1)式に従ってＩＶＣ以降におけるゾーン１内のガス質量Gz1を逐次取得することができる。なお、燃料噴射開始前では、上述したようにゾーン２内筒内ガス質量Gz2＝０であるから(1)式よりゾーン１内のガス質量Gz1は全筒内ガス質量Gcと等しくなる（図３、図４を参照）。

次に、ＩＶＣ以降におけるゾーン１の体積Vz1について説明する。ゾーン１の体積Vz1は下記(3)式にて表すことができる。

(3)式において、Vz2は上述したようにゾーン２の体積である。ゾーン２の体積Vz2は燃料噴射開始前では「０」に維持される。燃料噴射開始後におけるゾーン２の体積Vz2の取得方法については後述する。以上より、(3)式の右辺の各項の値が全て取得され得るから(3)式に従ってＩＶＣ以降におけるゾーン１の体積Vz1を逐次取得することができる。なお、燃料噴射開始前では、上述したようにゾーン２の体積Vz2＝０であるから(3)式よりゾーン１の体積Vz1は筒内容積Vc(CA)と等しくなる（図３、図４を参照）。

次に、ＩＶＣ以降におけるゾーン１内のガス温度Tz1について説明する。ゾーン１内のガス温度Tz1は下記(4)式にて表すことができる。

(4)式において、Cpz1はゾーン１内のガス（従って、筒内ガス）の定圧比熱である。以下、ゾーン１内のガスの定圧比熱Cpz1の求め方について説明する。一般に、或るガスについて、比熱比をκ、定圧比熱をCp、ガス定数をRとすると、下記(5)式が成立する。従って、ガスの定圧比熱はガスの比熱比に依存する。

ここで、ガスの比熱比は、そのガスを構成する成分の組成割合に依存する。機関の燃焼室内のガスは、燃料噴射開始前では筒内ガスのみから構成され、燃料噴射開始後では筒内ガスと燃料蒸気とから構成される。筒内ガスを構成する成分は主として、酸素Ｏ_２、二酸化炭素ＣＯ_２、窒素Ｎ_２、水Ｈ_２Ｏを含み、これらの組成割合は吸気中の酸素濃度、前回の排気行程中における排気中の酸素濃度、及び現時点でのガスの温度に大きく依存する。

従って、一般に、燃焼室内のガスの比熱比（従って、燃焼室内のガスの定圧比熱）は、吸気中の酸素濃度、前回の排気行程中における排気中の酸素濃度、現時点でのガスの温度、及びガス中の燃料濃度（今回噴射された燃料の濃度）に大きく依存する。以上より、一般に、燃焼室内のガスの定圧比熱Cpcは、吸気酸素濃度[O2]in，前回の排気行程中における所定時期での排気酸素濃度[O2]exb，現時点でのガス温度Tg，ガス中の燃料濃度[Fuel]を引数とする関数funcCpを用いて下記(6)式に従って表すことができる。

ここで、ゾーン１内には燃料蒸気が存在しないから、ゾーン１内での燃料濃度は「０」である。従って、ＩＶＣ以降におけるゾーン１内のガスの定圧比熱Cpz1は、下記(7)式に従って逐次求めることができる。(7)式において、吸気酸素濃度[O2]inとしては吸気酸素濃度センサ７６から得られる今回の吸気についての吸気酸素濃度RO2inが使用され得、排気酸素濃度[O2]exbとしては排気酸素濃度センサ７７から得られる前回の排気行程中における所定時期での排気酸素濃度RO2exが使用され得る。以上、ゾーン１内のガスの定圧比熱Cpz1の求め方について説明した。

再び、(4)式を参照すると、dWz1は、クランク角度CAが微小クランク角度dCAだけ進行する間においてゾーン１内のガス（従って、筒内ガス）が受ける微小エネルギーである。本例では、下記(8)式に示すように、微小エネルギーdWz1は、クランク角度CAが微小クランク角度dCAだけ進行する間においてピストンがゾーン１内のガスに対して行う仕事（微小ピストン仕事dWpistonz1）と等しいものとする。以下、微小ピストン仕事dWpistonz1の求め方について説明する。

ゾーン１の体積Vz1の筒内容積Vc(CA)に対する割合は「Vz1/Vc(CA)」と表すことができる。従って、クランク角度CAが微小クランク角度dCAだけ進行する間において、筒内容積がdV（＞０）だけ減少したとするとゾーン１の体積Vz1は「dV・Vz1/Vc(CA)」（＞０）だけ減少する。従って、微小ピストン仕事dWpistonz1は、現時点での燃焼室内圧力Pgを用いて下記(9)式に従って表すことができる。

(9)式において、燃焼室内圧力Pgは、本例では燃焼室内で均一とする。燃焼室がゾーン１のみから構成される燃料噴射開始前では、燃焼室内圧力Pgは、現時点でのゾーン１内のガス質量Gz1、現時点でのゾーン１内のガス温度Tz1、現時点でのゾーン１の体積Vz1と、現時点でのゾーン１内のガスについての状態方程式とに基づいて「Pg＝Gz1・R・Tz1/Vz1」なる関係から求めることができる。ゾーン２が発生する燃料噴射開始後における燃焼室内圧力Pgについては後にフローチャートを参照しながら説明する。

また、筒内容積Vcをクランク角度CAで微分した値(dVc/dCA)は機関１０の設計諸元に基づいてクランク角度CAの関数(dVc/dCA)(CA)として取得することができるから、dVは下記(10)式に従って表すことができる。従って、(9)式、及び(10)式により、微小ピストン仕事ΔWpistonz1は下記(11)式に従って逐次求めることができる。従って、(8)式に従ってＩＶＣ以降における微小エネルギーdWz1も逐次求めることができる。

以上より、(4)式の右辺の各変数の値が全て取得され得るから(4)式に従ってＩＶＣ以降におけるゾーン１内のガス温度Tz1を逐次取得することができる。以上のように取得されるゾーン１内のガスに関する各物理量は、上述したように上記熱発生率の推定に必要な値である。以上、ゾーン１内のガス状態の推定方法について説明した。

<ゾーン２内のガス状態>
次に、燃料噴射開始後（即ち、燃料噴射時期CAinj以降）において発生するゾーン２内のガス状態の推定方法について説明する。先ず、燃料噴射開始後における（即ち、噴射後経過時間tz2における）ゾーン２（従って、混合気が占める領域）の体積Vz2について説明する。

図４に示すように、燃料噴射開始後における混合気は噴霧角θをもって円錐状に広がっていく。図５に示すように、この円錐の体積を噴霧体積VOLtotalと称呼するものとすると、ゾーン２の体積Vz2は噴霧体積VOLtotalに等しい。以下、噴霧体積VOLtotalの求め方について説明する。

燃料噴射弁２１から噴射された燃料（燃料蒸気）の先頭部の噴射後経過時間tz2における到達距離Ｓは、例えば、下記(12)式に従って求めることができる。(12)式において、ΔP0は有効噴射圧力であって燃料噴射圧力（本例では、基本燃料噴射圧力Pcrbase）から燃料噴射開始時（噴射後経過時間tz2＝０）における燃焼室内圧力Pgを減じた値である。ρfは液体燃料の密度であり、ρaは空気の密度である。d0は燃料噴射弁２１の噴孔径である。(12)式は、自動車技術会論文集No.21,1980 「ディーゼル噴霧の到達距離と噴霧角」、広安博之、新井雅隆に紹介されている。

また、噴霧角θは、燃料噴射開始時点（噴射後経過時間tz2＝０）における筒内ガスの密度（従って、噴射後経過時間tz2＝０におけるゾーン１内のガス密度ρ0z1（＝Gz1/Vz1））、及び上記有効噴射圧力ΔP0に応じて変化すると考えられるから、上記密度ρ0z1及び有効噴射圧力ΔP0と、噴霧角θとの関係を予め規定したテーブルMapθに基づいて取得することができる。

このように、噴霧角θ、及び噴射後経過時間tz2における燃料蒸気先頭部の到達距離Ｓが取得できれば、噴射後経過時間tz2における噴霧体積VOLtotal（図５を参照）は、下記(13)式に従って求めることができる。従って、噴射後経過時間tz2におけるゾーン２の体積Vz2（＝VOLtotal）も求めることができる。

次に、燃料噴射開始から燃料噴射後経過時間tz2までにゾーン２内に取り込まれた筒内ガス（即ち、ゾーン２内筒内ガス）の質量Gz2について説明する。本例では、クランク角度CAが微小クランク角度dCAだけ進行する間においてゾーン２内に取り込まれる筒内ガスの体積は、クランク角度CAが微小クランク角度dCAだけ進行する間におけるゾーン２の体積Vz2の増加量dVz2に等しいと考える。加えて、クランク角度CAが微小クランク角度dCAだけ進行する間においてゾーン２内に取り込まれる筒内ガスの質量dGz2は、下記(14)式に従って表すことができる。

(14)式において、ρz1はゾーン１内のガス（即ち、筒内ガスそのもの）の密度（＝Gz1/Vz1）である。従って、燃料噴射後におけるゾーン２内筒内ガス質量Gz2は、下記(15)式に従って求めることができる。

次に、燃料噴射開始後におけるゾーン２内のガス質量Mz2について説明する。上述したように、噴射後経過時間tz2におけるゾーン２内のガス質量Mz2は、噴射後経過時間tz2までにおいて既に噴射された燃料の質量Qz2と、(15)式にて求められる噴射後経過時間tz2におけるゾーン２内筒内ガス質量Gz2の和となる。ここで、上記質量Qz2は、上記基本燃料噴射圧力Pcrbase、噴射期間TAU等に基づいて取得することができる。従って、噴射後経過時間tz2におけるゾーン２内のガス質量Mz2（＝Gz2＋Qz2）も求めることができる。

次に、燃料噴射開始後におけるゾーン２内のガス温度Tz2（即ち、混合気温度）について説明する。ゾーン２内のガス温度Tz2は下記(16)式にて表すことができる。(16)式において、T0z2は、燃料噴射開始時点（即ち、噴射後経過時間tz2＝０）におけるゾーン１内のガス温度Tz1であり、上記(4)式から求めることができる。

(16)式において、Cpz2はゾーン２内のガス（従って、混合気）の定圧比熱であり、上記(6)式に基づいて下記(17)式に従って逐次求めることができる。(17)式において、[Fuel]z2はゾーン２内の燃料濃度（噴射された燃料の質量濃度。ＥＧＲ等による筒内ガス中の残留燃料分は含まれない。）である。ゾーン２内の燃料濃度[Fuel]z2の求め方については後述する。

(16)式において、dWz2は、クランク角度CAが微小クランク角度dCAだけ進行する間においてゾーン２内のガス（従って、混合気）が受ける微小エネルギーである。本例では、微小エネルギーdWz2は下記(18)式にて表されるものとする。(18)式において、dWpistonz2は、クランク角度CAが微小クランク角度dCAだけ進行する間においてピストンがゾーン２内のガスに対して行う仕事（微小ピストン仕事）であり、Qlatentz2は、クランク角度CAが微小クランク角度dCAだけ進行する間において噴射された液体燃料がゾーン２内にて燃料蒸気に変化する際の潜熱Qlatentz2であり、Hrz2はゾーン２内における単位クランク角度あたりの燃料の反応に基づく熱発生率である。この熱発生率Hrz2の求め方については後述する。

(18)式において、微小ピストン仕事dWpistonz2は、上記(11)式にて求められる微小ピストン仕事dWpistonz1と同じ考え方に基づいて、下記(19)式に従って逐次求めることができる。

(18)式において、潜熱Qlatentz2は、クランク角度CAが微小クランク角度dCAだけ進行する間において噴射された液体燃料の質量qに応じた値になると考えられるから、qを引数とする関数funcQlatent(q)を用いて逐次求めることができる。なお、係る質量qは、上記基本燃料噴射圧力Pcrbase、噴射期間TAU等に基づいて取得することができる。以上より、(18)式の右辺の各項の値が全て取得され得るから、(18)式に従って燃料噴射開始後における微小エネルギーdWz2も逐次取得することができる。

以上より、(16)式の右辺の各変数の値が全て取得され得るから(16)式に従って燃料噴射開始後におけるゾーン２内のガス温度Tz2を逐次取得することができる。以上のように取得されるゾーン２内のガスに関する各物理量は、上述したように上記熱発生率の推定に必要な値である。以上、ゾーン２内のガス状態（即ち、混合気の状態）の推定方法について説明した。以上が、混合気状態の推定方法の概要である。

(ゾーン２内における燃料の反応に基づく熱発生率)
次に、上記ゾーン２内（即ち、混合気内）における単位クランク角度あたりの燃料（今回噴射された燃料。ＥＧＲ等により筒内ガスに含まれる燃料分を含まない。）の反応に基づく熱発生率Hrz2(J/deg)について説明する。本例では、上記ゾーン２内のガス温度Tz2が所定の着火温度Tig（一定）を超えた時点で混合気が着火すると仮定する。以下、ゾーン２のガス温度Tz2が着火温度Tig以下の段階を「着火前」、着火温度Tigを超えた段階を「着火後」と称呼することもある。

上記熱発生率Hrz2を求めるためには、ゾーン２における燃料の反応速度（燃焼速度）qrz2(g/sec)を求める必要がある。着火前と着火後では燃料の燃焼態様が異なるから、本例ではゾーン２における燃焼速度qrz2を着火前と着火後で分けて求める。

着火前では、ゾーン２内の燃料の燃焼態様として低温酸化反応（即ち、冷炎反応）に基づく予混合燃焼が支配的である。従って、着火前でのゾーン２における燃焼速度qrz2(g/sec)は、低温酸化反応に基づく予混合燃焼により決定されると考えられる。

低温酸化反応に基づく予混合燃焼は、所謂シェル・モデル（Shell Model M.P.Halstead,1977,Reitz,SAE950278）により精度良く模擬できることが広く知られている。以上より、本例では、着火前でのゾーン２における燃焼速度qrz2(g/sec)は、ゾーン２内における酸素濃度[O2]z2、ゾーン２内における燃料濃度[Fuel]z2、ゾーン２内のガス温度Tz2を引数とするシェル・モデルに対応する関数funcshell([O2]z2,[Fuel]z2,Tz2)を用いて求められる。[O2]z2，[Fuel]z2の求め方については後述する。

一方、着火後では、ゾーン２内の燃料の燃焼態様として高温酸化反応に基づく予混合燃焼と拡散燃焼とが発生し得ると考えられる。これは、以下の理由に基づく。即ち、ゾーン２内において着火時点にて噴射後比較的長い時間が経過している燃料に基づく領域（即ち、ゾーン２内における燃料噴射弁２１の噴孔から比較的遠い領域）では、燃料濃度が比較的均一になっていると考えられる。従って、この領域では、同領域内の燃料がほぼ同時に着火（燃焼）する予混合燃焼（或いは、予混合燃焼に近い燃焼）が発生すると考えられる。

一方、ゾーン２内において着火時点にて噴射後比較的短い時間しか経過していない燃料に基づく領域（即ち、ゾーン２内における燃料噴射弁２１の噴孔から比較的近い領域）では、燃料濃度が比較的不均一になっていると考えられる。従って、この領域では、噴射された燃料が、筒内ガスと混ざり合いながらその燃料濃度が可燃範囲内に入った段階で順次着火（燃焼）していく拡散燃焼が発生すると考えられる。

従って、着火後でのゾーン２における燃焼速度qrz2(g/sec)は、高温酸化反応に基づく予混合燃焼と拡散燃焼とにより決定されると考えられる。着火後でのゾーン２における燃焼速度qrz2(g/sec)は下記(20)式に従って求められる。

(20)式において、τcz2は着火後におけるゾーン２内の燃料の燃焼に係わる特性時間である。(20)式から理解できるように、特性時間τcz2が長いほどゾーン２における燃焼速度qrz2は小さくなる。即ち、特性時間τcz2は、ゾーン２における燃焼速度qrz2の遅さを指標する値である。特性時間τcz2は下記(21)式に従って求められる。

(21)式において、τaz2はゾーン２内における高温酸化反応に基づく予混合燃焼に係わる特性時間（層流特性時間）であり、本例では、所謂アレニウス（Arrehenius）の式に基づく下記(22)式で取得される。(22)式において、Eは活性化エネルギー（ここでは、定数）であり、C1は定数である。(22)式から理解できるように、ゾーン２内における酸素濃度[O2]z2と燃料濃度[Fuel]z2の積が大きいほど、或いはゾーン２内のガス温度Tz2が高いほど、層流特性時間τaz2は短くなる（従って、特性時間τcz2も短くなる）。

(21)式において、τmz2はゾーン２内における拡散燃焼に係わる特性時間（乱流特性時間）であり、本例では、所謂マグネッセン（Magnussenn）の式に基づく下記(23)式で取得される。(23)式において、C2は定数である。kは乱流エネルギー強度であり、εは乱流エネルギー散逸率である。値(k/ε)は、上記有効噴射圧力ΔP0と液体燃料の密度ρfとに大きく依存するから、上記有効噴射圧力ΔP0と、液体燃料密度ρfと、ΔP0，ρfを引数とする関数funck/ε(ΔP0,ρf)とに基づいて取得することができる。funcτmは、ゾーン２内における酸素濃度[O2]z2と燃料濃度[Fuel]z2とに基づく乱流特性時間τmの補正係数を求める関数である。

(21)式において、C3は拡散燃焼の特性時間τcz2への影響度合いを表す係数であり、係数C3は、燃料噴射量Qfinに対する燃焼した（消費された）燃料量（消費燃料量sumqrz2）の割合（＝sumqrz2/Qfin）が増大するほど大きい値に設定される。これにより、着火後の燃料の燃焼の進行に応じて（即ち、値(sumqrz2/Qfin)の増大に応じて）、乱流特性時間τmz2（従って、拡散燃焼）の特性時間τcz2への影響度合いが大きくなる。これは、着火後の燃焼は、当初は高温酸化反応に基づく予混合燃焼が支配的である一方、燃焼の進行に応じて拡散燃焼が徐々に支配的になっていく傾向があることに基づく。なお、消費燃料量sumqrz2の求め方については後述する。

以上により、着火前と着火後で燃料の燃焼態様が異なることを考慮してゾーン２における燃焼速度qrz2(g/sec)が着火前後を問わず精度良く逐次求められ得る。この燃焼速度qrz2を用いることで燃料噴射開始後における上記ゾーン２の熱発生率Hrz2(J/deg)は、下記(24)式に従って逐次求めることができる。

(24)式において、Hfは燃料の反応に基づく燃料の単位質量あたりの発熱量(J/g)である（一定値）。dt/dCAは単位時間あたりに進行するクランク角度CAの値の逆数であり、エンジン回転速度NEにより決定される値である。

<燃料濃度[Fuel]z2の取得>
次に、燃料噴射開始後におけるゾーン２内における燃料濃度（質量濃度）[Fuel]z2の取得方法について説明する。噴射後経過時間tz2におけるゾーン２内の燃料濃度[Fuel]z2は、噴射後経過時間tz2における上記ゾーン２内のガス質量Mz2に対する、「噴射後経過時間tz2においてゾーン２内に存在する燃料（蒸気）の質量」の割合である。

ここで、「噴射後経過時間tz2においてゾーン２内に存在する燃料の質量」は、噴射開始から噴射後経過時間tz2までにおいて既に噴射された燃料の質量（即ち、上記質量Qz2）から、噴射開始から噴射後経過時間tz2までにおいて化学反応（即ち、燃焼）で消費された燃料分（即ち、上記消費燃料量sumqrz2）を減じた値である。従って、ゾーン２内に亘って燃料濃度は均一であるものと仮定すれば、噴射後経過時間tz2におけるゾーン２内の燃料濃度[Fuel]z2は、下記(25)式で表すことができる。

(25)式において、消費燃料量sumqrz2は、上述のように求められるゾーン２における燃焼速度qrz2(g/sec)の積算であるから下記(26)式に従って求めることができる。従って、(25)式の右辺の各変数の値が取得され得る。この結果、(25)式に従って、ゾーン２内に亘って燃料濃度は均一であるものと仮定した場合における燃料噴射開始後におけるゾーン２内の燃料濃度[Fuel]z2を逐次求めることができる。

<酸素濃度[O2]z2の取得>
次に、燃料噴射開始後におけるゾーン２内における酸素濃度（質量濃度）[O2]z2の取得方法について説明する。噴射後経過時間tz2におけるゾーン２内の酸素濃度[O2]z2は、噴射後経過時間tz2における上記ゾーン２内のガス質量Mz2に対する、「噴射後経過時間tz2においてゾーン２内に存在する酸素の質量」の割合である。

「噴射後経過時間tz2においてゾーン２内に存在する酸素の質量」は、「噴射後経過時間tz2における上記ゾーン２内筒内ガス中の酸素の質量」から、「噴射開始から噴射後経過時間tz2までにおいて化学反応（即ち、燃焼）で消費された酸素分」を減じた値である。

ここで、「噴射後経過時間tz2における上記ゾーン２内筒内ガス中の酸素の質量」は、上記ゾーン２内筒内ガスの質量Gz2に吸気酸素濃度[O2]inを乗じた値に等しい。また、「噴射開始から噴射後経過時間tz2までにおいて化学反応で消費された酸素分」は、上記消費燃料量sumqrz2に、化学反応時の質量割合Ho（酸素質量／燃料質量、一定）を乗じた値に等しい。従って、燃料噴射開始後におけるゾーン２内の酸素濃度[O2]z2は、下記(27)式に従って逐次求めることができる。

以上が、燃料噴射開始後における上記ゾーン２の熱発生率Hrz2(J/deg)の取得方法の概要である。ただし、上述したように、これにより取得される上記ゾーン２の熱発生率Hrz2(J/deg)は、「混合気着火後において燃料噴射が継続する場合であっても前記着火前噴射部分と前記着火後噴射部分とを分けることなく混合気が占める領域はゾーン２のみからなるとの仮定（以下、「仮定１」と称呼する。）」、並びに、「ゾーン２における燃焼速度qrz2(g/sec)（従って、ゾーン２の熱発生率Hrz2(J/deg)）を求める際に使用されるゾーン２内の燃料濃度［Fuel]z2はゾーン２内で均一であるとの仮定（以下、「仮定２」と称呼する。）」のもとで取得される値である。

図６は、混合気着火後においても燃料噴射がなお継続する場合におけるクランク角度CAに対する熱発生率(J/deg)について、実験結果（破線を参照）と上述した取得方法により取得される計算結果（即ち、上記ゾーン２の熱発生率Hrz2、実線を参照）との比較結果の一例を示したグラフである。

図６から理解できるように、この計算結果は、実験結果と大略的には一致しているものの、以下の２つの観点で実験結果から大きく乖離している。先ず１つめの観点（Ａを参照）は、実験結果では１回目の比較的大きいピークの後で２回目の比較的小さいピークが発生している一方で計算結果では２回目のピークが発生しない点である（以下、「問題点１」と称呼する。）。２つめの観点（Ｂを参照）は、計算結果でのピーク（１回目のピーク）が実験結果での１回目のピークに比して過大となっている点である（以下、「問題点２」と称呼する。）。

本発明者は、上記問題点１は上記仮定１に起因し、上記問題点２は上記仮定２に起因することを見出した。即ち、上記問題点１は、混合気着火後において燃料噴射が継続する場合は前記着火前噴射部分と前記着火後噴射部分とを分けて取り扱うことで解消し、上記問題点２は、混合気内における燃料濃度の不均一性を考慮することで解消することを見出した。以下、本装置による上記問題点１，２の解消方法について具体的に説明する。

（問題点１の解消方法）
上述したように、本装置では、混合気着火後もなお燃料噴射が継続する場合（以下、「着火後噴射継続の場合」と称呼することもある。）、図７に示すように、混合気が、着火時点以前において噴射された燃料に基づく部分（前記着火前噴射部分）と着火時点以降において噴射された燃料に基づく部分（前記着火後噴射部分）とに分けて取り扱われる。着火前噴射部分が占める領域はそのまま「ゾーン２」と称呼され、着火後噴射部分が占める領域は「ゾーン３」と称呼される。

即ち、着火後噴射継続の場合、着火時点以降、ゾーン３が発生する。燃料の噴射は、着火時点以前では上述のようにゾーン２に対して行われる一方、着火時点以降ではゾーン２に代えてゾーン３に対して行われるものと仮定する。換言すれば、着火後噴射継続の場合における着火後では、それまで継続していたゾーン２内への燃料（蒸気）の供給が停止される。

そして、着火時点以降、ゾーン３内（即ち、混合気内）における単位クランク角度あたりの燃料（今回噴射された燃料。ＥＧＲ等により筒内ガスに含まれる燃料分を含まない。）の反応に基づく熱発生率Hrz3(J/deg)が、後述するように、上述したゾーン２の熱発生率Hrz2(J/deg)とは異なる計算手法により同ゾーン２の熱発生率Hrz2と並行して別個独立して逐次計算されていく。

上述のように、ゾーン３は着火時点以降から発生する。着火時点からの任意の経過時間（以下、「着火後経過時間tz3」と称呼する。）におけるゾーン３内の各物理量は、上記熱発生率Hrz3を除いて、上述した噴射後経過時間tz2におけるゾーン２内の各物理量の取得方法と同様の方法で取得され得る。ゾーン３内の各物理量の取得については後にフローチャートを参照しながら説明する。

なお、ゾーン３が発生した時点以降（即ち、着火後噴射継続の場合における着火時点以降）、ゾーン３の体積Vz3、及び、着火時点から着火後経過時間tz3までにゾーン３内に取り込まれた筒内ガス（以下、「ゾーン３内筒内ガス」と称呼する。）の質量Gz3が発生する。このことにより、ゾーン３が発生した時点以降、ゾーン２の体積Vz2は上述した「Vz2＝VOLtotal」に代えて下記(28)式に従って計算され、ゾーン１内のガス質量Gz1及びゾーン１の体積Vz1は、上記(1)式及び(3)式に代えて、下記(29)式及び(30)式に従ってそれぞれ計算される（図４と図７を参照）。

<ゾーン３内における燃料の反応に基づく熱発生率>
次に、上記ゾーン３（即ち、混合気）内の熱発生率Hrz3(J/deg)について説明する。上記熱発生率Hrz2と同様、熱発生率Hrz3を求めるためには、ゾーン３における燃料の反応速度（燃焼速度）qrz3(g/sec)を求める必要がある。

先に、ゾーン２内での着火後では、ゾーン２内の燃料の燃焼態様として高温酸化反応に基づく予混合燃焼と拡散燃焼とが発生し得ると述べた。しかしながら、着火後噴射継続の場合、着火時点以降において噴射された燃料に基づく部分（着火後噴射部分）が占める領域内（即ち、ゾーン３内）に限っては、拡散燃焼が支配的であると考えられる。これは以下の理由に基づく。

即ち、ゾーン３は着火時点以降に発生するから、ゾーン３内のガス温度Tz3はゾーン３発生当初から上記着火温度Tigを超えている。従って、ゾーン３内に噴射された燃料（蒸気）は、筒内ガスと混ざり合いながらその燃料濃度が可燃範囲内に入った段階で（噴射から極短時間後に）直ちに順次着火していく。

換言すれば、（着火後での）ゾーン３では拡散燃焼のみが発生し、ゾーン３における燃焼速度qrz3(g/sec)は、燃料が筒内ガスと混ざり合っていく速度のみに依存すると考えられる。従って、（着火後での）ゾーン３における燃焼速度qrz3(g/sec)は、上記(20)式に対応する下記(31)式に従って求められる。

(31)式において、τmz3はゾーン３内における拡散燃焼に係わる特性時間（乱流特性時間）である。(31)式から理解できるように、乱流特性時間τmz3が長いほどゾーン３における燃焼速度qrz3は小さくなる。乱流特性時間τmz3は、上記(23)式と同様、所謂マグネッセン（Magnussenn）の式に基づく下記(32)式で取得される。

このように、本装置では、着火後噴射継続の場合、ゾーン３が発生する着火後において、ゾーン２とゾーン３とを分けて取り扱いそれぞれの内部での燃焼速度（反応速度）qrz2,qrz3(g/sec)（従って、それぞれの熱発生率Hrz2,Hrz3(J/deg)）を別個独立に並行して順次求めていく。その際、計算の基礎となる燃焼形態（従って、計算に使用される式）をゾーン２とゾーン３とで異ならせる。

そして、着火後においては、クランク角度CAに対する熱発生率(J/deg)を(Hrz2＋Hrz3)として計算していく。これにより、本発明者は、上記問題点１が解消されることを見出した。この効果については後に図１４を参照しながら説明する。以上が、本装置による上記問題点１の解消方法である。

（問題点２の解消方法）
上記問題点２を解消するため、本装置では、混合気内における燃料濃度の不均一性が考慮される。ここで、混合気内における燃料濃度の不均一性を考慮するために、上述した燃焼速度qrz2,qrz3(g/sec)（従って、熱発生率Hrz2,Hrz3(J/deg)）の計算に使用される上記燃料濃度［Fuel]z2,[Fuel]z3そのものの不均一度を求める手法は、従来より種々の論文等にて紹介されてきている。

しかしながら、係る手法はいずれも、例えば、確率密度関数等を用いた複雑な微分方程式等を解く必要があり膨大な計算負荷を伴うものばかりだった。そこで、本発明では、上述した燃焼速度qrz2,qrz3(g/sec)（従って、熱発生率Hrz2,Hrz3(J/deg)）の計算に使用される燃料濃度［Fuel]z2,[Fuel]z3そのものとしては上述と同様に混合気内で均一であると仮定した場合の値を使用する一方で、複雑な微分方程式等を解く必要のない以下に示す比較的簡易な計算により混合気内における燃料濃度の不均一度（後述する標準偏差σ）を求める。

そして、この不均一度を上述した燃焼速度qrz2,qrz3(g/sec)（従って、熱発生率Hrz2,Hrz3(J/deg)）の計算に反映させることで混合気内における燃料濃度の不均一性が考慮される。以下、本発明による混合気内における燃料濃度の不均一度（標準偏差σ）を求める手法について説明する。

（本発明による混合気の不均一度の取得方法）
本発明では、混合気内における燃料濃度の不均一度（標準偏差σ）を求めるため、混合気は、燃料質量分率（燃料質量割合、燃料質量濃度）が異なる球形のガス塊の集合体であると仮定される。ガス塊の「燃料質量分率」とは、同ガス塊の全質量に対する、同ガス塊に含まれる燃料（蒸気）の質量の割合を意味する。ここにいう「燃料」とは、燃料噴射により新たに燃焼室内に供給された燃料を意味していて、筒内ガスに予め含まれる燃料（例えばＥＧＲ等により筒内ガスに含まれる未燃燃料）を含まないものとする。

また、混合気を構成する上記球形のガス塊の集合体には、筒内ガスそのもの（燃料蒸気を含んでいない）の塊と、燃料蒸気そのもの（筒内ガスを含んでいない）の塊とが含まれるものとする。筒内ガスの塊の燃料質量分率＝０であり、燃料蒸気の塊の燃料質量分率＝１である。

そして、混合気が形成されていく燃料噴射開始時点以降において、微小時間Δtの経過毎に（微小クランク角度ΔCAの進行毎に）、その時点で存在している（即ち、質量が「０」でない）燃料質量分率が異なる複数のガス塊のうちの任意の２つの組み合わせの全てについて以下に説明する「衝突反応」が混合気内で繰り返し行われると仮定される。これにより、混合気内のガス塊の質量分布が取得・更新されていくことで混合気の不均一度が取得・更新されていく。

ここで、「衝突反応」とは、以下の１〜３の一連の反応をいう。
１．燃料質量分率が異なる２つのガス塊が衝突する。
２．衝突した２つのガス塊のそれぞれの一部分が混合する。
３．混合した各一部分が対応するガス塊からそれぞれ離脱して前記２つのガス塊の燃料質量分率とは異なる燃料質量分率を有する他のガス塊の一部又は全部となる。

以下、衝突の対象となる２つのガス塊を「衝突対象ガス塊」と称呼し、「衝突反応」により質量が増大するガス塊（即ち、前記他のガス塊）を「混合形成ガス塊」と称呼することもある。

以下、図８、図９を参照しながら「衝突反応」を利用して混合気内のガス塊の質量分布を取得・更新していく具体的な方法について説明する。なお、以降、図中に示す円の直径の大きさ、及び棒グラフの長さは、その円、及び棒グラフに対応するガス塊の量（質量、体積等）の大きさを表しているものとする。

図８は、燃料噴射開始から微小時間Δtが経過した時点（噴射後経過時間t＝Δt）で発生する「衝突反応」の前後の様子を模式的に示している。図８の左側の図に示すように、噴射後経過時間t＝Δt（衝突前）では、混合気内（上述したゾーン２内に相当する。）では、燃料蒸気の塊（質量M1、体積V1、燃料質量分率＝１）と、筒内ガスの塊（質量M0、体積V0、燃料質量分率＝０）のみが存在する。

ここで、この時点での燃料蒸気の塊の質量M1は、燃料噴射開始時点からの微小時間Δtの間（噴射後経過時間t：０〜Δt）に噴射された燃料の質量に等しい。この時点での筒内ガスの塊の質量M0は、燃料噴射開始時点からの微小時間Δtの間（噴射後経過時間t：０〜Δt）における上記噴霧体積VOLtotal（上記(13)式を参照）の増加部分（この時点では、噴射後経過時間t＝Δtにおける噴霧体積VOLtotal全体に等しい）内に取り込まれた筒内ガスの質量であり、上記「噴霧体積VOLtotalの増加部分」の体積に噴射後経過時間t＝Δtにおける筒内ガスの密度（＝上記ゾーン１内のガスの密度ρz1）を乗じることで求めることができる。

また、本例では、混合気内の各ガス塊の密度は混合気の平均密度ROaveに等しいとの仮定のもと、混合気内の各ガス塊の体積（従って、直径）はそのガス塊の質量を混合気の平均密度ROaveで除することで求められる。従って、例えば、図８において、燃料蒸気の塊の体積V1＝M1/ROaveであり、筒内ガスの塊の体積V0＝M0/ROaveである。後述する新ガス塊１（質量：Mf1）の体積Vf1＝Mf1/ROaveである。混合気の平均密度ROaveは、混合気内に存在する各ガス塊の質量の和をMAStotalとすると、下記(33)式に従って求めることができる。このように、ガス塊の質量が取得できれば、そのガス塊の体積（従って、直径）も取得することができる。

図８の右側の図は、この状態（即ち、噴射後経過時間t＝Δt（衝突前）の状態）にて存在している燃料質量分率が異なる複数の（ここでは、２つの）ガス塊のうちの任意の２つの組み合わせの全てについて上記「衝突反応」が発生した後における（即ち、噴射後経過時間t＝Δt（衝突後）における）混合気内のガス塊の分布（質量分布）を示している。

この場合、衝突前にて混合気内には筒内ガスの塊と燃料蒸気の塊しか存在しないから、「筒内ガスの塊と燃料蒸気の塊についての衝突反応」のみが発生する。この「筒内ガスの塊と燃料蒸気の塊についての衝突反応」では、先ず、筒内ガスの塊のうちの一部分と燃料蒸気の塊のうちの一部分とが混合する。この混合に係わる部分の体積（以下、「混合体積VOLmix」と称呼する。）の計算方法については後述する。この混合体積VOLmixは衝突反応毎に計算される。

次に、図８の左側の図に示した筒内ガスの塊から上記混合体積VOLmixに対応する部分だけが離脱するとともに、図８の左側の図に示した燃料蒸気の塊から上記混合体積VOLmixに対応する部分だけが離脱する。以下、このように、衝突反応によりガス塊から離脱する部分を「離脱部分」とも称呼する。次いで、上記それぞれの（２つの）離脱部分が混合して燃料質量分率＝f1(０＜f1＜１)の新ガス１の塊（即ち、上記混合形成ガス塊）が新たに発生する。なお、この燃料質量分率f1のような、混合形成ガス塊の燃料質量分率の求め方については後述する（後述する図９に示した燃料質量分率＝f2,f3についても同様）。

従って、この衝突反応後における筒内ガスの塊の質量MO及び燃料蒸気の塊の質量M1（図８の右側の図を参照）は、衝突反応前におけるそれぞれの値（図８の左側の図を参照）に対して上記離脱部分に対応する分だけそれぞれ減少している。一方、新たに発生した燃料質量分率＝f1の新ガス１の塊の質量Mf1は、上記それぞれの離脱部分に対応する質量の和となる。以上より、図８の右側の図に示すように、噴射後経過時間t＝Δt（衝突後）では、混合気中において燃料質量分率が異なる３つのガス塊が存在している。

図９は、上述した図８に示した状態から更に微小時間Δtが経過した時点（噴射後経過時間t＝２Δt）で発生する「衝突反応」の前後の様子を模式的に示している。図９の左側の図に示すように、噴射後経過時間t＝２Δt（衝突前）では、噴射後経過時間t＝Δt（衝突後）と同様、混合気内では、燃料蒸気の塊と筒内ガスの塊に加えて、上述した燃料質量分率＝f1の新ガス１の塊が存在する。

この時点での新ガス１の塊の質量Mf1は上述した図８の右側の図に対応する値に等しい。この時点での燃料蒸気の塊の質量M1は、上述した図８の右側の図に対応する値に対して「噴射後経過時間t：Δt〜２Δtの間に噴射された燃料の質量」を加えた値に等しい。この時点での筒内ガスの塊の質量M0は、上述した図８の右側の図に対応する値に対して「噴射後経過時間t：Δt〜２Δtの間における上記噴霧体積VOLtotalの増加部分内に取り込まれた筒内ガスの質量（＝「噴霧体積VOLtotalの増加部分」の体積に噴射後経過時間t＝２Δtにおける筒内ガスの密度を乗じた値）」を加えた値に等しい。

図９の右側の図は、この状態（即ち、噴射後経過時間t＝２Δt（衝突前）の状態）にて存在している燃料質量分率が異なる複数の（ここでは、３つの）ガス塊のうちの任意の２つの組み合わせの全てについて上記「衝突反応」が発生した後における（即ち、噴射後経過時間t＝２Δt（衝突後）における）混合気内のガス塊の分布（質量分布）を示している。

この場合、衝突前にて混合気内には、筒内ガスの塊と燃料蒸気の塊と新ガス１の塊との３つのガス塊が存在するから、「筒内ガスの塊と燃料蒸気の塊についての衝突反応」、「筒内ガスの塊と新ガス１の塊についての衝突反応」、「燃料蒸気の塊と新ガス１の塊についての衝突反応」の３つの衝突反応がそれぞれ発生する。

係る３つの衝突反応のうち、「筒内ガスの塊と燃料蒸気の塊についての衝突反応」に関しては、上述した図８に示した「筒内ガスの塊と燃料蒸気の塊についての衝突反応」と同様、筒内ガスの塊の質量MO及び燃料蒸気の塊の質量M1は、対応するガス塊からの上記離脱部分に対応する分だけそれぞれ減少し、上記燃料質量分率＝f1の新ガス１の塊の質量Mf1は、筒内ガスの塊と燃料蒸気の塊からのそれぞれの離脱部分に対応する質量の和だけ増大する。

「筒内ガスの塊と新ガス１の塊についての衝突反応」では、筒内ガスの塊の離脱部分と、新ガス１の塊の離脱部分とが混合して燃料質量分率＝f2（０＜f2＜f1）の新ガス２の塊が新たに発生する。この結果、この衝突反応に関しては、筒内ガスの塊の質量MO及び新ガス１の塊の質量Mf1は、対応するガス塊からの離脱部分に対応する分だけそれぞれ減少し、新たに発生した燃料質量分率＝f2の新ガス２の塊の質量Mf2は、筒内ガスの塊と新ガス１の塊からのそれぞれの離脱部分に対応する質量の和となる。

「燃料蒸気の塊と新ガス１の塊についての衝突反応」では、燃料蒸気の塊の離脱部分と、新ガス１の塊の離脱部分とが混合して燃料質量分率＝f3（f1＜f3＜１）の新ガス３の塊が新たに発生する。この結果、この衝突反応に関しては、燃料蒸気の塊の質量M1及び新ガス１の塊の質量Mf1は、対応するガス塊からの離脱部分に対応する分だけそれぞれ減少し、新たに発生した燃料質量分率＝f3の新ガス３の塊の質量Mf3は、燃料蒸気の塊と新ガス１の塊からのそれぞれの離脱部分に対応する質量の和となる。以上、係る３つの衝突反応により、図９の右側の図に示すように、噴射後経過時間t＝２Δt（衝突後）では、混合気中において燃料質量分率が異なる５つのガス塊が存在している。

以降も、微小時間Δtの経過毎に（微小クランク角度ΔCAの進行毎に）、上記と同様、以下に示す一連の処理（処理１、処理２）が行われて行く。
（処理１）燃料蒸気の塊の質量M1に微小時間Δtの間において新たに噴射された燃料の質量が加えられるとともに、筒内ガスの塊の質量M0に微小時間Δtの間において混合気中に新たに取り込まれた筒内ガスの質量が加えられる。
（処理２）処理１の後、その時点で存在している燃料質量分率が異なる複数のガス塊（筒内ガスの塊と燃料蒸気の塊とを含む。）のうちの任意の２つの組み合わせの全てについて上記衝突反応が混合気内で順に行われる。

これにより、混合気内に存在する（質量が「０」でない）燃料質量分率が異なるガス塊の個数が増大していくとともに（例えば、図８の右側の図では３つ、図９の右側の図では５つ）、混合気内のガス塊の質量分布が取得・更新されていく。以上が、「衝突反応」を利用して混合気内のガス塊の質量分布を取得・更新していく具体的な方法である。

次に、係る混合気内のガス塊の質量分布の計算に必要となる、「衝突反応」に係わる種々の値の求め方について説明する。以下、説明の便宜上、或る衝突反応Ｚについて、混合気内における２つのガス塊Ａ，Ｂが衝突対象ガス塊となり、ガス塊Ｃが混合形成ガス塊となるものとする。加えて、衝突反応Ｚの前におけるガス塊Ａの質量をMASA、速度をVELA、体積をVOLA（直径をDIAA）、燃料質量分率をFRACmasAとし、衝突反応Ｚの前におけるガス塊Ｂの質量をMASB、速度をVELB、体積をVOLB（直径をDIAB）、燃料質量分率をFRACmasBとし、衝突反応Ｚの前におけるガス塊Ｃの質量をMASC、速度をVELC、体積をVOLC（直径をDIAC）、燃料質量分率をFRACmasC（「FRACmix」とも称呼する。）とする。

<混合形成ガス塊の燃料質量分率FRACmix>
先ず、上記衝突反応Ｚにおける混合形成ガス塊の燃料質量分率FRACmasC（＝FRACmix）の計算（特定）方法について説明する。この混合形成ガス塊の燃料質量分率FRACmixは、「衝突対象ガス塊Ａ，Ｂの離脱部分に含まれるそれぞれの燃料の質量の和が、衝突対象ガス塊Ａ，Ｂの離脱部分のそれぞれの質量の和に混合形成ガス塊Ｃの燃料質量分率FRACmixを乗じた値に等しい」との燃料についての質量保存則を利用して、下記(34)式に従って求めることができる。

(34)式において、dMASA，dMASBはそれぞれ、ガス塊Ａの離脱部分の質量，ガス塊Ｂの離脱部分の質量である。ガス塊Ａ，Ｂの離脱部分の質量dMASA，dMASBは、上記混合体積VOLmixに上記混合気の平均密度ROaveを乗じることで求めることができる。(34)式において、「dMASA・FRACmasA」はガス塊Ａの離脱部分に含まれる燃料の質量に相当し、「dMASB・FRACmasB」はガス塊Ｂの離脱部分に含まれる燃料の質量に相当する。(34)式により、衝突反応毎に混合形成ガス塊の燃料質量分率FRACmixが計算される。

<ガス塊の質量MAS>
次に、ガス塊の質量MASの更新方法について説明する。

<<衝突反応に係わるガス塊の質量>>
上記衝突反応Ｚにより、ガス塊Ａから上記質量dMASAの離脱部分が離脱し、ガス塊Ｂから上記質量dMASBの離脱部分が離脱する。一方、上記ガス塊Ａ，Ｂの離脱部分は共に、ガス塊Ｃの一部（又は全部）としてガス塊Ｃに加わる。従って、衝突反応Ｚの後におけるガス塊Ａ，Ｂ，Ｃの質量MASA’，MASB’，MASC’は下記(35)式に従って求めることができる。(35)式により、衝突反応Ｚに係わるガス塊の質量MAS（従って、体積VOL）を更新することができる。従って、(35)式を衝突反応毎に適用していくことで衝突反応に係わるガス塊の質量を更新していくことができる。

<<筒内ガスの塊と燃料蒸気の塊の質量>>
混合気中において、筒内ガスの塊（即ち、燃料質量分率FRACmas＝０のガス塊）と、燃料蒸気の塊（即ち、燃料質量分率FRACmas＝１のガス塊）を除いた他のガス塊（即ち、０＜燃料質量分率FRACmas＜１であるガス塊）の質量は、そのガス塊が衝突反応に係わることによってのみ上記(35)式に従って変化していく。

一方、筒内ガスの塊の質量は、筒内ガスの塊が衝突反応に係わる場合に加えて、微小時間Δt毎に（微小クランク角度ΔCAの進行毎に）新たに混合気中に取り込まれた筒内ガス（以下、「新規筒内ガス」と称呼する。）が筒内ガスの塊の一部に加わることによっても変化し得る。

即ち、新規筒内ガスが加えられる前における筒内ガスの塊の質量をMAS0とし、新規筒内ガスの質量をΔGmとし、新規筒内ガスが加えられた後における筒内ガスの塊の質量をMAS0’とすると、新規筒内ガスが加えられた後における筒内ガスの塊の質量MAS0’は、下記(36)式に従って求めることができる。従って、(36)式を微小時間Δt毎に適用していくことで、新規筒内ガスが筒内ガスの塊に加えられることに関して、筒内ガスの塊の質量を微小時間Δt毎に更新していくことができる。

他方、燃料蒸気の塊の質量も、筒内ガスの塊の質量と同様、燃料蒸気の塊が衝突反応に係わる場合に加えて、微小時間Δt毎に（微小クランク角度ΔCAの進行毎に）新たに噴射された燃料（以下、「新規燃料」と称呼する。）が燃料蒸気の塊の一部に加わることによっても変化し得る。

即ち、新規燃料が加えられる前における燃料蒸気の塊の質量をMAS1とし、新規燃料の質量をqとし、新規燃料が加えられた後における燃料蒸気の塊の質量をMAS1’とすると、新規燃料が加えられた後における燃料蒸気の塊の質量MAS1’は、下記(37)式に従って求めることができる。従って、(37)式を微小時間Δt毎に適用していくことで、新規燃料が燃料蒸気の塊に加えられることに関して、燃料蒸気の塊の質量を微小時間Δt毎に更新していくことができる。

以上のようにして、(36)式、(37)式を微小時間Δt毎に適用していくとともに、(35)式を衝突反応毎に適用していくことで、混合気内のガス塊の質量分布を順次取得・更新していくことができる。

<混合体積VOLmix>
次に、混合気内のガス塊の質量分布を取得するために必要となる混合体積VOLmixの計算方法について説明する。混合気体積VOLmixは、衝突した２つのガス塊Ａ，Ｂが係わりあった体積に等しいと考えることができる。いま、衝突する２つのガス塊Ａ，Ｂのうち体積（直径）が小さい方のガス塊（小ガス塊）の全体が体積（直径）が大きい方のガス塊（大ガス塊）の内部を通過するものと仮定する。

この場合、微小時間Δt（微小クランク角度ΔCA）の間に２つのガス塊Ａ，Ｂが係わりあう体積（以下、「通過体積VOLswp」と称呼する。）は、図１０に示すようにカプセル状の体積となり、下記(38)式に従って計算することができる。(38)式において、DIAsmallは小ガス塊の直径である。｜VELA−VELB｜はガス塊Ａ，Ｂの相対速度に相当する。

混合体積VOLmixは、上記通過体積VOLswpに大きく依存すると考えられる。他方、小ガス塊全体が大ガス塊の内部を通過する割合（確率）は、大ガス塊の体積VOLbig（＝VOLAとVOLBのうちの大きい方）の混合気全体の体積（即ち、噴霧体積VOLtotal）に対する割合に等しいと仮定することができる。

更には、エンジン回転速度NEが大きいほど（即ち、微小クランク角度ΔCAの進行に対応する微小時間Δtが短いほど）ガス塊Ａ，Ｂの混合に供され得る時間が短くなるから、混合体積VOLmixはエンジン回転速度NEが大きいほど小さくなると考えることができる。

以上のことから、混合気体積VOLmixは、下記(39)式に従って表すことができる。(39)式において、αは係数であり、エンジン回転速度NEが大きいほどより小さい値に設定される。下記(39)式に従って、衝突反応毎に混合体積VOLmixが計算される。

<ガス塊の速度VEL>
次に、上記混合体積VOLmixの算出に必要となるガス塊の速度VELの取得・更新方法について説明する。

<<混合形成ガス塊の速度>>
上記衝突反応Ｚにより、上記ガス塊Ａ，Ｂの離脱部分がガス塊Ｃに加わる。これに伴い、混合形成ガス塊であるガス塊Ｃの速度は、衝突反応Ｚの前における値VELCから変化する。衝突反応Ｚの後におけるガス塊Ｃの速度VELC’は、「衝突反応Ｚの前における混合形成ガス塊Ｃの運動量と衝突反応Ｚの前における衝突対象ガス塊Ａ，Ｂの離脱部分のそれぞれの運動量との和が、衝突反応Ｚの後における混合形成ガス塊Ｃの運動量に等しい」との運動量保存則を利用して、下記(40)式に従って求めることができる。

(40)式において、分母全体は、衝突反応Ｚの後におけるガス塊Ｃの質量（＝MASC’）である（上記(35)式を参照）。分子において、「MASC・VELC」は衝突反応Ｚの前におけるガス塊Ｃの運動量に相当し、「dMASA・VELA」は衝突反応Ｚの前におけるガス塊Ａの離脱部分の運動量に相当し、「dMASB・VELB」は衝突反応Ｚの前におけるガス塊Ｂの離脱部分の運動量に相当する。上記(40)式により、衝突反応Ｚにおける混合形成ガス塊Ｃの速度VELCを更新することができる。従って、(40)式を衝突反応毎に適用していくことで衝突反応毎に、混合形成ガス塊となったガス塊の速度を更新していくことができる。

<<筒内ガスの塊と燃料蒸気の塊の速度>>
混合気中において、筒内ガスの塊と、燃料蒸気の塊を除いた他のガス塊（即ち、０＜燃料質量分率FRACmas＜１であるガス塊）の速度は、そのガス塊が衝突反応における混合形成ガス塊となることによってのみ上記(40)式に従って変化していく。

一方、筒内ガスの塊の速度は、筒内ガスの塊が衝突反応における混合形成ガス塊となった場合に加えて、微小時間Δt毎に（微小クランク角度ΔCAの進行毎に）新たに混合気中に取り込まれた筒内ガス（即ち、上記新規筒内ガス）が筒内ガスの塊の一部に加わることによっても変化し得る。

ここで、新規筒内ガスが加えられる前における筒内ガスの塊の質量及び速度をそれぞれMAS0，VEL0とし、新規筒内ガスの質量及び速度をそれぞれΔGm，VELgasとし、新規筒内ガスが加えられた後における筒内ガスの塊の速度をVEL0’とすると、新規筒内ガスが加えられた後における筒内ガスの塊の速度VEL0’は、「新規筒内ガスが加えられる前における筒内ガスの塊の運動量と新規筒内ガスの運動量との和が、新規筒内ガスが加えられた後における筒内ガスの塊の運動量に等しい」との運動量保存則を利用して、下記(41)式に従って求めることができる。

(41)式において、分母全体は、新規筒内ガスが加えられた後における筒内ガスの塊の質量である。分子において、「MAS0・VEL0」は新規筒内ガスが加えられる前における筒内ガスの塊の運動量に相当し、「ΔGm・VELgas」は新規筒内ガスの運動量に相当する。なお、本例では、新規筒内ガスの速度VELgas＝０とする。これは、筒内ガスが混合気内に取り込まれていくのは、燃焼室内において静止している筒内ガス中において燃料噴霧の領域（従って、混合気の領域）が広がっていくからであると考えることができることに基づく。

以上、(41)式を微小時間Δt毎に適用していくことで、新規筒内ガスが筒内ガスの塊に加えられることに関して、筒内ガスの塊の速度を微小時間Δt毎に更新していくことができる。

他方、燃料蒸気の塊の速度も、筒内ガスの塊の速度と同様、燃料蒸気の塊が衝突反応における混合形成ガス塊となった場合に加えて、微小時間Δt毎に（微小クランク角度ΔCAの進行毎に）新たに噴射された燃料（即ち、上記新規燃料）が燃料蒸気の塊の一部に加わることによっても変化し得る。

ここで、新規燃料が加えられる前における燃料蒸気の塊の質量及び速度をそれぞれMAS1，VEL1とし、新規燃料の質量及び速度をそれぞれq，VELinjとし、新規燃料が加えられた後における燃料蒸気の塊の速度をVEL1’とすると、新規燃料が加えられた後における燃料蒸気の塊の速度VEL1’は、「新規燃料が加えられる前における燃料蒸気の塊の運動量と新規燃料の運動量との和が、新規燃料が加えられた後における燃料蒸気の塊の運動量に等しい」との運動量保存則を利用して、下記(42)式に従って求めることができる。

(42)式において、分母全体は、新規燃料が加えられた後における燃料蒸気の塊の質量である。分子において、「MAS1・VEL1」は新規燃料が加えられる前における燃料蒸気の塊の運動量に相当し、「q・VELinj」は新規燃料の運動量に相当する。ここで、新規燃料の速度VELinjは噴射された燃料の速度に等しいと考えることができ、噴射された燃料の速度は上記有効噴射圧力ΔP0と液体燃料の密度ρfに大きく依存する。従って、新規燃料の速度VELinjは、ΔP0とρfを引数とする新規燃料の速度VELinjを求める関数funcVELinjを利用して求めることができる。

以上、(42)式を微小時間Δt毎に適用していくことで、新規燃料が燃料蒸気の塊に加えられることに関して、燃料蒸気の塊の速度を微小時間Δt毎に更新していくことができる。以上、混合気内のガス塊の質量分布の計算に必要となる、「衝突反応」に係わる種々の値の求め方について説明した。

<燃料質量分率の離散的な設定に対する対処>
次に、混合気を構成するガス塊の燃料質量分率FRACmasの設定について説明する。上記(34)式から理解できるように、混合形成ガス塊の燃料質量分率FRACmixは、「０」から「１」の範囲内における任意の有理数に計算され得る値である。しかしながら、混合気を構成するガス塊の燃料質量分率FRACmasが「０」から「１」の範囲内で任意の有理数に設定されることを許容すると、混合気を構成するガス塊の個数が無限個に増大し得、この結果、計算負荷が膨大なものとなる可能性が発生する。

以上のことから、本装置では、混合気を構成するガス塊の燃料質量分率FRACmasが「０」から「１」の範囲内で一定の間隔をもって離散的に設定される。以下、このことを、図１１を参照しながら説明する。図１１は、本装置により以上のような処理を燃料噴射開始からの微小時間Δtの経過毎に（微小クランク角度ΔCAの進行毎に）繰り返した場合における混合気内のガス塊の質量分布の一例を示している。

図１１に示すように、本装置では、混合気内において存在し得るガス塊の個数が「Nd＋1」個（Ndは自然数）に限定され、「Nd＋1」個のガス塊のそれぞれに対して区分番号i（i＝0,1,2,・・・,Nd-1,Nd）が順に付される。以下、区分番号iのガス塊を「ガス塊(i)」と称呼し、ガス塊(i)の燃料質量分率を「FRACmas(i)」と称呼するものとすると、ガス塊(i)の燃料質量分率FRACmas(i)は、値「i/Nd」に設定される（i＝0,1,2,・・・,Nd-1,Nd）。このように、混合気を構成するガス塊の燃料質量分率FRACmasが「０」から「１」の範囲内で値「1/Nd」毎に離散的に設定される。

このように、混合気を構成するガス塊の燃料質量分率FRACmasが「０」から「１」の範囲内で離散的に設定されると、上記(34)式に従って計算される混合形成ガス塊の燃料質量分率FRACmixの値と完全に一致する燃料質量分率FRACmasを有するガス塊が存在しない事態が発生し得る。

係る事態が発生する場合、本装置では、上記(34)式に従って計算される混合形成ガス塊の燃料質量分率FRACmixの値を挟む２つの燃料質量分率FRACmas(md)，FRACmas(mu)が特定される(mu＝md＋1)。そして、混合形成ガス塊として、ガス塊(md)とガス塊(mu)の２つのガス塊が使用される。以下、このことを、図１２を参照しながらより具体的に説明する。

図１２は、燃料噴射開始後の或る時点で発生する「ガス塊(mb)とガス塊(ms)についての衝突反応」の前後の様子を模式的に示している（mb,ms：０以上Nd以下の整数）。この衝突反応について上記(34)式に従って計算される混合形成ガス塊の燃料質量分率FRACmixの値は、燃料質量分率FRACmas(md)とFRACmas(mu)の間に入るものとする。なお、理解を容易にするため、図１２では、この衝突反応に係わる４つのガス塊（ガス塊(mb)、ガス塊(ms)、ガス塊(md)、ガス塊(mu)）以外のガス塊は表示されていない。MAS(i)，VOL(i)はそれぞれ、ガス塊(i)の質量、体積である（i＝0,1,2,・・・,Nd-1,Nd）。

この場合、図１２に示すように、ガス塊(mb)及びガス塊(ms)からの各離脱部分は、混合形成ガス塊となるガス塊(md)とガス塊(mu)の２つのガス塊に所定の配分割合をもって移行される。この配分割合（ガス塊(md)への配分と、ガス塊(mu)への配分の割合）は、（値FRACmas(mu)−値FRACmix）：（値FRACmix−値FRACmas(md)）に決定される。係る処理の詳細については後にフローチャートを参照しながら説明する。

<混合気の不均一度を表す燃料質量分率の標準偏差σ>
本装置は、上述した手法により、燃料噴射開始時点以降において、微小時間Δtの経過毎に（微小クランク角度ΔCAの進行毎に）、混合気内のガス塊の質量分布を順次取得・更新していく。これに伴い、本装置は、下記(43)式に従って、燃料噴射開始時点以降における微小時間Δtの経過毎に、「混合気内における燃料質量分率の標準偏差σ」を、上記取得した混合気内のガス塊の質量分布に基づいて順次取得・更新していく。

(43)式において、FRACmasaveは、混合気内における燃料質量分率の平均値であり、下記（44）式に従って求めることができる。

このようにして求められる標準偏差σは、混合気中における燃料濃度の不均一性を表す値となる。従って、上記問題点２を解消するためには、この標準偏差σにより表される混合気の不均一度を上述した燃焼速度qrz2,qrz3(g/sec)（従って、熱発生率Hrz2,Hrz3(J/deg)）の計算に反映させればよい。

このため、本装置は、図１３に示すように設定される補正係数Ａを導入する。この補正係数Ａは、標準偏差σが値σ1（定数）以下では「１」に固定され、標準偏差σが値σ1より大きい場合、標準偏差σが大きいほどより小さい値に設定される。

そして、本装置は、qrz2(g/sec)（従って、熱発生率Hrz2(J/deg)）を、上記(20)式において補正係数Ａが考慮された下記(45)式に従って計算し、qrz3(g/sec)（従って、熱発生率Hrz3(J/deg)）を、上記(31)式において補正係数Ａが考慮された下記(46)式に従って計算する。

これにより、標準偏差σが大きくなるほど（従って、混合気の不均一度が大きくなるほど）、qrz2,qrz3(g/sec)（従って、熱発生率Hrz2,Hrz3(J/deg)）が小さめに計算される。これにより、本発明者は、上記問題点２が解消されることを見出した。

図１４は、上述した図６に対応するグラフであって、上記問題点１を「着火後噴射継続の場合において前記着火前噴射部分（ゾーン２）と前記着火後噴射部分（ゾーン３）とを分けて取り扱うこと」で解消し、上記問題点２を「混合気内における燃料濃度の不均一性（標準偏差σ）を考慮すること」で解消した後における計算結果（実線を参照）を示している。

図１４から理解できるように、この計算結果では、実験結果と同様、１回目の比較的大きいピークの後で２回目の比較的小さいピーク
が発生している（Ａを参照）。即ち、上記問題点１が解消されている。この計算結果において、１回目の大きいピークは、ゾーン２における燃焼速度qrz2(g/sec)（従って、熱発生率Hrz2(J/deg)）のピークに対応し、２回目の小さいピークは、ゾーン３における燃焼速度qrz3(g/sec)（従って、熱発生率Hrz3(J/deg)）のピークに対応しているものと考えられる。

加えて、図１４から理解できるように、この計算結果における１回目のピークは実験結果における１回目のピークに比して過大とはなっていない（Ｂを参照）。即ち、上記問題点２が解消されている。これは、上記標準偏差σにより表される混合気の不均一度が燃焼速度qrz2,qrz3(g/sec)（従って、熱発生率Hrz2,Hrz3(J/deg)）の計算に反映されたことによるものであると考えられる。以上、本装置による混合気状態の推定方法、クランク角度CAに対する熱発生率の推定方法、並びに混合気の不均一度の推定方法の概要について説明した。

（実際の作動）
次に、上記のように構成された内燃機関の混合気の不均一度取得装置を含んだ混合気状態取得装置の実際の作動について説明する。

ＣＰＵ６１は、図１５〜図１６にフローチャートにより示した一連のルーチン（メインルーチン）を、吸気弁Vinが閉弁する毎に（即ち、ＩＶＣが到来する毎に）、ＩＶＣが到来した気筒について実行するようになっている。このメインルーチンの実行により、クランク角度CAinj〜クランク角度CAend（例えば、膨張下死点近傍）の間におけるクランク角度CA（CA＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)(j＝1,2,・・・)が計算される。このメインルーチンの実行はＩＶＣの直後に完了するから、上記熱発生率Hr(j)
(j＝1,2,・・・)もＩＶＣの直後（従って、実際の燃料噴射時期CAinjの到来前）に取得され得る。以下、係る作動について詳細に説明していく。

或る気筒についてＩＶＣが到来すると、ＣＰＵ６１はメインルーチンの処理を開始してステップ１５０５に進み、同ステップ１５０５を経由して図１７〜図１８にフローチャートにより示した「初期値の設定」を行う一連のルーチン（サブルーチン）の処理をステップ１７００から開始する。

ＣＰＵ６１は、ステップ１７００からステップ１７０５に進むと、ＩＶＣ時クランク角度CAivcをクランクポジションセンサ７４から取得される現時点での実クランク角度CAactの値に設定し、ＩＶＣ時筒内ガス圧力Pgivcを吸気管圧力センサ７３から得られる現時点での吸気管圧力Pbの値に設定し、ＩＶＣ時筒内ガス温度Tgivcを吸気温センサ７２から得られる現時点での吸気温度Tbの値に設定し、吸気酸素濃度[O2]inを吸気酸素濃度センサ７６から得られる現時点での吸気酸素濃度RO2inの値に設定する。

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

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

次に、ＣＰＵ６１はステップ１７２０に進み、燃料噴射量Qfin、エンジン回転速度NE、及び図３１に示したテーブルMapCAinjから燃料噴射時期CAinjを決定する。テーブルMapCAinjは、燃料噴射量Qfin及びエンジン回転速度NEと燃料噴射時期CAinjとの関係を規定するテーブルであり、ＲＯＭ６２内に格納されている。

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

次に、ＣＰＵ６１はステップ１７３０に進み、上記求めた燃料噴射期間TAUと、エンジン回転速度NEとから、燃料噴射期間TAUを微小クランク角度ΔCA（一定）に対応する微小時間Δt（エンジン回転速度NEに依存する）で分割した場合の分割数qdivを決定する。

続いて、ＣＰＵ６１はステップ１７３５に進んで、上記分割数qdivと、上記基本燃料噴射圧力Pcrbaseとから、燃料噴射期間TAUに亘る微小クランク角度ΔCA毎の燃料噴射圧力の分布、具体的にはPcr(j)(j＝1,・・・,qdiv)を求める。ここで、Pcr(j)
(j＝1,・・・,qdiv)は、クランク角度（CAinj＋(j-1)・ΔCA）〜クランク角度（CAinj＋j・ΔCA）の間に噴射される燃料の噴射圧力の平均値である。

次に、ＣＰＵ６１はステップ１７４０に進んで、上記分割数qdivと、上記基本燃料噴射圧力Pcrbaseとから、燃料噴射期間TAUに亘る微小クランク角度ΔCA毎の燃料噴射量の分布、具体的にはq(j)(j＝1,・・・,qdiv)を求める。ここで、q(j)
(j＝1,・・・,qdiv)は、クランク角度（CAinj＋(j-1)・ΔCA）〜クランク角度（CAinj＋j・ΔCA）の間に噴射される燃料の量(g)である。なお、「Qfin＝q(1)＋・・・＋q(qdiv)」の関係が成立する。

次に、ＣＰＵ６１はステップ１７４５に進み、エンジン回転速度NEと、微小クランク角度ΔCAとから、微小クランク角度ΔCAに対応する微小時間Δt（微小クランク角度ΔCAの進行に要する時間）を決定する。

続いて、ＣＰＵ６１はステップ１７５０に進んで、燃焼室内圧力Pg（の初期値）を上記ＩＶＣ時筒内ガス圧力Pgivcと等しい値に設定するとともに、ゾーン１内のガス温度Tz1（の初期値）を上記ＩＶＣ時筒内ガス温度Tgivcと等しい値に設定し、続くステップ１７５５にてガス塊(i)の燃料質量分率FRACmas(i)の値を値「i/Nd」に設定する(i＝0,1,・・・,Nd-1,Nd)。

次に、ＣＰＵ６１は図１８のステップ１８０５に進み、ゾーン２の体積の前回値Vbz2（の初期値）、ゾーン３の体積の前回値Vbz3（の初期値）、及び噴霧体積の前回値VOLbtotal（の初期値）を「０」に設定し、続くステップ１８１０にて、ガス塊(i)の質量MAS(i)（の初期値）、及びガス塊(i)の速度VEL(i)（の初期値）を「０」に設定する(i＝0,1,・・・,Nd-1,Nd)。

次いで、ＣＰＵ６１はステップ１８１５に進んで、ゾーン２内筒内ガス質量Gz2（の初期値）、及びゾーン３内筒内ガス質量Gz3（の初期値）を「０」に設定し、続くステップ１８２０にてゾーン２内において既に噴射された燃料の量（ゾーン２内噴射燃料量）Qz2（の初期値）、及びゾーン３内において既に噴射された燃料の量（ゾーン３内噴射燃料量）Qz3（の初期値）を「０」に設定する。

次に、ＣＰＵ６１はステップ１８２５に進み、ゾーン２内の消費燃料量sumqrz2（の初期値）、及びゾーン３内の消費燃料量sumqrz3（の初期値）を「０」に設定し、続くステップ１８３０にて標準偏差σ（の初期値）を「０」に設定した後、ステップ１８９５を経由して図１５のステップ１５１０に復帰する。以上により、計算に使用される変数等の初期値の設定が完了する。

ＣＰＵ６１はステップ１５１０に進むと、計算上のクランク角度CA（以下、単に「クランク角度CA」とも称呼する。）を上記ＩＶＣ時クランク角度CAivcに設定し、変数jの値（自然数）を「１」に設定する。ここで、変数jの値は、計算上のクランク角度CAにおける燃料噴射時期CAinjからの微小クランク角度ΔCAの進行回数を表す（CA＝CAinj＋j・ΔCA）。

次いで、ＣＰＵ６１はステップ１５１５に進み、フラグＸＩＧの値、及びフラグＸＺ３の値を共に「０」に設定する。ここで、フラグＸＩＧの値は、その値が「１」のとき「着火後」を表し、その値が「０」のとき「着火前」を表す。フラグＸＺ３の値は、その値が「１」のとき着火後においてゾーン３が発生する状態（即ち、着火後噴射継続の場合）を表し、その値が「０」のとき着火後においてゾーン３が発生しない状態（即ち、着火後噴射継続でない場合）を表す。

続いて、ＣＰＵ６１はステップ１５２０に進んで、同ステップ１５２０を経由して図１９にフローチャートにより示した「ゾーン１の計算」を行うルーチン（サブルーチン）の処理をステップ１９００から開始する。

ＣＰＵ６１は、ステップ１９００からステップ１９０５に進むと、ゾーン１内のガス質量Gz1を、先のステップ１７１０にて求めた筒内ガスの全質量Gcと、ゾーン２内筒内ガス質量Gz2及びゾーン３内筒内ガス質量Gz3（現時点では、これらは先のステップ１８１５の処理により「０」である）と、上記(29)式とに基づいて求める。従って、現時点では、ゾーン１内のガス質量Gz1は筒内ガスの全質量Gcと等しい値に設定される。

次に、ＣＰＵ６１はステップ１９１０に進み、ゾーン１の体積Vz1を、ゾーン２の体積の前回値Vbz2及びゾーン３の体積の前回値Vbz3（現時点では、これらは先のステップ１８０５の処理により「０」である）と、上記(30)式に相当する式とに基づいて求める。従って、現時点では、ゾーン１の体積Vz1は筒内容積Vc(CA)と等しい値に設定される。

続いて、ＣＰＵ６１はステップ１９１５に進んで、上記ゾーン１内のガス質量Gz1を上記ゾーン１の体積Vz1で除してゾーン１内のガス密度ρz1を求める。次いで、ＣＰＵ６１はステップ１９２０に進み、ゾーン１に対する微小ピストン仕事ΔWpistonz1を、燃焼室内圧力Pg（現時点では、先のステップ１７５０の処理により上記ＩＶＣ時筒内ガス圧力Pgivcである。）と、上記ゾーン１の体積Vz1と、上記(11)式に相当する式とに基づいて求める。

次いで、ＣＰＵ６１はステップ１９２５に進んで、ゾーン１に対する微小エネルギーΔWz1を上記ゾーン１に対する微小ピストン仕事ΔWpistonz1と等しい値に設定する。次に、ＣＰＵ６１はステップ１９３０に進み、ゾーン１内のガスの定圧比熱Cpz1を、上記ステップ１７０５にて設定された吸気酸素濃度[O2]inと、排気酸素濃度センサ７７から得られる前回の排気行程中における所定時期での排気酸素濃度[O2]exbと、ゾーン１内のガス温度Tz1と、上記(7)式とに基づいて求める。Tz1としては、今回（初回）に限り先のステップ１７５０の処理によりＩＶＣ時筒内ガス温度Tgivcと等しい値が使用されるが、次回からは、後述するステップ１９４０にて１回前の本ルーチン実行時にて更新されている値が使用される。

続いて、ＣＰＵ６１はステップ１９３５に進んで、ゾーン１内のガス温度の増大量ΔTz1を、先のステップ１９２５にて得られた微小エネルギーΔWz1と、先のステップ１９０５にて求めたゾーン１内のガス質量Gz1と、上記ゾーン１内のガスの定圧比熱Cpz1とに基づいて求める。このゾーン１内のガス温度の増大量ΔTz1は、上記(4)式の右辺第２項に対応している。次いで、ＣＰＵ６１はステップ１９４０に進み、ゾーン１内のガス温度Tz1を、現時点での値に上記増大量ΔTz1を加えることで更新する。

次に、ＣＰＵ６１はステップ１９４５に進んで、ゾーン１内のガス圧力Pz1を、上記ゾーン１内のガス質量Gz1と、上記更新されたゾーン１内のガス温度Tz1と、先のステップ１９１０にて求めたゾーン１の体積Vz1と、ゾーン１内のガスについての気体の状態方程式（ステップ１９４５内の式を参照）とに基づいて求める。

そして、ＣＰＵ６１はステップ１９５０に進み、ゾーン１の体積の前回値Vbz1を上記ゾーン１の体積Vz1に設定した後、ステップ１９９５を経由して図１５のステップ１５２５に復帰する。以上により、クランク角度CA（現時点では、先のステップ１５１０の処理によりＩＶＣ時クランク角度CAivcと等しい）におけるゾーン１内のガスの各物理量が計算・更新される。

ＣＰＵ６１はステップ１５２５に進むと、クランク角度CAが先のステップ１７２０にて設定された燃料噴射時期CAinj以降となっているか否かを判定する。現時点では、クランク角度CAは、ＩＶＣ時クランク角度CAivcと等しいから燃料噴射時期CAinj以前である。従って、ＣＰＵ６１はステップ１５２５にて「Ｎｏ」と判定してステップ１５３０に進み、燃焼室内圧力Pgを先のステップ１９４５にて求めたゾーン１内のガス圧力Pz1と等しい値に設定する。

次に、ＣＰＵ６１は図１６のステップ１６２０に進み、クランク角度CAを現時点での値（＝CAivc）に微小クランク角度ΔCAを加えることで更新する。続いて、ＣＰＵ６１はステップ１６２５に進んで、ステップ１６２０にて更新したクランク角度CAが計算終了クランク角度CAend（例えば、膨張下死点近傍）以降となっているか否かを判定する。

現時点では、クランク角度CAは、ＩＶＣ時クランク角度CAivcと等しいから計算終了クランク角度CAend以前である。従って、ＣＰＵ６１はステップ１６２５にて「Ｎｏ」と判定して図１５のステップ１５２０に戻り、図１９のルーチンを再び実行する。これにより、クランク角度CA（＝CAivc＋ΔCA）におけるゾーン１内のガスの各物理量が計算・更新される。

このような処理（ステップ１５２０、１５２５、１５３０、１６２０、１６２５の処理）は、ステップ１６２０にて更新されていくクランク角度CAが燃料噴射時期CAinjに達しない限りにおいて繰り返し実行される。この結果、図１９のルーチンの繰り返し実行により、ゾーン１内のガスの各物理量がクランク角度CAの更新に応じて計算・更新されていく。

次に、この状態にてステップ１６２０にて更新されていくクランク角度CAが燃料噴射時期CAinjに達した場合について説明する。この場合、ＣＰＵ６１はステップ１５２５に進んだとき「Ｙｅｓ」と判定してステップ１５３５に進むようになり、同ステップ１５３５を経由して図２０〜図２２にフローチャートにより示した「ゾーン２の計算」を行う一連のルーチン（サブルーチン）の処理をステップ２０００から開始する。

ＣＰＵ６１は、ステップ２０００からステップ２００５に進むと、変数jの値が「１」であるか否かを判定する。現時点では、先のステップ１５１０の処理により変数j＝１となっている。従って、ＣＰＵ６１はステップ２００５にて「Ｙｅｓ」と判定してステップ２０１０に進み、ゾーン２内のガス温度Tz2（の初期値）を先のステップ１９４０にて更新されているゾーン１内のガス温度Tz1と等しい値（即ち、燃料噴射開始時点でのゾーン１内のガス温度）に設定する。

続いて、ＣＰＵ６１はステップ２０１５に進んで、燃料噴射開始時点でのゾーン１内のガス密度ρ0z1（即ち、燃料噴射開始時点での筒内ガスの密度）を先のステップ１９１５にて更新されているゾーン１内のガス密度ρz1と等しい値に設定する。

次いで、ＣＰＵ６１はステップ２０２０に進み、有効噴射圧力ΔP0を、先のステップ１７２５にて求めた基本燃料噴射圧力Pcrbaseから燃焼室内圧力Pg（現時点では、先のステップ１５３０の処理によりゾーン１内のガス圧力Pz1と等しい）を減じることで求める。

次に、ＣＰＵ６１はステップ２０２５に進んで、噴霧角θを、上記有効噴射圧力ΔP0と、上記ゾーン１内のガス密度ρ0z1と、上記テーブルMapθとに基づいて求める。続いて、ＣＰＵ６１はステップ２０３０に進み、噴射後経過時間tz2を微小時間Δtに設定し、続くステップ２０３５にて、上記有効噴射圧力ΔP0と、液体燃料の密度ρfと、上記関数funck/εとに基づいて値(k/ε)を求める。以上、ステップ２０１０〜ステップ２０３５にてゾーン２の計算を開始するにあたり必要な各種値が計算される。

次に、ＣＰＵ６１はステップ２０４０に進み、上記有効噴射圧力ΔP0と、噴射後経過時間tz2（現時点では、微小時間Δt）と、上記(12)式に対応する関数funcとに基づいて到達距離Sz2を求め、続くステップ２０４５にて、上記求めた噴霧角θと、到達距離Sz2と、上記(13)式に相当する式とに基づいて噴霧体積VOLtotalを求める。

次いで、ＣＰＵ６１はステップ２０５０に進んで、フラグＸＺ３の値が「０」であるか否かを判定する。現時点では、先のステップ１５１５の処理によりフラグＸＺ３＝０（即ち、ゾーン３が発生しない状態）となっている。従って、ＣＰＵ６１はステップ２０５０にて「Ｙｅｓ」と判定してステップ２０５５に進み、ゾーン２の体積Vz2を上記噴霧体積VOLtotalと等しい値に設定する。

続いて、ＣＰＵ６１はステップ２０６０に進んで、上記求めたゾーン２の体積Vz2からゾーン２の体積の前回値Vbz2を減じることで、微小時間Δtの間におけるゾーン２の体積の増加量ΔVz2を求める。Vbz2としては、今回（初回）に限り先のステップ１８０５の処理により「０」が使用されるが、次回からは、後述するステップ２２４０にて１回前の本ルーチン実行時にて更新されている値が使用される。

次に、ＣＰＵ６１はステップ２０６５に進み、先のステップ１９１５にて更新されているゾーン１内のガス密度ρz1に上記ゾーン２の体積の増加量ΔVz2を乗じることで、微小時間Δtの間におけるゾーン２内筒内ガス質量の増加量ΔGz2を求める。続いて、ＣＰＵ６１はステップ２０７０に進んで、ゾーン２内筒内ガス質量Gz2を、現時点での値（今回（初回）は、先のステップ１８１５の処理により「０」）に上記増加量ΔGz2を加えることで更新する。

続いて、ＣＰＵ６１は図２１のステップ２１０５に進み、フラグＸＺ３の値が「０」であるか否かを判定し、上述のとおり、現時点では「Ｙｅｓ」と判定してステップ２１１０に進んで変数jの値が分割数qdiv以下か否かを判定する。ここで、「j≦qdiv」は、燃料の噴射継続中に対応し、「j＞qdiv」は、燃料の噴射終了後に対応する。

現時点では、変数j＝１であり「j≦qdiv」が成立するから（燃料噴射継続中であるから）、ＣＰＵ６１はステップ２１１０にて「Ｙｅｓ」と判定してステップ２１１５に進み、ゾーン２内噴射燃料量Qz2を、現時点での値（今回（初回）は、先のステップ１８２０の処理により「０」）に先のステップ１７４０にて求めた微小時間Δtの間の燃料量q(j)（現時点では、q(1)）を加えることで更新する。

次いで、ＣＰＵ６１はステップ２１２０に進んで、潜熱Qlatentz2を、上記燃料量q(j)と、q(j)を引数とする潜熱を求める関数funcQlatentとに基づいて求め、続くステップ２１２５にて、ゾーン２内のガス質量Mz2を、先のステップ２０７０にて求めたゾーン２内筒内ガス質量Gz2に先のステップ２１１５にて求めたゾーン２内噴射燃料量Qz2を加えることで求める。

続いて、ＣＰＵ６１はステップ２１３０に進み、ゾーン２内の燃料濃度[Fuel]z2を、上記ゾーン２内噴射燃料量Qz2と、ゾーン２内の消費燃料量sumqrz2と、上記ゾーン２内のガス質量Mz2と、上記(25)式とに基づいて求める。sumqrz2としては、今回（初回）に限り先のステップ１８２５の処理により「０」が使用されるが、次回からは、後述するステップ２２３５にて１回前の本ルーチン実行時にて更新されている値が使用される。

次いで、ＣＰＵ６１はステップ２１３５に進んで、ゾーン２内の酸素濃度[O2]z2を、上記ゾーン２内筒内ガス質量Gz2と、先のステップ１７０５にて設定されている吸気酸素濃度[O2]inと、上記消費燃料量sumqrz2と、上記ゾーン２内のガス質量Mz2と、上記(27)式とに基づいて求める。sumqrz2としては、ステップ２１３０に使用される値と同じ値が使用される。

次に、ＣＰＵ６１はステップ２１４０に進み、フラグＸＩＧの値が「０」であるか否か（従って、「着火前」であるか否か）を判定する。現時点では、先のステップ１５１５の処理により「０」となっている。従って、ＣＰＵ６１はステップ２１４０にて「Ｙｅｓ」と判定してステップ２１４５に進み、上記ゾーン２内の酸素濃度[O2]z2と、上記ゾーン２内の燃料濃度[Fuel]z2と、ゾーン２内のガス温度Tz2と、上記シェル・モデルに対応する関数funcshellとに基づいて、低温酸化反応に基づく予混合燃焼によるゾーン２における燃焼速度qrz2(g/sec)を求める。Tz2としては、今回（初回）に限り先のステップ２０１０の処理により「Tz1」が使用されるが、次回からは、後述するステップ２２４５にて１回前の本ルーチン実行時にて更新されている値が使用される。

続いて、ＣＰＵ６１はステップ２１５０に進み、ゾーン２の熱発生率Hrz2(j)を、上記求めた燃焼速度qrz2と、上記(24)式に相当する式とに基づいて求める。この熱発生率Hrz2(j)は、クランク角度CA（＝CAinj＋j・ΔCA）における値である。

次に、ＣＰＵ６１は図２２のステップ２２０５に進み、ゾーン２に対する微小ピストン仕事ΔWpistonz2を、燃焼室内圧力Pgと、上記ゾーン２の体積Vz2と、上記(19)式に相当する式とに基づいて求める。

次いで、ＣＰＵ６１はステップ２２１０に進んで、ゾーン２に対する微小エネルギーΔWz2を上記ゾーン２に対する微小ピストン仕事ΔWpistonz2と、先のステップ２１２０にて求めた潜熱Qlatentz2と、先のステップ２１５０にて求めた熱発生率Hrz2(j)と、上記(18)式に相当する式とに基づいて求める。

次に、ＣＰＵ６１はステップ２２１５に進み、ゾーン２内のガスの定圧比熱Cpz2を、上記吸気酸素濃度[O2]inと、上記排気酸素濃度[O2]exbと、ゾーン２内のガス温度Tz2と、先のステップ２１３０にて求めたゾーン２内の燃料濃度[Fuel]z2と、上記(17)式とに基づいて求める。Tz2としては、先のステップ２１４５と同じ値が使用される。

続いて、ＣＰＵ６１はステップ２２２０に進んで、ゾーン２内のガス温度の増大量ΔTz2を、先のステップ２２１０にて得られた微小エネルギーΔWz2と、先のステップ２１２５にて求めたゾーン２内のガス質量Mz2と、上記ゾーン２内のガスの定圧比熱Cpz2とに基づいて求める。このゾーン２内のガス温度の増大量ΔTz2は、上記(16)式の右辺第２項に対応している。次いで、ＣＰＵ６１はステップ２２２５に進み、ゾーン２内のガス温度Tz2を、現時点での値に上記増大量ΔTz2を加えることで更新する。

次に、ＣＰＵ６１はステップ２２３０に進んで、ゾーン２内のガス圧力Pz2を、上記ゾーン２内のガス質量Mz2と、上記更新されたゾーン２内のガス温度Tz2と、先のステップ２０５５にて求めたゾーン２の体積Vz2と、ゾーン２内のガスについての気体の状態方程式（ステップ２２３０内の式を参照）とに基づいて求める。このゾーン２内のガス圧力Pz2は先のステップ１９４５にて更新されているゾーン１内のガス圧力Pz1とは原則的に異なる値に計算される。

次いで、ＣＰＵ６１はステップ２２３５に進み、ゾーン２内の消費燃料量sumqrz2を、現時点での値に値「qrz2・Δt」を加えることで更新する。ここで、qrz2としては、先のステップ２１４５（或るいは、後述するステップ２１９０）にて求めた値が使用される。このステップの計算は上記(26)式の計算に相当する。

そして、ＣＰＵ６１はステップ２２４０に進み、ゾーン２の体積の前回値Vbz2を上記ゾーン２の体積Vz2に設定し、続くステップ２２４５にて噴射後経過時間tz2を現時点での値（＝Δt）に微小時間Δtを加えることで更新した後、ステップ２２９５を経由して図１５のステップ１５４０に復帰する。

以上により、クランク角度CAが燃料噴射時期CAinjに達した時点以降において、クランク角度CA（＝CAinj＋j・ΔCA）におけるゾーン２内のガスの各物理量が計算・更新される。この段階では、ゾーン２内における低温酸化反応に基づく予混合燃焼による燃焼速度qrz2(g/sec)（従って、熱発生率Hrz2(j)）が計算されていく。

ＣＰＵ６１はステップ１５４０に進むと、フラグＸＺ３の値が「１」であるか否かを判定する。上述したように、現時点では、先のステップ１５１５の処理によりフラグＸＺ３の値は「０」である。従って、ＣＰＵ６１はステップ１５４０にて「Ｎｏ」と判定してステップ１５４５に進み、クランク角度CA（＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)を上記ステップ２１５０にて求めたゾーン２の熱発生率Hrz2(j)と等しい値に設定し、その値を変数jの値に対応させながらバックアップＲＡＭ６４に記憶する。

続いて、ＣＰＵ６１はステップ１５５０に進み、先のステップ１９４５にて更新されているゾーン１内のガス圧力Pz1と、先のステップ２２３０にて更新されているゾーン２内のガス圧力Pz2とから、「Pz1＝Pz2＝Pg」の関係が成立するように、燃焼室内圧力Pg、ゾーン１の体積Vz1（従って、その前回値Vbz1）、ゾーン２の体積Vz2（従って、その前回値Vbz2）を調整・設定する。この処理は、燃焼室内の圧力は燃焼室内で均一であるとの前提に基づくものである。

次に、ＣＰＵ６１はステップ１５５５に進んで、同ステップ１５５５を経由して図２６〜図２９にフローチャートにより示した「不均一度（σ）の計算」を行う一連のルーチン（サブルーチン）の処理をステップ２６００から開始する。

ＣＰＵ６１は、ステップ２６００からステップ２６０５に進むと、先のステップ２０４５にて更新されている噴霧体積VOLtotalから噴霧体積の前回値VOLbtotalを減じることで噴霧体積の増加量ΔVmを求める。VOLbtotalとしては、今回（初回）に限り先のステップ１８０５の処理により「０」が使用されるが、次回からは、後述するステップ２９５０にて１回前の本ルーチン実行時にて更新されている値が使用される。

次に、ＣＰＵ６１はステップ２６１０に進み、微小時間Δtの間に混合気内に取り込まれた筒内ガス（即ち、上記新規筒内ガス）の質量ΔGmを、先のステップ１９１５にて更新されているゾーン１内のガス密度ρz1に上記噴霧体積の増加量ΔVmを乗じることで求める。

続いて、ＣＰＵ６１はステップ２６１５に進んで、先のステップ２１１０と同様、変数jの値が上記分割数qdiv以下であるか否か（即ち、燃料の噴射継続中であるか否か）を判定する。上述のように、現時点では、変数j＝１であり「j≦qdiv」が成立するから（燃料噴射継続中であるから）、ＣＰＵ６１はステップ２６１５にて「Ｙｅｓ」と判定してステップ２６２０に進み、先のステップ１７３５にて設定されている燃料噴射圧力Pcr(j)から燃焼室内圧力Pgを減じることで有効噴射圧力ΔPmを求める。

次いで、ＣＰＵ６１はステップ２６２５に進んで、上記有効噴射圧力ΔPmと、液体燃料の密度ρfと、上記関数funcVELinjとに基づいて微小時間Δtの間に噴射された燃料（即ち、上記新規燃料）の速度VELinjを求める。

次に、ＣＰＵ６１はステップ２６３０に進み、ガス塊(0)（即ち、筒内ガスの塊）の質量及び速度MAS(0),VEL(0)と、上記新規筒内ガスの質量ΔGmと、新規筒内ガスの速度VELgas（本例では「０」）と、上記(41)式に相当する式とに基づいて筒内ガスの塊の速度VEL(0)を更新する。MAS(0)としては、今回（初回）に限り先のステップ１８１０の処理により「０」が使用されるが、次回からは、後述するステップ２７６５、或いはステップ２８５５にて１回前の本ルーチン実行時にて更新されている値が使用される。VEL(0)としては、今回（初回）に限り先のステップ１８１０の処理により「０」が使用されるが、次回からは、ステップ２６３０、或いは後述するステップ２８４５にて１回前の本ルーチン実行時にて更新されている値が使用される。

続いて、ＣＰＵ６１はステップ２６３５に進み、ガス塊(Nd)（即ち、燃料蒸気の塊）の質量及び速度MAS(Nd),VEL(Nd)と、先のステップ１７４０にて設定されている燃料量q(j)と、先のステップ２６２５にて求めた新規燃料の速度VELinjと、上記(42)式に相当する式とに基づいて燃料蒸気の塊の速度VEL(Nd)を更新する。MAS(Nd)としては、今回（初回）に限り先のステップ１８１０の処理により「０」が使用されるが、次回からは、後述するステップ２７６５、或いはステップ２８５５にて１回前の本ルーチン実行時にて更新されている値が使用される。VEL(Nd)としては、今回（初回）に限り先のステップ１８１０の処理により「０」が使用されるが、次回からは、ステップ２６３５、或いは後述するステップ２８５０にて１回前の本ルーチン実行時にて更新されている値が使用される。

次に、ＣＰＵ６１はステップ２６４０に進んで、上記(36)式に従って、筒内ガスの塊の質量MAS(0)を、現時点での値に上記新規筒内ガスの質量ΔGmを加えることで更新し、続くステップ２６４５にて、上記(37)式に従って、燃料蒸気の塊の質量MAS(Nd)を、現時点での値に上記燃料量q(j)を加えることで更新する。

続いて、ＣＰＵ６１はステップ２６５０に進み、混合気内のガス塊(i)の質量MAS(i)(i＝0,1,・・・,Nd-1,Nd)の総和（ガス塊質量和MAStotal）を求め、続くステップ２６５５にて、上記ガス塊質量和MAStotalを上記噴霧体積VOLtotalで除することで上記(33)式に従って混合気の平均密度ROaveを求める。

次いで、ＣＰＵ６１はステップ２６６０に進んで、ガス塊(i)の質量MAS(i)を上記混合気の平均密度ROaveで除することでガス塊(i)の体積VOL(i)（i＝0,1,・・・,Nd-1,Nd）を求め、続くステップ２６６５にて、上記体積VOL(i)と、ステップ２６６５内に記載の式とからガス塊(i)の直径DIA(i)
（i＝0,1,・・・,Nd-1,Nd）を求める。

次に、ＣＰＵ６１は図２７のステップ２７０５に進み、区分番号m1を「０」に、区分番号m2を「m1+1」に設定する。続いて、ＣＰＵ６１はステップ２７１０に進んで、ガス塊(m1)の体積VOL(m1)もガス塊(m2)の体積VOL(m2)も「０」でないか否かを判定し、「Ｎｏ」と判定する場合、後述する図２９のステップ２９０５に直ちに進む。「ガス塊(m1)とガス塊(m2)についての衝突反応」において少なくとも一方のガス塊の体積（従って、質量）が「０」の場合、衝突反応が発生しない（観念し得ない）。従って、この場合、衝突反応に係わる不必要な計算が省略され得る。

いま、ガス塊(m1)の体積VOL(m1)もガス塊(m2)の体積VOL(m2)も「０」でないものとして説明を続ける。この場合、ＣＰＵ６１はステップ２７１０にて「Ｙｅｓ」と判定してステップ２７１５に進み、体積VOL(m1)が体積VOL(m2)よりも大きいか否かを判定し、「Ｙｅｓ」と判定する場合、ステップ２７２０に進んで区分番号mbを値m1に、区分番号msを値m2に設定する。一方、ステップ２７１５にて「Ｎｏ」と判定される場合、ステップ２７２５に進んで区分番号mbを値m2に、区分番号msを値m1に設定する。これにより、衝突対象ガス塊がガス塊(mb)とガス塊(ms)とに設定されるとともに、ガス塊(mb)が「大ガス塊」となり、ガス塊(ms)が「小ガス塊」となる。

次いで、ＣＰＵ６１はステップ２７３０に進んで、小ガス塊(ms)の直径DIA(ms)と、ガス塊(ms)及びガス塊(mb)の速度VEL(ms)，VEL(mb)と、上記(38)式に相当する式とに基づいて通過体積VOLswpを求める。

続いて、ＣＰＵ６１はステップ２７３５に進み、係数αを、エンジン回転速度NEと、NEを引数とする関数funcαとに基づいて求める。これにより、係数αは、エンジン回転速度NEが大きいほどより小さい値に設定される。

次に、ＣＰＵ６１はステップ２７４０に進んで、上記係数αと、上記通過体積VOLswpと、大ガス塊(mb)の体積VOL(mb)と、上記噴霧体積VOLtotalと、(39)式に相当する式とに基づいて仮の混合体積VOLmix1を求め、続くステップ２７４５にて、混合体積VOLmixを、上記仮の混合体積VOLmix1と小ガス塊(ms)の体積VOL(ms)のうち小さい方の値に設定する。これにより、次のステップ２７５０にてガス塊(ms)及びガス塊(mb)の体積VOL(ms)，VOL(mb)が負となる事態（即ち、衝突反応後において衝突対象ガス塊の体積が負となる事態）の発生が防止される。

続いて、ＣＰＵ６１はステップ２７５０に進んで、小ガス塊(ms)の現時点での体積VOL(ms)から上記混合体積VOLmixを減じることで衝突反応後における小ガス塊(ms)の体積VOL(ms)を求め、大ガス塊(mb)の現時点での体積VOL(mb)から上記混合体積VOLmixを減じることで衝突反応後における大ガス塊(mb)の体積VOL(mb)を求める。

次に、ＣＰＵ６１はステップ２７５５に進み、上記衝突反応後における体積VOL(mb)及び体積VOL(ms)と、ステップ２７５５内に記載の式とから、衝突反応後におけるガス塊(mb)及びガス塊(ms)の直径DIA(mb)及びDIA(ms)をそれぞれ求める。

次いで、ＣＰＵ６１はステップ２７６０に進んで、上記混合体積VOLmixに先のステップ２６５５にて求めた混合気の平均密度ROaveを乗じることで、大ガス塊(mb)の離脱部分の質量dMASmbと、小ガス塊(ms)の離脱部分の質量dMASmsとをそれぞれ求める。

続いて、ＣＰＵ６１はステップ２７６５に進み、上記(35)式の第１、第２式に従って、小ガス塊(ms)の現時点での質量MAS(ms)から上記離脱部分の質量dMASmsを減じることで衝突反応後における小ガス塊(ms)の質量MAS(ms)を求め、大ガス塊(mb)の現時点での質量MAS(mb)から上記離脱部分の質量dMASmbを減じることで衝突反応後における大ガス塊(mb)の質量MAS(mb)を求める。

次いで、ＣＰＵ６１は図２８のステップ２８０５に進んで、混合形成ガス塊の燃料質量分率FRACmixを、先のステップ２７６０にて求めた大ガス塊(mb)及び小ガス塊(ms)のそれぞれの離脱部分の質量dMASmb,dMASmsと、先の１７５５にて設定されている大ガス塊(mb)及び小ガス塊(ms)のそれぞれの燃料質量分率FRACmas(mb),FRACmas(ms)と、上記(34)式に相当する式とに基づいて求める。

次に、ＣＰＵ６１はステップ２８１０に進み、値Ｘを、上記求めた混合形成ガス塊の燃料質量分率FRACmixに値Ndを乗じた値に設定し、続くステップ２８１５にて、値Ｘを挿む２つの整数mdとmu(＝md+1)を特定する。これにより、「ガス塊(m1)とガス塊(m2)についての衝突反応」における混合形成ガス塊がガス塊(md)とガス塊(mu)とに設定される。

続いて、ＣＰＵ６１はステップ２８２０に進んで、値frを、値Ｘから値mdを減じた値に設定する。この値frは、上記それぞれの離脱部分のガス塊(md)とガス塊(mu)への配分割合を決定する際に使用される。即ち、ＣＰＵ６１はステップ２８２５に進み、ガス塊(md)の質量の増加分dMASmd及びガス塊(mu)の質量の増加分dMASmuをステップ２８２５内に記載の式に従って求める。これにより、上記それぞれの離脱部分（の質量和(dMASmb+dMASms)）のうち割合(1-fr)がガス塊(md)に配分され、割合frがガス塊(mu)に配分される。

次に、ＣＰＵ６１はステップ２８３０に進み、上記ガス塊(md)の質量の増加分dMASmd及び上記ガス塊(mu)の質量の増加分dMASmuを上記混合気の平均密度ROaveでそれぞれ除することで、上記ガス塊(md)の体積の増加分dVOLmd及び上記ガス塊(mu)の体積の増加分dVOLmuをそれぞれ求める。

続いて、ＣＰＵ６１はステップ２８３５に進んで、衝突反応後における上記ガス塊(md)の体積VOL(md)を、現時点での値に上記ガス塊(md)の体積の増加分dVOLmdを加えることで求め、衝突反応後における上記ガス塊(mu)の体積VOL(mu)を、現時点での値に上記ガス塊(mu)の体積の増加分dVOLmuを加えることで求める。

次に、ＣＰＵ６１はステップ２８４０に進み、上記衝突反応後における体積VOL(md)及び体積VOL(mu)と、ステップ２８４０内に記載の式とから、衝突反応後におけるガス塊(md)及びガス塊(mu)の直径DIA(md)及びDIA(mu)をそれぞれ求める。

次いで、ＣＰＵ６１はステップ２８４５に進んで、上記(40)式に対応するステップ２８４５内に記載の式に従って、ガス塊(md)の速度VEL(md)を更新する。この式は、「衝突反応の前における混合形成ガス塊(md)の運動量と、衝突反応の前における衝突対象ガス塊(mb)及び(ms)の離脱部分における混合形成ガス塊(md)へ配分される部分のそれぞれの運動量との和が、衝突反応の後における混合形成ガス塊(md)の運動量に等しい」との運動量保存則を利用して得られる。

同様に、ＣＰＵ６１はステップ２８５０に進んで、上記(40)式に対応するステップ２８５０内に記載の式に従って、ガス塊(mu)の速度VEL(mu)を更新する。この式は、「衝突反応の前における混合形成ガス塊(mu)の運動量と、衝突反応の前における衝突対象ガス塊(mb)及び(ms)の離脱部分における混合形成ガス塊(mu)へ配分される部分のそれぞれの運動量との和が、衝突反応の後における混合形成ガス塊(mu)の運動量に等しい」との運動量保存則を利用して得られる。

次いで、ＣＰＵ６１はステップ２８５５に進み、上記(35)式の第３式に従って、ガス塊(md)の現時点での質量MAS(md)に上記ステップ２８２５にて求めたガス塊(md)の質量の増加分dMASmdを加えることで衝突反応後におけるガス塊(md)の質量MAS(md)を求め、ガス塊(mu)の現時点での質量MAS(mu)に上記ステップ２８２５にて求めたガス塊(mu)の質量の増加分dMASmuを加えることで衝突反応後におけるガス塊(mu)の質量MAS(mu)を求める。以上により、「ガス塊(m1)とガス塊(m2)についての衝突反応」により変化する種々の物理量の計算が完了する。

次に、ＣＰＵ６１は図２９のステップ２９０５に進み、区分番号m2が値Ndと等しいか否かを判定する。現時点では、区分番号m2は先のステップ２７０５の処理により「１」となっていて値Ndよりも小さい。従って、ＣＰＵ６１はステップ２９０５にて「Ｎｏ」と判定してステップ２９１０に進み、区分番号m2の値を「１」だけインクリメントした後、図２７のステップ２７１０に戻る。

そして、ＣＰＵ６１はステップ２７１０にて「Ｙｅｓ」と判定する場合、新たな組み合わせに係る「ガス塊(m1)とガス塊(m2)についての衝突反応」について、ステップ２７１５〜ステップ２８５５までの処理を実行する。これにより、新たな組み合わせに係る「ガス塊(m1)とガス塊(m2)についての衝突反応」により変化する種々の物理量の計算が完了する。このような処理は、ステップ２９１０の繰り返し実行により増大していく区分番号m2が値Ndに達するまで繰り返し実行される。

この状態にて区分番号m2が値Ndに達したものとすると、ＣＰＵ６１はステップ２９０５にて「Ｙｅｓ」と判定してステップ２９１５に進み、区分番号m1が値「Nd-1」と等しいか否かを判定する。現時点では、区分番号m1は先のステップ２７０５の処理により「０」となっていて値「Nd-1」よりも小さい。従って、ＣＰＵ６１はステップ２９１５にて「Ｎｏ」と判定してステップ２９２０に進み、区分番号m1の値を「１」だけインクリメントし、続くステップ２９２５にて区分番号m2を値「m1+1」に設定した後、図２７のステップ２７１０に戻る。

このような処理は、ステップ２９２０の繰り返し実行により増大していく区分番号m1が値「Nd-1」に達するまで繰り返し実行される。この結果、混合気内において現時点にて存在している（体積VOL(i)が「０」でない）複数のガス塊(i)のうちの２つの組み合わせの全てについて上記衝突反応に係わる計算が順に実行されていく。

そして、この状態にて区分番号m1が値「Nd-1」に達したものとすると、ＣＰＵ６１はステップ２９１５にて「Ｙｅｓ」と判定してステップ２９３０に進み、混合気内における燃料質量分率の平均値FRACmasaveを、現時点でのガス塊(i)の質量MAS(i)(i＝0,1,・・・,Nd)と、先のステップ２６５０にて求めたガス塊質量和MAStotalと、上記(44)式とに基づいて求める。

次いで、ＣＰＵ６１はステップ２９３５に進んで、ガス塊(i)の質量MAS(i)の上記ガス塊質量和MAStotalに対する割合（質量分率pdf(i)）を求め(i＝0,1,・・・,Nd)、続くステップ２９４０にて、上記燃料質量分率の平均値FRACmasaveと、上記質量分率pdf(i)(i＝0,1,・・・,Nd)と、ステップ２９４０内に記載の式とに基づいて燃料質量分率の分散xrmsを求める。

そして、ＣＰＵ６１はステップ２９４５に進み、上記分散xrmsの平方根をとることで燃料質量分率の標準偏差σを求めた後、ステップ２９５０に進んで噴霧体積の前回値VOLbtotalを先のステップ２０４５にて更新されている噴霧体積VOLtotalと等しい値に設定した後、ステップ２９９５を経由して図１６のステップ１６０５に進む。

ここで、ステップ２９３５、２９４０、２９４５の一連の計算は上記(43)式の計算に対応している。これにより、現時点での混合気内における燃料質量分率の標準偏差σが計算される。この標準偏差σは、後述するゾーン２の着火後の燃焼速度qrz2の計算（図２１のステップ２１８０、２１８５）、及び後述するゾーン３の燃焼速度qrz3の計算（図２４のステップ２４４０）に反映されていく。

ＣＰＵ６１は図１６のステップ１６０５に進むと、フラグＸＩＧの値が「０」であるか否かを判定する。現時点では、先のステップ１５１５の処理によりフラグＸＩＧの値は「０」となっている。従って、ＣＰＵ６１はステップ１６０５にて「Ｙｅｓ」と判定してステップ１６１０に進み、先のステップ２２２５にて更新・増大されているゾーン２内のガス温度Tz2が所定の着火温度Tigより大きいか否か（即ち、着火したか否か）を判定する。

現時点では、クランク角度CAは燃料噴射時期CAinjの直後であるから、ゾーン２内のガス温度Tz2は着火温度Tigより十分に小さい。従って、ＣＰＵ６１はステップ１６１０にて「Ｎｏ」と判定してステップ１６１５に進み、変数jの値を「１」だけインクリメントした後、ステップ１６２０〜１６２５の処理を行う。これにより、変数jのインクリメントに合わせて計算上のクランク角度CAが微小クランク角度ΔCAだけ進行していく。

そして、ＣＰＵ６１は、ステップ１６２５にて「Ｎｏ」と判定する場合、図１５のステップ１５２０に戻り、ステップ１５２０（図１９のルーチン）、１５２５、１５３５（図２０〜図２２のルーチン）、１５４０、１５４５、１５５０、１５５５（図２６〜図２９のルーチン）、図１６のステップ１６０５、１６１０、１６１５、１６２０、１６２５の処理を再び順に実行していく。

これにより、クランク角度CA（＝CAinj＋j・ΔCA）の進行に合わせて、ゾーン１内のガスの各物理量、ゾーン２内のガスの各物理量、及び混合気内における燃料質量分率の標準偏差σが計算・更新されていく。また、この段階では、クランク角度CA（＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)は、ゾーン２内における低温酸化反応に基づく予混合燃焼の燃焼速度qrz2(g/sec)から得られる熱発生率Hrz2(j)と等しい値に計算されていく。

このような処理は、図２１のステップ２１４０にて「Ｙｅｓ」と判定される限りにおいて（従って、ゾーン２内のガス温度Tz2が着火温度Tig以下である限りにおいて）繰り返し実行されていく。

以下、先ずは、ゾーン２内のガス温度Tz2が着火温度Tig以下である状態で燃料噴射が終了する場合（即ち、上記着火後噴射継続でない場合。ゾーン３が発生しない場合）について説明する。この場合、ＣＰＵ６１は図２１のステップ２１１０に進んだとき「Ｎｏ」と判定してステップ２１５５に進むようになり、以降、（図２２のステップ２２１０にて使用される）潜熱Qlatentz2を「０」に設定する。加えて、ステップ２１１５の処理が実行されないから、以降、（ステップ２１３０にて使用される）ゾーン２内噴射燃料量Qz2は現時点での値（＝q(1)＋q(2)＋・・・＋q(qdiv)＝Qfin）に固定される。

また、ＣＰＵ６１は図２６のステップ２６１５に進んだとき「Ｎｏ」と判定してステップ２６７０に進むようになり、燃料量q(j)の値を「０」に設定し、続くステップ２６７５にて新規燃料の速度VELinjを「０」に設定するようになる。従って、以降、ステップ２６３５の処理を行っても、ガス塊(Nd)（従って、燃料蒸気の塊）の速度VEL(Nd)は変化しない。

この状態にて、ゾーン２内のガス温度Tz2が着火温度Tigを超えると（即ち、着火すると）、ＣＰＵ６１は図１６のステップ１６１０に進んだとき「Ｙｅｓ」と判定してステップ１６３０に進み、フラグＸＩＧの値を「０」から「１」に変更し、続くステップ１６３５にて変数jの値が分割数qdivより小さいか否か（即ち、燃料噴射継続中か否か）を判定する。

現時点では、上述のごとく燃料噴射が終了している。従って、ＣＰＵ６１はステップ１６３５にて「Ｎｏ」と判定してステップ１６４０を実行することなくステップ１６１５に進む。これにより、以降、ＣＰＵ６１は図１６のステップ１６０５に進んだとき「Ｎｏ」と判定してステップ１６１５に直ちに進むようになる。この結果、フラグＸＩＧの値は「１」に維持される一方で、フラグＸＺ３の値は「０」に維持される。

従って、以降も、ＣＰＵ６１は図１５のステップ１５４０に進む毎に「Ｎｏ」と判定し続ける。この結果、ステップ１５６０（図２３〜図２５のルーチン）の「ゾーン３の計算」が実行されない。

一方、これ以降、ＣＰＵ６１は図２１のステップ２１４０に進んだとき「Ｎｏ」と判定してステップ２１６０に進むようになり、層流特性時間τaz2を、先のステップ２１３５にて求めたゾーン２内の酸素濃度[O2]z2と、先のステップ２１３０にて求めたゾーン２内の燃料濃度[Fuel]z2と、前回の本ルーチン実行時において図２２のステップ２２２５にて更新されているゾーン２内のガス温度Tz2と、上記(22)式とに基づいて求める。

続いて、ＣＰＵ６１はステップ２１６５に進み、乱流特性時間τmz2を、図２０のステップ２０３５にて求めた値(k/ε)と、上記ゾーン２内の酸素濃度[O2]z2と、上記ゾーン２内の燃料濃度[Fuel]z2と、上記関数funcτmと、上記(23)式とに基づいて求める。

次いで、ＣＰＵ６１はステップ２１７０に進んで、係数C3を、前回の本ルーチン実行時において図２２のステップ２２３５にて更新されているゾーン２の消費燃料量sumqrz2と、図１７のステップ１７１５にて計算されている燃料噴射量Qfinと、関数funcC3とに基づいて求める。これにより、係数C3は値(sumqrz2/Qfin)が大きくなるほど大きい値に設定される。

次に、ＣＰＵ６１はステップ２１７５に進み、上記求めた層流特性時間τaz2と、乱流特性時間τmz2と、係数C3と、上記(21)式とに基づいてゾーン２内の燃焼に係わる特性時間τcz2を求める。次いで、ＣＰＵ６１はステップ２１８０に進んで、補正係数Ａを、図２９のステップ２９４５にて更新されている標準偏差σと、図１３に示したテーブルとに基づいて求める。

続いて、ＣＰＵ６１はステップ２１８５に進み、燃料濃度減少速度d[Fuel]z2を、上記ゾーン２内の燃料濃度[Fuel]z2と、上記特性時間τcz2と、上記補正係数Ａと、ステップ２１８５内に記載の式とに基づいて求め、続くステップ２１９０にてゾーン２の燃焼速度qrz2(g/sec)を、上記求めた燃料濃度減少速度d[Fuel]z2と、ステップ２１２５にて求めたゾーン２内のガス質量Mz2と、ステップ２１９０内に記載の式とに基づいて求める。このステップ２１８５とステップ２１９０の計算は、上記(45)式の計算に相当する。

このように、上記着火後噴射継続でない場合（即ち、ゾーン３が発生しない場合）における着火後（即ち、フラグＸＩＧ＝１、フラグＸＺ３＝０）は、クランク角度CA（＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)は、ゾーン２内における高温酸化反応に基づく予混合燃焼と拡散燃焼とを考慮した燃焼速度qrz2(g/sec)から得られる熱発生率Hrz2(j)と等しい値に計算されていく（図２１のステップ２１５０、図１５のステップ１５４５を参照）。

そして、図１６のステップ１６２０の繰り返し実行により微小クランク角度ΔCA毎に進行していく計算上のクランク角度CAが上記計算終了クランク角度CAendに達すると、ＣＰＵ６１はステップ１６２５に進んだとき「Ｙｅｓ」と判定して図１５〜図１６の一連のメインルーチンを終了する。

上述したように、メインルーチンに係る処理は、ＩＶＣの直後に完了する。従って、この場合、上記着火後噴射継続でない場合について、今回燃料が噴射される気筒についてのクランク角度CAinj〜クランク角度CAendの間におけるクランク角度CA（CA＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)(j＝1,2,・・・)がＩＶＣの直後（従って、実際の燃料噴射時期CAinjの到来前）に取得され得る。

従って、今回の燃料噴射による燃焼に基づいて発生するであろう燃焼騒音や発生トルク等がＩＶＣの直後に予測され得るから、その予測結果に基づいて燃焼騒音や発生トルク等が適切な値に維持されるように、今回の燃料噴射についての燃料噴射形態（燃料噴射圧力、燃料噴射時期等）をフィードバック制御することができる。

次に、上記着火後噴射継続の場合（ゾーン３が発生する場合）について説明する。この場合、燃料噴射継続中（即ち、変数j＜qdivの間）にゾーン２内のガス温度Tz2が着火温度Tigを超える。従って、ゾーン２内のガス温度Tz2が着火温度Tigを超えると、ＣＰＵ６１は図１６のステップ１６１０に進んだとき、「Ｙｅｓ」と判定してステップ１６３０の処理を実行し、続くステップ１６３５でも「Ｙｅｓ」と判定してステップ１６４０に進み、フラグＸＺ３の値を「０」から「１」に変更する。

この結果、以降、フラグＸＩＧの値に加えてフラグＸＺ３の値も「１」に維持される。従って、ＣＰＵ６１は図１５のステップ１５４０に進んだとき「Ｙｅｓ」と判定してステップ１５６０（従って、図２３〜図２５の「ゾーン３の計算」ルーチン）を実行するようになる。即ち、ステップ１５２０（図１９のルーチン）、１５２５、１５３５（図２０〜図２２のルーチン）、１５４０、１５６０（図２３〜図２５のルーチン）、１５４５、１５５０、１５５５（図２６〜図２９のルーチン）、図１６のステップ１６０５、１６１５、１６２０、１６２５の処理が繰り返し実行されていく。図２３〜図２５の「ゾーン３の計算」ルーチンの詳細については後述する。

これにより、クランク角度CA（＝CAinj＋j・ΔCA）の進行に合わせて、ゾーン１内のガスの各物理量、ゾーン２内のガスの各物理量、ゾーン３内のガスの各物理量、及び混合気内における燃料質量分率の標準偏差σが計算・更新されていく。また、この場合、図２１のステップ２１４０にて「Ｎｏ」と判定されて上述したステップ２１６０〜２１９０が実行されるから、クランク角度CA（＝CAinj＋j・ΔCA）に対するゾーン２の熱発生率Hrz2(j)は、ゾーン２内における高温酸化反応に基づく予混合燃焼と拡散燃焼とを考慮した燃焼速度qrz2(g/sec)から計算されていく（ステップ２１５０を参照）。

また、後述するように、クランク角度CA（＝CAinj＋j・ΔCA）に対するゾーン３の熱発生率Hrz3(j)は、ゾーン３内における拡散燃焼のみを考慮した燃焼速度qrz3(g/sec)から計算されていく（後述する図２４のステップ２４５０を参照）。クランク角度CA（＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)は、ゾーン２の熱発生率Hrz2(j)とゾーン３の熱発生率Hrz3(j)の和（Hrz2(j)＋Hrz3(j)）と等しい値に計算されていく（図１５のステップ１５４５を参照）。

なお、この場合、着火時点以降（即ち、フラグＸＺ３の値が「１」になった時点以降）、ＣＰＵ６１は図２１のステップ２１０５に進んだとき「Ｎｏ」と判定してステップ２１５５に直ちに進む。即ち、着火時点以降、ゾーン２に対する潜熱Qlatentz2が「０」に維持されるとともに、ゾーン２内噴射燃料量Qz2が現時点での値（＜Qfin）に固定される。これは、着火時点以降、燃料噴射はゾーン３に対してのみ実行されゾーン２に対しては実行されなくなることに対応した処理である。

加えて、この場合、着火時点以降（即ち、フラグＸＺ３の値が「１」になった時点以降）、ＣＰＵ６１は図２０のステップ２０５０に進んだとき「Ｎｏ」と判定してステップ２０７５に進んで、上記(28)式に従って、ステップ２０４５にて求めた噴霧体積VOLtotalから後述する図２３のステップ２３２５にて更新されているゾーン３の体積Vz3を減じることでゾーン２の体積Vz2を求めるようになる。これは、着火時点以降、噴霧体積VOLtotal中においてゾーン３が占める部分が発生し始めることに対応している。

次に、図２３〜図２５の「ゾーン３の計算」ルーチンの詳細に説明する。上述したように、上記着火後噴射継続の場合、着火時点においてフラグＸＺ３の値が「０」から「１」に変更されるから、着火時点以降、ＣＰＵ６１は図１５のステップ１５４０に進んだとき「Ｙｅｓ」と判定してステップ１５６０に進むようになり、同ステップ１５６０を経由して図２３〜図２５にフローチャートにより示した「ゾーン３の計算」を行う一連のルーチン（サブルーチン）の処理をステップ２３００から開始する。

ＣＰＵ６１は、ステップ２３００からステップ２３０５に進むと、フラグＸＺ３の値が「０」から「１」に変更された直後であるか否かを判定する。いま、図１６のステップ１６４０の処理によりフラグＸＺ３の値が「０」から「１」に変更された直後であるものとすると、ＣＰＵ６１はステップ２３０５にて「Ｙｅｓ」と判定してステップ２３１０に進み、ゾーン３内のガス温度Tz3（の初期値）を先のステップ２２２５にて更新されているゾーン２内のガス温度Tz2と等しい値（即ち、着火開始時点でのゾーン２内のガス温度）に設定する。

続いて、ＣＰＵ６１はステップ２３１５に進んで、着火後経過時間tz3を微小時間Δtに設定する。以上、ステップ２３１０及びステップ２３１５にてゾーン３の計算を開始するにあたり必要な各種値が計算される。

次に、ＣＰＵ６１はステップ２３２０に進み、上記有効噴射圧力ΔP0と、着火後経過時間tz3と、上記(12)式に対応する関数funcとに基づいて到達距離Sz3を求め、続くステップ２３２５にて、先のステップ２０２５にて求めた噴霧角θと、到達距離Sz3と、上記(13)式に相当する式とに基づいてゾーン３の体積Vz3を求める。

次いで、ＣＰＵ６１はステップ２３３０に進んで、上記求めたゾーン３の体積Vz3からゾーン３の体積の前回値Vbz3を減じることで、微小時間Δtの間におけるゾーン３の体積の増加量ΔVz3を求める。Vbz3としては、今回（初回）に限り先のステップ１８０５の処理により「０」が使用されるが、次回からは、後述するステップ２５４０にて１回前の本ルーチン実行時にて更新されている値が使用される。

次に、ＣＰＵ６１はステップ２３３５に進み、先のステップ１９１５にて更新されているゾーン１内のガス密度ρz1に上記ゾーン３の体積の増加量ΔVz3を乗じることで、微小時間Δtの間におけるゾーン３内筒内ガス質量の増加量ΔGz3を求める。続いて、ＣＰＵ６１はステップ２３４０に進んで、ゾーン３内筒内ガス質量Gz3を、現時点での値（今回（初回）は、先のステップ１８１５の処理により「０」）に上記増加量ΔGz3を加えることで更新する。

続いて、ＣＰＵ６１は図２４のステップ２４０５に進み、変数jの値が分割数qdiv以下か否かを判定する。上述と同様、「j≦qdiv」は、燃料の噴射継続中に対応し、「j＞qdiv」は、燃料の噴射終了後に対応する。

燃料噴射継続中の場合（j≦qdiv）、ＣＰＵ６１はステップ２４０５にて「Ｙｅｓ」と判定してステップ２４１０に進み、ゾーン３内噴射燃料量Qz3を、現時点での値（今回（初回）は、先のステップ１８２０の処理により「０」）に先のステップ１７４０にて求めた微小時間Δtの間の燃料量q(j)を加えることで更新する。次いで、ＣＰＵ６１はステップ２４１５に進んで、潜熱Qlatentz3を、上記燃料量q(j)と、q(j)を引数とする潜熱を求める関数funcQlatentとに基づいて求める。

一方、燃料噴射終了後の場合（j＞qdiv）、ＣＰＵ６１はステップ２４０５にて「Ｎｏ」と判定してステップ２４５５に進み、以降、潜熱Qlatentz3を「０」にする。

ＣＰＵ６１はステップ２４２０に進むと、ゾーン３内のガス質量Mz3を、先のステップ２３４０にて求めたゾーン３内筒内ガス質量Gz3に先のステップ２４１０にて求めたゾーン３内噴射燃料量Qz3を加えることで求める。

続いて、ＣＰＵ６１はステップ２４２５に進み、ゾーン３内の燃料濃度[Fuel]z3を、上記ゾーン３内噴射燃料量Qz3と、ゾーン３内の消費燃料量sumqrz3と、上記ゾーン３内のガス質量Mz3と、上記(25)式に相当する式とに基づいて求める。sumqrz3としては、今回（初回）に限り先のステップ１８２５の処理により「０」が使用されるが、次回からは、後述するステップ２５３５にて１回前の本ルーチン実行時にて更新されている値が使用される。

次いで、ＣＰＵ６１はステップ２４３０に進んで、ゾーン３内の酸素濃度[O2]z3を、上記ゾーン３内筒内ガス質量Gz3と、先のステップ１７０５にて設定されている吸気酸素濃度[O2]inと、上記消費燃料量sumqrz3と、上記ゾーン３内のガス質量Mz3と、上記(27)式に相当する式とに基づいて求める。sumqrz3としては、ステップ２４２５にて使用される値と同じ値が使用される。

次に、ＣＰＵ６１はステップ２４３５に進み、乱流特性時間τmz3を、図２０のステップ２０３５にて求めた値(k/ε)と、上記ゾーン３内の酸素濃度[O2]z3と、上記ゾーン３内の燃料濃度[Fuel]z3と、上記関数funcτmと、上記(32)式とに基づいて求める。

次いで、ＣＰＵ６１はステップ２４４０に進んで、燃料濃度減少速度d[Fuel]z3を、上記ゾーン３内の燃料濃度[Fuel]z23と、上記乱流特性時間τmz3と、上記補正係数Ａと、ステップ２４４０内に記載の式とに基づいて求め、続くステップ２４４５にてゾーン３の燃焼速度qrz3(g/sec)を、上記求めた燃料濃度減少速度d[Fuel]z3と、ステップ２４２０にて求めたゾーン３内のガス質量Mz3と、ステップ２４４５内に記載の式とに基づいて求める。このステップ２４４０とステップ２４４５の計算は、上記(46)式の計算に相当する。

続いて、ＣＰＵ６１はステップ２４５０に進み、ゾーン３の熱発生率Hrz3(j)を、上記求めた燃焼速度qrz3と、上記(24)式に相当する式とに基づいて求める。この熱発生率Hrz3(j)は、クランク角度CA（＝CAinj＋j・ΔCA）における値である。このように、ゾーン３の熱発生率Hrz3(j)は、ゾーン３内における拡散燃焼のみを考慮した燃焼速度qrz3(g/sec)から計算されていく。

次に、ＣＰＵ６１は図２５のステップ２５０５に進み、ゾーン３に対する微小ピストン仕事ΔWpistonz3を、燃焼室内圧力Pg（図１５のステップ１５５０を参照）と、上記ゾーン３の体積Vz3と、上記(19)式に相当する式とに基づいて求める。

次いで、ＣＰＵ６１はステップ２５１０に進んで、ゾーン３に対する微小エネルギーΔWz3を上記ゾーン３に対する微小ピストン仕事ΔWpistonz3と、先のステップ２４１５にて求めた潜熱Qlatentz3と、先のステップ２４５０にて求めた熱発生率Hrz3(j)と、上記(18)式に相当する式とに基づいて求める。

次に、ＣＰＵ６１はステップ２５１５に進み、ゾーン３内のガスの定圧比熱Cpz3を、上記吸気酸素濃度[O2]inと、上記排気酸素濃度[O2]exbと、ゾーン３内のガス温度Tz3と、先のステップ２４２５にて求めたゾーン３内の燃料濃度[Fuel]z3と、上記(17)式に相当する式とに基づいて求める。

続いて、ＣＰＵ６１はステップ２５２０に進んで、ゾーン３内のガス温度の増大量ΔTz3を、先のステップ２５１０にて得られた微小エネルギーΔWz3と、先のステップ２４２０にて求めたゾーン３内のガス質量Mz3と、上記ゾーン３内のガスの定圧比熱Cpz3とに基づいて求める。次いで、ＣＰＵ６１はステップ２５２５に進み、ゾーン３内のガス温度Tz3を、現時点での値に上記増大量ΔTz3を加えることで更新する。

次に、ＣＰＵ６１はステップ２５３０に進んで、ゾーン３内のガス圧力Pz3を、上記ゾーン３内のガス質量Mz3と、上記更新されたゾーン３内のガス温度Tz3と、先のステップ２３２５にて求めたゾーン２の体積Vz3と、ゾーン３内のガスについての気体の状態方程式（ステップ２５３０内の式を参照）とに基づいて求める。このゾーン３内のガス圧力Pz3は先のステップ１９４５にて更新されているゾーン１内のガス圧力Pz1とも、先のステップ２２３０にて更新されているゾーン２内のガス圧力Pz2とも原則的に異なる値に計算される。

次いで、ＣＰＵ６１はステップ２５３５に進み、ゾーン３内の消費燃料量sumqrz3を、現時点での値に値「qrz3・Δt」を加えることで更新する。ここで、qrz3としては、先のステップ２４４５にて求めた値が使用される。

そして、ＣＰＵ６１はステップ２５４０に進み、ゾーン３の体積の前回値Vbz3を上記ゾーン３の体積Vz3に設定し、続くステップ２５４５にて着火後経過時間tz3を現時点での値に微小時間Δtを加えることで更新した後、ステップ２５９５を経由して図１５のステップ１５４５に復帰する。

以上により、上記着火後噴射継続の場合、着火後において、クランク角度CA（＝CAinj＋j・ΔCA）におけるゾーン３内のガスの各物理量が計算・更新されていく。上記着火後噴射継続の場合も、上述した上記着火後噴射継続でない場合と同様、今回燃料が噴射される気筒についてのクランク角度CAinj〜クランク角度CAendの間におけるクランク角度CA（CA＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)(j＝1,2,・・・)がＩＶＣの直後（従って、実際の燃料噴射時期CAinjの到来前）に取得され得る。

従って、今回の燃料噴射による燃焼に基づいて発生するであろう燃焼騒音や発生トルク等がＩＶＣの直後に予測され得るから、その予測結果に基づいて燃焼騒音や発生トルク等が適切な値に維持されるように、今回の燃料噴射についての燃料噴射形態（燃料噴射圧力、燃料噴射時期等）をフィードバック制御することができる。

なお、前記不均一度取得手段は図２６〜図２９のルーチンに対応し、前記燃料質量割合特定手段は図２８のステップ２８０５に対応し、前記ガス塊量取得手段は図２６のステップ２６４０，２６４５、図２７のステップ２７６５、図２８のステップ２８５５に対応し、前記混合体積取得手段は図２７のステップ２７４５に対応し、前記ガス塊速度取得手段は図２６のステップ２６３０，２６３５、図２８のステップ２８４５，２８５０に対応し、前記混合気状態取得手段は図１９〜図２５のルーチンに対応し、前記混合気温度取得手段は図２２のステップ２２２５、図２５のステップ２５２５に対応し、前記燃料反応速度取得手段は図２１のステップ２１４５，２１９０、図２４のステップ２４４５に対応している。

以上、説明したように、本発明による不均一度取得装置を含んだ混合気状態取得装置の実施形態では、混合気は、燃料質量分率FRACmas(i)が異なる球形のガス塊(i)の集合体であると仮定される。混合気が形成されていく燃料噴射開始時期CAinj以降、微小時間Δtの経過毎（微小クランク角度ΔCAの進行毎）に、その時点で存在しているガス塊(i)のうちの任意の２つの組み合わせの全てについて上述した「衝突反応」が順に行われる。これにより、混合気内のガス塊(i)の質量MAS(i)の分布が微小時間Δt毎に取得・更新され、混合気の不均一度（燃料質量分率の標準偏差σ）が複雑な微分方程式等を解くことなく取得・更新されていく。

この実施形態では、燃料噴射開始時期CAinj以降におけるクランク角度CA（CA＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)(j＝1,2,・・・)が計算される。この熱発生率Hr(j)を計算する際、上記問題点２の解消のため（図６、図１４のＢを参照）、上述のように得られる混合気の不均一度（上記標準偏差σ）が反映される。

加えて、混合気着火後もなお燃料噴射が継続する場合（着火後噴射継続の場合）における上記問題点１の解消のため（図６、図１４のＡを参照）、混合気が、着火時点以前において噴射された燃料に基づく部分（ゾーン２）と着火時点以降において噴射された燃料に基づく部分（ゾーン３）とに分けて取り扱われる。具体的には、ゾーン３の熱発生率Hrz3(j)が、ゾーン２の熱発生率Hrz2(j)とは異なる計算手法により同ゾーン２の熱発生率Hrz2(j)と並行して別個独立して逐次計算されていく（Hr(j)＝Hrz2(j)＋Hrz3(j)）。

以上の手法を採用すると、計算されたクランク角度CA（CA＝CAinj＋j・ΔCA）に対する熱発生率Hr(j)(j＝1,2,・・・)は、実験結果と非常によく一致することが判明した（図１４を参照）。従って、燃焼騒音や発生トルク等が精度良く予測できるから、その予測結果に基づいて燃焼騒音や発生トルク等が適切な値に維持されるように、燃料噴射形態（燃料噴射圧力、燃料噴射時期等）を精度良くフィードバック制御することができる。

本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、上記実施形態では、混合気の不均一度（具体的には、燃料質量分率の標準偏差σ）を熱発生率Hr(j)の計算に反映させるため、標準偏差σに応じた補正係数Ａ（図１３を参照）が考慮された上記(45)式に従ってゾーン２の燃焼速度qrz2が計算されるようになっているが、ゾーン２の燃焼速度qrz2を、上記補正係数Ａが考慮された下記(47)式により得られる特性時間τcz2を用いて、上記(20)式に従って計算してもよい。

また、上記実施形態においては、混合体積VOLmixを求める際の基礎となる通過体積VOLswpを図１０に示したカプセル状の体積を求める式である上記(38)式に従って求めているが、上記カプセル状の体積から両端の半球部分を取り除いた体積（円柱状の体積）を求める式である下記(48)式（上記(38)式において右辺第２項を省略した式）に従って通過体積VOLswpを求めても良い。

また、上記実施形態においては、微小クランク角度ΔCAの進行の間に噴射された燃料の量q(j)（図１７のステップ１７４０を参照）の全てが直ちに燃料蒸気となって燃料蒸気の塊（ガス塊(Nd)）に加わるようになっているが、噴射された燃料の量q(j)の一部が燃料蒸気となるものとしてもよい。この場合、例えば、燃料の量q(j)に係数C（０＜C＜１）を乗じた値が直ちに燃料蒸気となって燃料蒸気の塊（ガス塊(Nd)）に加わるように構成してもよい。

また、上記実施形態においては、混合気中にて発生する上記「衝突反応」の対象となる２つのガス塊の組み合わせを、燃料質量分率FRACmasの小さいものから順に選択しているが、例えば、ガス塊の質量の大きいもの（或いは、小さいもの）から順に選択してもよい。

また、上記実施形態においては、「衝突反応」の対象となる２つのガス塊の離脱部分の体積が共に上記混合体積VOLmix（上記(39)式を参照）と等しいものとして計算されているが、「衝突反応」の対象となる２つのガス塊の離脱部分の体積を異ならせてもよい。

また、上記実施形態においては、混合気の不均一度を表す燃料質量分率の標準偏差σを、混合気内のガス塊(i)の質量MAS(i)の分布から求めているが、混合気内のガス塊(i)の体積VOL(i)の分布から求めてもよい。

また、上記実施形態においては、混合気の不均一度を表す指標値として、燃料質量分率の標準偏差σを使用しているが、燃料質量分率の分散（＝σ²）を使用してもよい。

また、上記実施形態においては、上記(34)式に従って計算される混合形成ガス塊の燃料質量分率FRACmixの値（ステップ２８０５を参照）と完全に一致する燃料質量分率FRACmas(i)を有するガス塊(i)が存在しない場合、上記混合形成ガス塊の燃料質量分率FRACmixの値を挟む２つの燃料質量分率FRACmas(md)，FRACmas(mu)を有する２つのガス塊(md)，ガス塊(mu)が混合形成ガス塊として使用されるが、この場合、上記混合形成ガス塊の燃料質量分率FRACmixの値に最も近い燃料質量分率FRACmas(mz)を有する１つのガス塊(mz)を混合形成ガス塊として使用してもよい。

また、上記実施形態においては、通過体積VOLswpを求める際に使用する上記(38)式において、２つの衝突対象ガス塊の相対速度の値として、単純にそれぞれのガス塊の速度の差の値を使用している（即ち、全てのガス塊の運動方向が同じであると扱っている）が、ガス塊の運動の方向をも考慮して２つの衝突対象ガス塊の相対速度の値を求めてもよい。

また、上記実施形態においては、混合気は、燃料質量分率FRACmas(i)が異なる球形のガス塊(i)の集合体であると仮定されているが、球形とは異なる形状（例えば、立方体形状等）のガス塊(i)の集合体であると仮定してもよい。

加えて、上記実施形態においては、燃焼室内のガスの温度（ゾーン１，２，３内のガス温度Tz1,Tz2,Tz3）を求める際に燃焼室内のガスから燃焼室壁への熱損失が考慮されていないが、燃焼室内のガスから燃焼室壁への熱損失が考慮されてもよい。

本発明の実施形態に係る内燃機関の混合気の不均一度取得装置を含んだ混合気状態取得装置を４気筒内燃機関（ディーゼル機関）に適用したシステム全体の概略構成図である。 或る一つの気筒のシリンダ内(筒内)に吸気マニホールドからガスが吸入され、筒内に吸入されたガスが排気マニホールドへ排出される様子を模式的に示した図である。 吸気弁閉弁後、燃料噴射開始前において筒内ガスのみが存在する燃焼室内の様子を示した図である。 燃料噴射開始後の或る時点において燃焼室内にて燃料蒸気（ゾーン２）が広がっている様子を模式的に示した図である。 噴霧体積の模式的に示した図である。 混合気着火後においても燃料噴射がなお継続する場合におけるクランク角度に対する熱発生率について、実験結果（破線）と従来の装置による計算結果（実線）との比較結果の一例を示したグラフである。 混合気が、着火時点以前において噴射された燃料に基づく部分（ゾーン２）と着火時点以降において噴射された燃料に基づく部分（ゾーン３）とに分けて取り扱われた場合における、着火後の或る時点において燃焼室内にてゾーン２とゾーン３とが広がっている様子を模式的に示した図である。 燃料噴射開始から微小時間Δtが経過した時点（噴射後経過時間t＝Δt）で発生する「衝突反応」の前後の様子を模式的に示した図である。 図８に示した状態から更に微小時間Δtが経過した時点（噴射後経過時間t＝２Δt）で発生する「衝突反応」の前後の様子を模式的に示した図である。 カプセル状の通過体積を模式的に示した図である。 図１に示した装置により燃料噴射開始からの微小時間Δtの経過毎に「衝突反応」に係わる計算を繰り返した場合における混合気内のガス塊の質量分布の一例を示したグラフである。 ガス塊の燃料質量分率を離散的に設定した場合であって、混合形成ガス塊として２つのガス塊が使用される場合の一例を示した図である。 図１に示したＣＰＵが参照する、ガス塊の燃料質量分率の標準偏差と、熱発生率を補正するための補正係数との関係を規定したテーブルを示したグラフである。 混合気着火後においても燃料噴射がなお継続する場合におけるクランク角度に対する熱発生率について、実験結果（破線）と図１に示した装置による計算結果（実線）との比較結果の一例を示したグラフである。 図１に示したＣＰＵが実行するメインルーチンの前半部を示したフローチャートである。 図１に示したＣＰＵが実行するメインルーチンの後半部を示したフローチャートである。 図１に示したＣＰＵが実行する初期値の設定を行うためのサブルーチンの前半部を示したフローチャートである。 図１に示したＣＰＵが実行する初期値の設定を行うためのサブルーチンの後半部を示したフローチャートである。 図１に示したＣＰＵが実行するゾーン１の計算を行うためのサブルーチンを示したフローチャートである。 図１に示したＣＰＵが実行するゾーン２の計算を行うためのサブルーチンの第１番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行するゾーン２の計算を行うためのサブルーチンの第２番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行するゾーン２の計算を行うためのサブルーチンの第３番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行するゾーン３の計算を行うためのサブルーチンの第１番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行するゾーン３の計算を行うためのサブルーチンの第２番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行するゾーン３の計算を行うためのサブルーチンの第３番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する不均一度の計算を行うためのサブルーチンの第１番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する不均一度の計算を行うためのサブルーチンの第２番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する不均一度の計算を行うためのサブルーチンの第３番目の部分を示したフローチャートである。 図１に示したＣＰＵが実行する不均一度の計算を行うためのサブルーチンの第４番目の部分を示したフローチャートである。 図１に示したＣＰＵが図１７に示したルーチンを実行する際に参照する燃料噴射量を決定するためのテーブルである。 図１に示したＣＰＵが図１７に示したルーチンを実行する際に参照する燃料噴射時期を決定するためのテーブルである。 図１に示したＣＰＵが図１７に示したルーチンを実行する際に参照する基本燃料噴射圧力を決定するためのテーブルである。

符号の説明

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

Claims (16)

1. 内燃機関の燃焼室内にて同燃焼室内に吸入されている筒内ガスと噴射された燃料とが混ざり合って混合気が形成されていく過程において、燃料質量割合が異なる２つのガス塊の衝突により前記２つのガス塊のそれぞれの一部分が混合し、同混合した前記各一部分が対応する前記ガス塊からそれぞれ離脱して前記２つのガス塊の燃料質量割合とは異なる燃料質量割合を有する他のガス塊の一部又は全部となる衝突反応、が前記混合気内で繰り返し行われると仮定して、前記混合気内における燃料濃度の不均一度を取得する不均一度取得手段を備えた内燃機関の混合気の不均一度取得装置。
2. 請求項１に記載の内燃機関の混合気の不均一度取得装置において、
   前記不均一度取得手段は、
   燃料質量割合が「０」の前記筒内ガスの塊、及び燃料質量割合が「１」の燃料蒸気の塊を前記衝突反応の対象となるガス塊として扱うように構成された内燃機関の混合気の不均一度取得装置。
3. 請求項２に記載の内燃機関の混合気の不均一度取得装置において、
   前記不均一度取得手段は、
   前記混合気内における燃料質量割合が異なる２つ以上の前記ガス塊のうちの任意の２つの組み合わせの全てについて前記衝突反応が順に行われるものと仮定して前記不均一度を取得するように構成された内燃機関の混合気の不均一度取得装置。
4. 請求項２又は請求項３に記載の内燃機関の混合気の不均一度取得装置であって、
   前記衝突した２つのガス塊の前記一部分に含まれるそれぞれの燃料の質量の和が、前記２つのガス塊の前記一部分のそれぞれの質量の和に前記他のガス塊の燃料質量割合を乗じた値に等しいという燃料についての質量保存則を利用して、前記他のガス塊の燃料質量割合を特定する燃料質量割合特定手段を更に備えた内燃機関の混合気の不均一度取得装置。
5. 請求項２乃至請求項４の何れか一項に記載の内燃機関の混合気の不均一度取得装置であって、
   前記混合気内における燃料質量割合が異なる２つ以上の前記ガス塊の量をそれぞれ取得するガス塊量取得手段と、
   前記衝突した２つのガス塊の前記それぞれの一部分が混合する体積である混合体積を取得する混合体積取得手段と、
   を更に備え、
   前記ガス塊量取得手段は、
   前記取得された混合体積を用いて前記衝突した２つのガス塊の量における前記一部分に相当するそれぞれの減少分と前記他のガス塊の量における増加分とを算出することで、前記２つのガス塊のそれぞれの量、及び前記他のガス塊の量を更新・取得するように構成されていて、
   前記不均一度取得手段は、
   前記更新・取得されたそれぞれのガス塊の量の分布に基づいて前記不均一度を取得するように構成された内燃機関の混合気の不均一度取得装置。
6. 請求項５に記載の内燃機関の混合気の不均一度取得装置において、
   前記不均一度取得手段は、
   前記不均一度として、前記更新・取得されたそれぞれのガス塊の量の分布に基づいて得られる燃料質量割合についての標準偏差を使用するように構成された内燃機関の混合気の不均一度取得装置。
7. 請求項５又は請求項６に記載の内燃機関の混合気の不均一度取得装置において、
   前記混合体積取得手段は、
   前記衝突した２つのガス塊のそれぞれの量のうち小さい方に基づいて前記混合体積を取得するように構成された内燃機関の混合気の不均一度取得装置。
8. 請求項７に記載の内燃機関の混合気の不均一度取得装置において、
   前記混合体積取得手段は、
   前記衝突した２つのガス塊のうち量の大きい方の体積の前記混合気全体の体積に対する割合に基づいて前記混合体積を取得するように構成された内燃機関の混合気の不均一度取得装置。
9. 請求項５乃至請求項８の何れか一項に記載の内燃機関の混合気の不均一度取得装置であって、
   前記混合気内における燃料質量割合が異なる２つ以上の前記ガス塊の速度をそれぞれ取得するガス塊速度取得手段を更に備え、
   前記混合体積取得手段は、
   前記衝突した２つのガス塊の速度の差に基づいて前記混合体積を取得するように構成された内燃機関の混合気の不均一度取得装置。
10. 請求項９に記載の内燃機関の混合気の不均一度取得装置において、
    前記ガス塊速度取得手段は、
    前記衝突反応前における前記他のガス塊の運動量と前記衝突反応前における前記２つのガス塊の前記一部分のそれぞれの運動量との和が前記衝突反応後における前記他のガス塊の運動量に等しいという運動量保存則を利用して、前記他のガス塊の速度を更新・取得するように構成された内燃機関の混合気の不均一度取得装置。
11. 請求項９又は請求項１０に記載の内燃機関の混合気の不均一度取得装置において、
    前記ガス塊速度取得手段は、
    前記混合気内に新たに取り込んだ前記筒内ガスの速度を「０」として前記筒内ガスの塊の速度を取得するように構成された内燃機関の混合気の不均一度取得装置。
12. 請求項９乃至請求項１１の何れか一項に記載の内燃機関の混合気の不均一度取得装置において、
    前記ガス塊速度取得手段は、
    前記噴射された燃料の速度を利用して前記燃料蒸気の塊の速度を取得するように構成された内燃機関の混合気の不均一度取得装置。
13. 請求項１乃至請求項１２の何れか一項に記載の内燃機関の混合気の不均一度取得装置により取得された前記混合気内における燃料濃度の不均一度を利用して前記混合気の状態を取得する混合気状態取得手段を備えた内燃機関の混合気状態取得装置。
14. 請求項１３に記載の内燃機関の混合気状態取得装置において、
    前記混合気状態取得手段は、
    前記混合気の温度を取得する混合気温度取得手段と、
    前記混合気の温度が所定の着火温度を超えたか否かによって前記混合気内での燃料の反応速度を取得する手法を変更して前記燃料の反応速度を取得する燃料反応速度取得手段と、
    を備えた内燃機関の混合気状態取得装置。
15. 請求項１４に記載の内燃機関の混合気状態取得装置において、
    前記燃料反応速度取得手段は、
    前記混合気の温度が前記所定の着火温度を超えていない場合には低温酸化反応に基づく予混合燃焼を考慮して前記燃料の反応速度を取得するとともに、前記混合気の温度が前記所定の着火温度を超えた場合には高温酸化反応に基づく予混合燃焼と拡散燃焼とを考慮して前記燃料の反応速度を取得するように構成された内燃機関の混合気状態取得装置。
16. 請求項１５に記載の内燃機関の混合気状態取得装置において、
    前記燃料反応速度取得手段は、
    前記混合気の温度が前記所定の着火温度を超えた場合であって燃料の噴射がなお継続する場合、前記混合気を、前記混合気の温度が前記所定の着火温度を超える前に噴射された燃料に基づく着火前噴射部分と前記混合気の温度が前記所定の着火温度を超えた後に噴射された燃料に基づく着火後噴射部分とに区分して個別に扱うとともに、前記着火前噴射部分についての燃料の反応速度を前記高温酸化反応に基づく予混合燃焼と前記拡散燃焼とを考慮して取得し、前記着火後噴射部分についての燃料の反応速度を前記拡散燃焼のみを考慮して取得するように構成された内燃機関の混合気状態取得装置。
    

Images