(57)【特許請求の範囲】
【請求項１】
内燃機関の気筒内の燃焼による熱発生をモデル化するためのWiebe関数の少なくとも一つのパラメータの値を決定する方法であって、
前記内燃機関について、少なくとも燃焼期間が含まれる区間において筒内圧をクランク角毎に計測する筒内圧計測ステップと、
決定すべきパラメータに暫定値を代入したWiebe関数により熱発生率を算出し、該熱発生率に基づいて燃焼期間中のクランク角毎のモデル筒内圧を算出するモデル筒内圧算出ステップと、
前記筒内圧計測ステップで計測された実筒内圧に対する、前記モデル筒内圧の誤差の大きさを表す誤差評価値を算出する誤差評価値算出ステップと、
前記筒内圧計測ステップで計測された実筒内圧に基づいて、実熱発生率を算出するステップと、
前記実熱発生率に基づいて、前記暫定値の初期値を算出する初期暫定値算出ステップと、
を備え、
前記暫定値を変化させて前記モデル筒内圧算出ステップおよび前記誤差評価値算出ステップを繰り返し行い、前記誤差評価値が最小となったときの前記暫定値を前記決定すべきパラメータの確定値とすることを特徴とするWiebe関数パラメータの決定方法。
【請求項２】
前記モデル筒内圧算出ステップでは、熱発生開始点でのモデル筒内圧が前記筒内圧計測ステップで計測された実筒内圧に一致するようにモデル筒内圧を算出することを特徴とする請求項１記載のWiebe関数パラメータの決定方法。
【請求項３】
前記決定すべきパラメータに、形状パラメータｍと、燃焼期間θｐ[degCA]と、筒内に供給された燃料の熱量が筒内での熱発生量に変換される効率ｋと、点火時期とWiebe関数の熱発生開始点とのズレ量を表す熱発生開始点ズレ量θｂ[degCA]とのうちの少なくとも一つを含むことを特徴とする請求項１または２記載のWiebe関数パラメータの決定方法。
【請求項４】
内燃機関の気筒内の燃焼による熱発生をモデル化するためのWiebe関数の少なくとも一つのパラメータの値を決定する装置であって、
前記内燃機関について、少なくとも燃焼期間が含まれる区間において筒内圧をクランク角毎に計測する筒内圧計測手段と、
決定すべきパラメータに暫定値を代入したWiebe関数により熱発生率を算出し、該熱発生率に基づいて燃焼期間中のクランク角毎のモデル筒内圧を算出するモデル筒内圧算出手段と、
前記筒内圧計測手段により計測された実筒内圧に対する、前記モデル筒内圧の誤差の大きさを表す誤差評価値を算出する誤差評価値算出手段と、
前記筒内圧計測手段により計測された実筒内圧に基づいて、実熱発生率を算出する手段と、
前記実熱発生率に基づいて、前記暫定値の初期値を算出する初期暫定値算出手段と、
前記暫定値を変化させて、前記モデル筒内圧算出手段によるモデル筒内圧の算出と、前記誤差評価値算出手段による誤差評価値の算出とを繰り返し行い、誤差評価値が最小となったときの前記暫定値を前記決定すべきパラメータの確定値として取得する確定値取得手段と、
を備えることを特徴とするWiebe関数パラメータの決定装置。
【請求項５】
Wiebe関数を下記（Ｉ）式で表したとき（ただし、下記（Ｉ）式中、Ｑは気筒内の発熱量[J]、Ｑtotalは気筒内に供給された燃料の熱量[J]、θは熱発生開始後の経過クランク角度[degCA]、θｐは燃焼期間[degCA]、ｍは形状パラメータ、ｋは効率、をそれぞれ表す）、前記決定すべきパラメータに前記形状パラメータｍを含み、
前記形状パラメータｍの初期暫定値を求める工程として、
前記筒内圧計測ステップで計測された実筒内圧に基づいて、クランク角毎の実熱発生率を算出するステップと、
前記実熱発生率が最大となったクランク角度での実燃焼割合αを算出するステップと、
前記実燃焼割合αに基づいて、下記（II）式またはこれと等価な式に従い、前記形状パラメータｍの初期暫定値を算出するステップと、
を更に備えることを特徴とする請求項１乃至３の何れか１項記載のWiebe関数パラメータの決定方法。
【数１】
0004577239.tif 000014
【請求項６】
前記決定すべきパラメータに前記燃焼期間θｐを更に含み、
前記燃焼期間θｐの初期暫定値を求める工程として、
前記実熱発生率の積分値を総発熱量ｋＱtotalとして求めるステップと、
前記実熱発生率の最大値をβとしたとき、この実熱発生率最大値βと、前記総発熱量ｋＱtotalと、前記形状パラメータｍの初期暫定値と、所定の定数ａとに基づいて、下記（III）式またはこれと等価な式に従い、前記燃焼期間θｐの初期暫定値を算出するステップと、
を更に備えることを特徴とする請求項５記載のWiebe関数パラメータの決定方法。
【数２】
0004577239.tif 000015
【請求項７】
前記内燃機関の点火時期をSA[degBTDC]とし、この点火時期SAとWiebe関数の熱発生開始点とのズレ量を熱発生開始点ズレ量θｂ[degCA]としたとき、前記決定すべきパラメータに前記熱発生開始点ズレ量θｂを更に含み、
前記熱発生開始点ズレ量θｂの初期暫定値を求める工程として、
前記実熱発生率が最大となったときのクランク角度をθATDC＊[degATDC]としたとき、前記クランク角度θATDC＊と、前記点火時期SAと、前記形状パラメータｍの初期暫定値と、前記燃焼期間θｐの初期暫定値と、所定の定数ａとに基づいて、下記（IV）および（Ｖ）式、またはこれらと等価な式に従い、前記熱発生開始点ズレ量θｂの初期暫定値を算出するステップを更に備えることを特徴とする請求項６記載のWiebe関数パラメータの決定方法。
【数３】
0004577239.tif 000016
【請求項８】
前記決定すべきパラメータに前記効率ｋを更に含み、
前記効率ｋの初期暫定値を求める工程として、
前記実熱発生率の積分値を総発熱量ｋＱtotalとして求めるステップと、
気筒内に供給された燃料の熱量Ｑtotalの値を求めるステップと、
前記総発熱量ｋＱtotalを前記Ｑtotalの値で除することにより、前記効率ｋの初期暫定値を算出するステップと、
を更に備えることを特徴とする請求項５乃至７の何れか１項記載のWiebe関数パラメータの決定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、Wiebe関数パラメータの決定方法および決定装置に係り、特に、内燃機関の気筒内の熱発生をモデル化するためのWiebe関数のパラメータの値を決定するのに好適なWiebe関数パラメータの決定方法および決定装置に関する。
【背景技術】
【０００２】
従来より、内燃機関の気筒内の熱発生をWiebe関数によってモデル化する試みがなされている。Wiebe関数によれば、クランク角度毎の熱発生率や燃焼割合等を算出することができる。そして、クランク角度毎の熱発生率等の値が分かれば、筒内圧や図示トルク等を算出することができる。ゆえに、気筒内の熱発生をWiebe関数によって正確にモデル化することができれば、内燃機関の開発や内燃機関の制御などの分野でWiebe関数を有効に利用することが可能となる。
【０００３】
気筒内の熱発生をWiebe関数によって正確にシミュレートするためには、Wiebe関数に含まれるいくつかのパラメータの値を、対象とする内燃機関の機種に応じて、精度良く決定する必要がある。また、それらのパラメータの値は、空燃比、機関回転数、負荷率、点火時期などの運転条件によっても異なる。
【０００４】
特開２００４−３３２６５９号公報には、Wiebe関数を利用して点火時期をＭＢＴ（Minimum advance for the Best Torque）付近に制御する技術が開示されている。この公報では、同一エンジンのＭＢＴにおいては、燃焼期間の長さにかかわらず熱発生率の形状（波形）が同じになり、かつ熱発生率が最大となるときのクランク角度および燃焼割合も一定になることを仮定している。そして、Wiebe関数パラメータのうち、熱発生率形状に関与する形状パラメータｍを定数とし、それらの値を実験で求めることとしている（上記公報の段落番号００３０参照）。なお、上記公報でのパラメータ「ｎ」が本明細書における「ｍ」に相当する。
【０００５】
【特許文献１】特開２００４−３３２６５９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、上記公報には、形状パラメータｍその他のWiebe関数パラメータを実験データに基づいて決定する具体的な方法については何ら開示されていない。また、本発明者の知見によれば、同一エンジンのＭＢＴにおいても、形状パラメータｍの値は運転条件によって変化し、一定とはならない。また、種々の制御が行われる現代のエンジンにおいては、必ずしも点火時期をＭＢＴに合わせればよいとは限らない。このため、ＭＢＴ以外の点火時期の下でのトルク等を推定する要請も存在し、そのためには点火時期の変化によるWiebe関数パラメータの値の変化傾向を把握する必要がある。
【０００７】
このように、Wiebe関数を内燃機関の制御等に実際に利用するには、様々な運転条件の下でWiebe関数パラメータを精度良く決定し、それらの値が運転条件に応じてどのように変化するかを正確に把握することが必要である。しかしながら、従来、Wiebe関数パラメータを精度良く決定する手法は確立されていないのが実状である。
【０００８】
この発明は、上述のような課題を解決するためになされたもので、内燃機関の気筒内の熱発生をモデル化するためのWiebe関数のパラメータを精度良く決定することのできるWiebe関数パラメータの決定方法および決定装置を提供することを目的とする。
【課題を解決するための手段】
【０００９】
第１の発明は、上記の目的を達成するため、Wiebe関数パラメータの決定方法であって、
内燃機関の気筒内の燃焼による熱発生をモデル化するためのWiebe関数の少なくとも一つのパラメータの値を決定する方法であって、
前記内燃機関について、少なくとも燃焼期間が含まれる区間において筒内圧をクランク角毎に計測する筒内圧計測ステップと、
決定すべきパラメータに暫定値を代入したWiebe関数により熱発生率を算出し、該熱発生率に基づいて燃焼期間中のクランク角毎のモデル筒内圧を算出するモデル筒内圧算出ステップと、
前記筒内圧計測ステップで計測された実筒内圧に対する、前記モデル筒内圧の誤差の大きさを表す誤差評価値を算出する誤差評価値算出ステップと、
を備え、
前記暫定値を変化させて前記モデル筒内圧算出ステップおよび前記誤差評価値算出ステップを繰り返し行い、前記誤差評価値が最小となったときの前記暫定値を前記決定すべきパラメータの確定値とすることを特徴とする。
【００１０】
また、第２の発明は、第１の発明において、
前記モデル筒内圧算出ステップでは、熱発生開始点でのモデル筒内圧が前記筒内圧計測ステップで計測された実筒内圧に一致するようにモデル筒内圧を算出することを特徴とする。
【００１１】
また、第３の発明は、第１または第２の発明において、
前記決定すべきパラメータに、形状パラメータｍと、燃焼期間θｐ[degCA]と、筒内に供給された燃料の熱量が筒内での熱発生量に変換される効率ｋと、点火時期とWiebe関数の熱発生開始点とのズレ量を表す熱発生開始点ズレ量θｂ[degCA]とのうちの少なくとも一つを含むことを特徴とする。
【００１２】
また、第４の発明は、第１乃至第３の発明の何れかにおいて、
前記筒内圧計測ステップで計測された実筒内圧に基づいて、実熱発生率を算出するステップと、
前記実熱発生率に基づいて、前記暫定値の初期値を算出する初期暫定値算出ステップと、
を更に備えることを特徴とする。
【００１３】
また、第５の発明は、Wiebe関数パラメータの決定装置であって、
内燃機関の気筒内の燃焼による熱発生をモデル化するためのWiebe関数の少なくとも一つのパラメータの値を決定する装置であって、
前記内燃機関について、少なくとも燃焼期間が含まれる区間において筒内圧をクランク角毎に計測する筒内圧計測手段と、
決定すべきパラメータに暫定値を代入したWiebe関数により熱発生率を算出し、該熱発生率に基づいて燃焼期間中のクランク角毎のモデル筒内圧を算出するモデル筒内圧算出手段と、
前記筒内圧計測手段により計測された実筒内圧に対する、前記モデル筒内圧の誤差の大きさを表す誤差評価値を算出する誤差評価値算出手段と、
前記暫定値を変化させて、前記モデル筒内圧算出手段によるモデル筒内圧の算出と、前記誤差評価値算出手段による誤差評価値の算出とを繰り返し行い、誤差評価値が最小となったときの前記暫定値を前記決定すべきパラメータの確定値として取得する確定値取得手段と、
を備えることを特徴とする。
【００１４】
また、第６の発明は、第１乃至第４の発明の何れかにおいて、
Wiebe関数を下記（Ｉ）式で表したとき（ただし、下記（Ｉ）式中、Ｑは気筒内の発熱量[J]、Ｑtotalは気筒内に供給された燃料の熱量[J]、θは熱発生開始後の経過クランク角度[degCA]、θｐは燃焼期間[degCA]、ｍは形状パラメータ、ｋは効率、をそれぞれ表す）、前記決定すべきパラメータに前記形状パラメータｍを含み、
前記形状パラメータｍの初期暫定値を求める工程として、
前記筒内圧計測ステップで計測された実筒内圧に基づいて、クランク角毎の実熱発生率を算出するステップと、
前記実熱発生率が最大となったクランク角度での実燃焼割合αを算出するステップと、
前記実燃焼割合αに基づいて、下記（II）式またはこれと等価な式に従い、前記形状パラメータｍの初期暫定値を算出するステップと、
を更に備えることを特徴とする。
【数１】
0004577239.tif 000002
【００１５】
また、第７の発明は、第６の発明において、
前記決定すべきパラメータに前記燃焼期間θｐを更に含み、
前記燃焼期間θｐの初期暫定値を求める工程として、
前記実熱発生率の積分値を総発熱量ｋＱtotalとして求めるステップと、
前記実熱発生率の最大値をβとしたとき、この実熱発生率最大値βと、前記総発熱量ｋＱtotalと、前記形状パラメータｍの初期暫定値と、所定の定数ａとに基づいて、下記（III）式またはこれと等価な式に従い、前記燃焼期間θｐの初期暫定値を算出するステップと、
を更に備えることを特徴とする。
【数２】
0004577239.tif 000003
【００１６】
また、第８の発明は、第７の発明において、
前記内燃機関の点火時期をSA[degBTDC]とし、この点火時期SAとWiebe関数の熱発生開始点とのズレ量を熱発生開始点ズレ量θｂ[degCA]としたとき、前記決定すべきパラメータに前記熱発生開始点ズレ量θｂを更に含み、
前記熱発生開始点ズレ量θｂの初期暫定値を求める工程として、
前記実熱発生率が最大となったときのクランク角度をθATDC＊[degATDC]としたとき、前記クランク角度θATDC＊と、前記点火時期SAと、前記形状パラメータｍの初期暫定値と、前記燃焼期間θｐの初期暫定値と、所定の定数ａとに基づいて、下記（IV）および（Ｖ）式、またはこれらと等価な式に従い、前記熱発生開始点ズレ量θｂの初期暫定値を算出するステップを更に備えることを特徴とする。
【数３】
0004577239.tif 000004
【００１７】
また、第９の発明は、第６乃至第８の発明の何れかにおいて、
前記決定すべきパラメータに前記効率ｋを更に含み、
前記効率ｋの初期暫定値を求める工程として、
前記実熱発生率の積分値を総発熱量ｋＱtotalとして求めるステップと、
気筒内に供給された燃料の熱量Ｑtotalの値を求めるステップと、
前記総発熱量ｋＱtotalを前記Ｑtotalの値で除することにより、前記効率ｋの初期暫定値を算出するステップと、
を更に備えることを特徴とする。
【発明の効果】
【００１８】
第１の発明によれば、実筒内圧と、Wiebe関数に基づいて算出されるモデル筒内圧との誤差がなるべく小さくなるように、Wiebe関数パラメータの値を決定することができる。このため、第１の発明によれば、内燃機関の筒内圧やトルクをWiebe関数モデルによって予測する上で、内燃機関の機種に対するWiebe関数パラメータの依存性や、運転条件とWiebe関数パラメータとの関係を精度良く把握することができる。その結果、Wiebe関数モデルを用いて、内燃機関の筒内圧やトルクを正確に予測することが可能となる。
【００１９】
第２の発明によれば、Wiebe関数パラメータの決定過程において、熱発生開始点でのモデル筒内圧が実筒内圧に一致するようにモデル筒内圧を算出することができる。これにより、Wiebe関数パラメータの値をより高い精度で決定することができる。
【００２０】
第３の発明によれば、Wiebe関数パラメータとして、形状パラメータｍと、燃焼期間θｐ[degCA]と、筒内に供給された燃料の熱量が筒内での熱発生量に変換される効率ｋと、点火時期とWiebe関数の熱発生開始点とのズレ量を表す熱発生開始点ズレ量θｂ[degCA]とのうちの少なくとも一つを精度良く決定することができる。
【００２１】
第４の発明によれば、Wiebe関数パラメータの決定過程において用いるWiebe関数パラメータの暫定値の初期値を、実熱発生率に基づいて精度良く算出することができる。このため、誤差評価値を最小とするようなWiebe関数パラメータの確定値（適合値）を容易に探し出すことができ、計算負荷を軽減することができる。
【００２２】
第５の発明によれば、実筒内圧と、Wiebe関数に基づいて算出されるモデル筒内圧との誤差がなるべく小さくなるように、Wiebe関数パラメータの値を決定することができる。このため、第５の発明によれば、内燃機関の筒内圧やトルクをWiebe関数モデルによって予測する上で、内燃機関の機種に対するWiebe関数パラメータの依存性や、運転条件とWiebe関数パラメータとの関係を精度良く把握することができる。その結果、Wiebe関数モデルを用いて、内燃機関の筒内圧やトルクを正確に予測することが可能となる。
【００２３】
第６の発明によれば、Wiebe関数パラメータの一つである形状パラメータｍの初期暫定値を精度良く算出することができる。このため、形状パラメータｍの最終的な確定値を探索する上で、精度の良い初期暫定値を出発点とすることができる。よって、最終的な確定値を容易に見つけ出すことができ、計算負荷を軽減することができる。
【００２４】
第７の発明によれば、Wiebe関数パラメータの一つである燃焼期間θｐの初期暫定値を精度良く算出することができる。このため、燃焼期間θｐの最終的な確定値を探索する上で、精度の良い初期暫定値を出発点とすることができる。よって、最終的な確定値を容易に見つけ出すことができ、計算負荷を軽減することができる。
【００２５】
第８の発明によれば、Wiebe関数パラメータの一つである熱発生開始点ズレ量θｂの初期暫定値を精度良く算出することができる。このため、熱発生開始点ズレ量θｂの最終的な確定値を探索する上で、精度の良い初期暫定値を出発点とすることができる。よって、最終的な確定値を容易に見つけ出すことができ、計算負荷を軽減することができる。
【００２６】
第９の発明によれば、Wiebe関数パラメータの一つである効率ｋの初期暫定値を精度良く算出することができる。このため、効率ｋの最終的な確定値を探索する上で、精度の良い初期暫定値を出発点とすることができる。よって、最終的な確定値を容易に見つけ出すことができ、計算負荷を軽減することができる。
【発明を実施するための最良の形態】
【００２７】
実施の形態１．
［システム構成の説明］
図１は、本発明の実施の形態１で用いられるシステム構成を説明するための図である。図１に示すように、本実施形態のシステムは、火花点火式の内燃機関１０を備えている。内燃機関１０の各気筒には、吸気弁１４、排気弁１６および点火プラグ３０が設けられている。
【００２８】
また、内燃機関１０には、クランク軸の回転位置（クランク角度）を検出するクランク角センサ１２が設けられている。クランク角センサ１２は、クランク軸が所定回転角だけ回転する毎に、Hi出力とLo出力を反転させるセンサである。クランク角センサ１２の出力によれば、機関回転数NEを検知することもできる。
【００２９】
また、内燃機関１０には、筒内圧センサ１８が組み込まれている。筒内圧センサ１８によれば、気筒内（燃焼室内）に生ずる圧力を検出することができる。
【００３０】
内燃機関１０の吸気通路１９の途中には、サージタンク２０が設けられている。サージタンク２０には、その内部の圧力、すなわち吸気管圧力を検出する吸気圧センサ２１が設置されている。また、吸気通路１９には、その内部を流通する吸入空気量GAを検出するエアフローメータ２２が配置されている。吸気圧センサ２１の出力によれば、内燃機関１０の負荷率KL［％］などを取得することができる。なお、負荷率KLは、吸入空気量GAと機関回転数NEとに基づいて算出してもよい。
【００３１】
エアフローメータ２２の下流には、スロットル弁２４が配置されている。スロットル弁２４の近傍には、スロットル開度TAを検出するスロットル開度センサ２６が組み付けられている。
【００３２】
内燃機関１０の吸気ポートには、ガソリン等の燃料を噴射するためのインジェクタ２８が配設されている。なお、本発明における内燃機関は、このようなポート噴射式の内燃機関に限定されるものではなく、筒内直接噴射式の内燃機関や、ポート噴射と筒内噴射とを併用する内燃機関であってもよい。
【００３３】
内燃機関１０の排気通路３２には、排気ガスの空燃比を検出する空燃比センサ３３が設置されている。また、排気通路３２には、排気ガスを浄化するための触媒３４が組み込まれている。
【００３４】
本実施形態のシステムは、制御装置として、ECU(Electronic Control Unit)５０を備えている。ECU５０には、上述した各種のセンサからセンサ信号が供給されている。ECU５０は、それらのセンサ信号に基づいて、スロットル弁２４や、インジェクタ２８、点火プラグ３０などの各種アクチュエータを制御することができる。更に、ECU５０は、クランク角センサ１２、筒内圧センサ１８等のセンサ信号に基づいて、後述するWiebe関数パラメータの決定処理を行うことができる。
【００３５】
（Wiebe関数）
本実施形態のWiebe関数は、下記（１）式あるいは下記（２）式で表される。
【数４】
0004577239.tif 000005
【００３６】
上記（１）および（２）式中の各記号の意味は、次の通りである。
Ｑ：内燃機関１０の気筒内の発熱量[J]
Ｑtotal：気筒内に供給された燃料の熱量[J]
θ：熱発生開始後の経過クランク角度[degCA]
θｐ：燃焼期間[degCA]
ｍ：形状パラメータ
ｋ：効率
【００３７】
上記（１）式は、熱発生率dＱ/dθを算出するためのWiebe関数である。上記（２）式の左辺のＱ/ｋＱtotalは、後述するように、燃焼割合を表す。つまり、上記（２）式は、燃焼割合Ｑ/ｋＱtotalを算出するためのWiebe関数である。
【００３８】
図２（Ａ）は、横軸にクランク角、縦軸に熱発生率をとったグラフであり、図２（Ｂ）は、横軸にクランク角、縦軸に筒内圧をとったグラフである。図２中、細線は実機での計測値であり、太線はWiebe関数に基づいて算出される値である。なお、図２（Ａ）と図２（Ｂ）とは別々のグラフであり、クランク角を表す横軸の位置やスケールは一致していない。
【００３９】
図１に示すシステムでは、筒内圧センサ１８により、図２（Ｂ）中の細線で示すような実機の筒内圧をクランク角毎に計測することができる。以下の説明では、この計測値を「実筒内圧」と称する。後述するように、実筒内圧からは、エネルギー保存則の式を用いて、クランク角毎の熱発生率を算出することができる。これが図２（Ａ）中の細線で示す実機の熱発生率である。以下の説明では、実筒内圧から算出された熱発生率を「実熱発生率」と称する。
【００４０】
一方、上記（１）式のWiebe関数からは、図２（Ａ）中の太線で示すようなクランク角毎の熱発生率を算出することができる。そして、Wiebe関数により算出された熱発生率（以下、「モデル熱発生率」と称する）からは、エネルギー保存則の式を用いて、クランク角毎の筒内圧を算出することができる。これが図２（Ｂ）中の太線で示す筒内圧である。以下の説明では、このようにしてWiebe関数に基づいて算出される筒内圧を「モデル筒内圧」と称する。
【００４１】
以下、Wiebe関数の各パラメータについて説明する。
（燃焼期間θｐ）
図２（Ａ）に示すように、燃焼期間θｐは、燃焼による熱の発生が継続する期間、つまり燃焼の開始から終了までの期間をクランク角度で表したものとしての物理的意味を有している。上記（１）式および（２）式のWiebe関数において、θ＝０の点は後述するように熱発生開始点（燃焼開始点）を意味しており、これらのWiebe関数は０≦θ≦θｐの範囲で定義される。
【００４２】
（形状パラメータｍ）
形状パラメータｍは、図２（Ａ）のようなグラフにおけるモデル熱発生率dＱ/dθの形状に大きく関与するパラメータであり、特に、熱発生率が最大値をとるクランク角度に大きく関与するパラメータである。
【００４３】
（熱発生開始点ズレ量θｂ）
図２（Ａ）に示すように、上記（１）式で表されるWiebe関数においては、θ＝０のときにモデル熱発生率dＱ/dθ＝０であり、θが０を超えた時点からdＱ/dθ＞０となり、熱発生が開始することとなる。つまり、本実施形態のWiebe関数では、θ＝０が熱発生開始点（燃焼開始点）であり、θは熱発生開始後（燃焼開始後）の経過クランク角度を表すことになる。
【００４４】
本発明者の知見によれば、Wiebe関数の熱発生開始点を実際の点火時期SAに一致させることとすると、Wiebe関数によるモデル熱発生率の波形を実熱発生率の波形に精度良く適合させることができない。そこで、本実施形態では、Wiebe関数の熱発生開始点（θ＝０）を実際の点火時期SAと異なる点とすることとした。そして、その両者のズレ量を熱発生開始点ズレ量θｂとして表し、これをWiebe関数のパラメータとして導入することとした。すなわち、熱発生開始点ズレ量θｂ[degCA]は、Wiebe関数の熱発生開始点と、点火時期SAとのズレ量をクランク角度で表すパラメータである（図２（Ａ）参照）。以下、本実施形態では、上記（１）式等のWiebe関数を用いる場合、点火時期SAから熱発生開始点ズレ量θｂ分だけずれた点を基準点（θ＝０）として上記（１）式等を適用するものとする。
【００４５】
（効率ｋ）
上記（１）および（２）式中、気筒内に供給された燃料の熱量Ｑtotalは、筒内に供給された燃料量に、その燃料の低位発熱量を乗ずることにより算出することができる。なお、低位発熱量は、真発熱量とも呼ばれる物性値である。低位発熱量とは、単位量の燃料が完全燃焼したときに発生する熱量から、燃料中に含まれる水分および燃焼によって生じる水分を蒸発させるのに必要な熱量（潜熱）を差し引いた残りの熱量を意味する。図１に示すシステムにおいては、Ｑtotalの値は、燃料の低位発熱量が既知であるものとして、インジェクタ２８からの燃料噴射量に基づいて算出することができる。あるいは、空燃比A/Fおよび筒内空気量（負荷率KL）からＱtotalを算出することもできる。
【００４６】
内燃機関における燃焼では、冷却損失や燃料の燃え残りなどに起因して、何らかの熱損失を伴うのが普通である。このため、気筒内に供給された燃料の熱量Ｑtotalが全部そのまま気筒内の発熱量Ｑに変換されることは実際上はあり得ない。本実施形態では、このことをWiebe関数に反映させるパラメータとして、効率ｋを導入することとした。すなわち、効率ｋは、気筒内に供給された燃料の熱量Ｑtotalが発熱量Ｑに変換される効率としての物理的意味を有しており、０＜ｋ＜１なる範囲の数である。
【００４７】
効率ｋを用いると、気筒内での総発熱量は、Ｑtotalにｋを乗じたｋＱtotalとして表すことができる。総発熱量ｋＱtotalは、図２（Ａ）において、実熱発生率を示す細線と、dＱ/dθ＝０の直線とがグラフ上で囲む面積に相当する。つまり、総発熱量ｋＱtotalの値は、実熱発生率を積分することによって算出することが可能である。そのようにして算出された総発熱量ｋＱtotalを、上述の方法で算出したＱtotalで除することにより、効率ｋの値を求めることができる。
【００４８】
上記のように総発熱量がｋＱtotalとして表されるので、あるクランク角度θにおける燃焼割合は、そのクランク角度θまでの発熱量Ｑと、総発熱量ｋＱtotalとの比Ｑ/ｋＱtotalとして表すことができる。この燃焼割合Ｑ/ｋＱtotalは、Wiebe関数によれば、上記（２）式で表される。上記（２）式をθで微分すると、上記（１）式が得られる。つまり、上記（２）式は、上記（１）式と等価である。
【００４９】
（パラメータａ）
燃焼割合の本来の定義からすると、燃焼割合の値は、燃焼終了時に１００％、すなわち１となるべきである。これに対し、Wiebe関数では、その性質上、燃焼終了時に燃焼割合の値が１とならない。すなわち、Wiebe関数では、燃焼終了時の燃焼割合は、上記（２）式にθ＝θｐを代入することにより、下記式で表される。
Ｑ/ｋＱtotal（燃焼終了時）＝１−exp（−ａ）
【００５０】
上記式のように、Wiebe関数上での燃焼終了時の燃焼割合Ｑ/ｋＱtotalの値は、パラメータａのみで表される。そこで、上記式の値｛１−exp（−ａ）｝を１に近い所定の値と等しいと置けば、ａの値を定めることができる。例えば、｛１−exp（−ａ）｝＝０．９９９と置けば、ａ≒６．９が得られる。そこで、以下の説明では、パラメータａを定数６．９として扱うものとする。ただし、本発明は、パラメータａを定数として扱う場合に限定されるものではなく、ａの値を他のパラメータとともに後述する手法で決定するようにしてもよい。
【００５１】
［実施の形態１の特徴］
本実施形態では、パラメータａを定数として扱うので、決定すべきWiebe関数パラメータは、形状パラメータｍ、効率ｋ、熱発生開始点ズレ量θｂ、燃焼期間θｐの４つとなる。以下、本実施形態のWiebe関数パラメータ決定方法の作用・効果を分かり易くするため、まず、比較例の方法について説明する。
【００５２】
（比較例のWiebe関数パラメータ決定方法）
比較例としては、図２（Ａ）に示すような実熱発生率とモデル熱発生率との誤差の大きさを表す誤差評価値が最小となるように、Wiebe関数パラメータｍ、ｋ、θｂ、θｐの各々の値を最適化する方法が考えられる。例えば最小二乗法を用いる場合には、クランク角毎の実熱発生率とモデル熱発生率との誤差（偏差）の二乗和が誤差評価値とされる。そして、この誤差評価値を、ｍ、ｋ、θｂ、θｐを変数とする多変数関数と考えて、ｍ、ｋ、θｂ、θｐの各々の値を少しずつ変化させながら誤差評価値を繰り返し計算し、誤差評価値が最小となるような一組のWiebe関数パラメータｍ、ｋ、θｂ、θｐの値を探し出し、それらの値を最終的な値として決定する。
【００５３】
上記比較例の方法によれば、実熱発生率とモデル熱発生率との誤差を最小とすることができる。しかしながら、実熱発生率の形状（波形）とモデル熱発生率の形状（波形）とは、完全に一致するものではない。このため、上記比較例の方法で決定したWiebe関数パラメータｍ、ｋ、θｂ、θｐの値を用いてモデル筒内圧を計算した場合、図２（Ｂ）に示すように、実筒内圧との誤差が十分に小さくならない事態が生じ易い。
【００５４】
ところで、内燃機関１０の開発段階や内燃機関１０の制御装置において、Wiebe関数を利用する上では、内燃機関１０の筒内圧やトルクを精度良く予測することが求められる場合が多い。そして、内燃機関１０の筒内圧が正確に予測できれば、内燃機関１０のトルクを正確に予測することが可能である。このため、Wiebe関数の利用実態を考えた場合、実熱発生率とモデル熱発生率との誤差を最小とするよりも、むしろ、実筒内圧とモデル筒内圧との誤差を最小とするようにWiebe関数パラメータを決定した方が合理的である。
【００５５】
（本実施形態のWiebe関数パラメータ決定方法）
そこで、本実施形態では、実筒内圧とモデル筒内圧との誤差の大きさを表す誤差評価値が最小となるように、Wiebe関数パラメータｍ、ｋ、θｂ、θｐの各々の値を最適化することとした。この誤差評価値は、例えば、クランク角毎の実筒内圧とモデル筒内圧との誤差（偏差）の二乗和とされる。そして、その誤差評価値を、ｍ、ｋ、θｂ、θｐを変数とする多変数関数と考えて、ｍ、ｋ、θｂ、θｐの各々の値を少しずつ変化させながら誤差評価値を繰り返し計算し、誤差評価値が最小となるような一組のWiebe関数パラメータｍ、ｋ、θｂ、θｐの値を探し出すこととした。
【００５６】
実筒内圧とモデル筒内圧との間の誤差評価値を算出するためには、まず、Wiebe関数パラメータｍ、ｋ、θｂ、θｐに暫定的な値（暫定値）を代入した上記（１）式のWiebe関数を用いて、モデル筒内圧を算出する必要がある。そして、その誤差評価値を最小とするようなWiebe関数パラメータｍ、ｋ、θｂ、θｐの適合値を探し出すには、そのような適合値が見つかるまで、ｍ、ｋ、θｂ、θｐの暫定値を徐々に変化させながら、モデル筒内圧の算出処理と、誤差評価値の算出処理とを繰り返す必要がある。
【００５７】
このような処理を行う場合、最終的な適合値を探し出すまでに要する計算負荷を軽減する観点からは、Wiebe関数パラメータｍ、ｋ、θｂ、θｐの暫定値の初期値（以下、「初期暫定値」と称する）が、最終的な適合値になるべく近い、精度の良い値であることが重要である。そこで、本実施形態では、実熱発生率に基づいて、Wiebe関数パラメータｍ、ｋ、θｂ、θｐの初期暫定値を精度良く算出することとした。以下、本実施形態におけるWiebe関数パラメータｍ、ｋ、θｂ、θｐの初期暫定値の算出方法について説明する。
【００５８】
（形状パラメータｍの初期暫定値の算出方法）
形状パラメータｍの初期暫定値は、熱発生率が最大となったときの燃焼割合に着目することにより、算出することができる。図３は、Wiebe関数により算出されるモデル熱発生率のグラフである。モデル熱発生率が最大となったときの燃焼割合は、図３のモデル熱発生率のグラフが囲む面積のうち、モデル熱発生率が最大となったクランク角度までの面積を、全部の面積（図３中のハッチング部分）で除した値として算出される。
【００５９】
この値を算出するには、まず、上記（１）式をθで微分し、d２Ｑ/dθ２＝０とおくことにより、モデル熱発生率dＱ/dθが最大となるクランク角度θ＊を求める。下記の式展開を行うことにより、クランク角度θ＊を表す式が得られる。
【数５】
0004577239.tif 000006
【００６０】
上記（３）式で表されるクランク角度θ＊を上記（２）式に代入することにより、モデル熱発生率dＱ/dθが最大となるクランク角度θ＊での燃焼割合は、下記（４）式で表すことができる。
【数６】
0004577239.tif 000007
【００６１】
上記（４）式に示されるように、モデル熱発生率が最大となるクランク角度での燃焼割合は、形状パラメータｍのみの関数で表される。そこで、これに対応する実測値、つまり、実熱発生率が最大となるクランク角度での燃焼割合を求め、その値が上記（４）式の右辺と等しいとおけば、ｍを算出することができる。
【００６２】
ここで、実熱発生率が最大となるクランク角度での燃焼割合（以下、「実燃焼割合」と称する）をαとする。本実施形態では、この実燃焼割合αを上記（４）式の右辺と等しいとおくことにより、形状パラメータｍの初期暫定値を下記（５）式で算出することができる。
【数７】
0004577239.tif 000008
【００６３】
（燃焼期間θｐの初期暫定値の算出方法）
上記（３）式で表されるクランク角度θ＊を上記（１）式に代入することにより、モデル熱発生率の最大値（dＱ/dθ）maxは、下記（６）式で表すことができる。
【数８】
0004577239.tif 000009
【００６４】
よって、実熱発生率の最大値をβとすると、この最大値βが上記（６）式の下段と等しいとおくことにより、下記式が得られる。
【数９】
0004577239.tif 000010
【００６５】
上記（７）式の右辺において、ａ/θｐm+1以外の部分は、総発熱量ｋＱtotalおよび形状パラメータｍのみで表されている。総発熱量ｋＱtotalの値は、前述したように、実熱発生率を積分することによって算出することが可能である。そして、ｍの初期暫定値は、上記（５）式により算出することができる。よって、それらの値と、実熱発生率最大値βとを上記（７）式に代入することにより、ａ/θｐm+1の値を算出することができる。また、本実施形態では、前述したように、ａ＝６．９とすることにしている。よって、ａ/θｐm+1の値が算出できれば、燃焼期間θｐの値を算出することができる。本実施形態では、以上のようにして算出される値を燃焼期間θｐの初期暫定値にすることとした。
【００６６】
（熱発生開始点ズレ量θｂの初期暫定値の算出方法）
本実施形態では、モデル熱発生率が最大となる位置と、実熱発生率が最大となる位置とを一致させることにより、熱発生開始点ズレ量θｂの初期暫定値を算出する。すなわち、実熱発生率が最大となったときの圧縮上死点基準のクランク角度をθATDC＊[degATDC]としたとき、モデル熱発生率が最大となるクランク角度θ＊の位置がθATDC＊に一致するようにする。
【００６７】
図３を参照して分かる通り、点火時期をSA[degBTDC]とし、θ＊とθATDC＊との位置が一致するものとすると、SA、θATDC＊、θｂ、およびθ＊の間に、下記（９）式が成り立つ。また、上記（３）式は、下記（８）式のように変形することができる。
【数１０】
0004577239.tif 000011
【００６８】
上記（８）式によれば、先に算出されたｍの値とａ/θｐm+1の値とに基づいてクランク角度θ＊を算出することができる。そして、そのθ＊と、点火時期SAと、実熱発生率データから得られるθATDC＊とを上記（９）式に代入することにより、熱発生開始点ズレ量θｂの初期暫定値を算出することができる。
【００６９】
以上の方法により、Wiebe関数パラメータｍ、θｐ、θｂの各々の初期暫定値を高い精度で算出することができる。なお、効率ｋの初期暫定値は、効率ｋの説明をした箇所で既述した方法によって算出するものとする。ところで、上記（１）式のWiebe関数は、下記式のように変形することができる。
【数１１】
0004577239.tif 000012
【００７０】
上記（１０）式によれば、ｍ、ｋ、ａ/θｐm+1、およびθｂの値が決定されれば、Wiebe関数を具体的に記述することができることが分かる。本実施形態では、前述したように、ａを定数として扱い、ｍ、ｋ、θｐ、θｂの４つをWiebe関数パラメータとしているが、Wiebe関数パラメータの組み合わせはこれに限られるものではない。例えば、Wiebe関数を上記（１０）式の形で使用することにより、ａおよびθｐの個々の値を決定することなく、両者をａ/θｐm+1と一まとめにした上で、ｍ、ｋ、ａ/θｐm+1、θｂの４つの値をWiebe関数パラメータとして扱うようにしてもよい。
【００７１】
［実施の形態１における具体的処理］
図４は、上述した原理に基づいてWiebe関数パラメータの初期暫定値を算出するために、本実施形態においてECU５０が実行するルーチンのフローチャートである。図５は、実熱発生率と実筒内圧との関係を示す図である。
【００７２】
図４に示すルーチンによれば、まず、内燃機関１０の運転中に、クランク角センサ１２および筒内圧センサ１８の出力に基づいて、図５に示すようなクランク角度毎の実筒内圧Ｐdata[Pa]が計測される（ステップ１００）。このステップ１００において、実筒内圧Ｐdataを計測するクランク角度の間隔は、特に限定されず、例えば１degCA間隔とすることができる。また、実筒内圧Ｐdataの計測は、全クランク角度範囲に渡って行わなくてもよく、少なくとも実際の燃焼期間が含まれる範囲で行えばよい。
【００７３】
次に、上記ステップ１００で実測されたクランク角度毎の実筒内圧Ｐdataに基づいて、次のような方法により、クランク角度毎の実熱発生率が算出される（ステップ１０２）。ここで、熱力学上のエネルギー保存則は、熱発生率dＱ/dθと、筒内圧Ｐと、筒内容積Ｖ[m3]と、クランク角度θとを用いると、次式で表される。なお、ここでのクランク角度θは、上記（１）式等のWiebe関数におけるθとは異なり、上死点基準の絶対的なクランク角度を表すものとする。
【数１２】
0004577239.tif 000013
【００７４】
上記（１１）式中、κは比熱比であり、燃焼ガスの組成などに基づいて定まる既知の値である。また、筒内容積Ｖおよびその変化率dＶ/dθは、クランク角度θに応じて幾何学的に決まる値である。よって、これらの既知の値と、上記ステップ１００で実測されたクランク角度毎の実筒内圧Ｐdataの値とを上記（１１）式に代入することにより、図５に示すようなクランク角度毎の実熱発生率を算出することができる。
【００７５】
次いで、上記ステップ１０２で算出された実熱発生率に基づいて、総発熱量ｋＱtotalが算出される（ステップ１０４）。具体的には、上記ステップ１０２で算出されたクランク角度毎の実熱発生率を積分（積算）することにより、総発熱量ｋＱtotalが算出される。
【００７６】
次いで、効率ｋの初期暫定値が算出される（ステップ１０６）。具体的には、まず、インジェクタ２８からの燃料噴射量等に基づいて、気筒内に供給された燃料の熱量Ｑtotalが算出される。そして、上記ステップ１０４で算出された総発熱量ｋＱtotalをこのＱtotalで除することにより、効率ｋの初期暫定値が算出される。
【００７７】
次いで、上記ステップ１０２で算出されたクランク角度毎の実熱発生率のうちの最大値およびそのときのクランク角度を見つけ出すことにより、実熱発生率最大値βと、そのときのクランク角度θATDC＊とが取得される（ステップ１０８）。
【００７８】
次いで、実熱発生率が最大となったクランク角度θATDC＊での実燃焼割合αが算出される（ステップ１１０）。具体的には、上記ステップ１０２で算出されたクランク角度毎の実熱発生率をクランク角度θATDC＊まで積分（積算）した値を、上記ステップ１０４で算出された総発熱量ｋＱtotalで除することにより、実燃焼割合αが算出される。
【００７９】
次いで、形状パラメータｍの初期暫定値が算出される（ステップ１１２）。具体的には、上記ステップ１１０で算出された実燃焼割合αを上記（５）式に代入することにより、ｍの初期暫定値が算出される。
【００８０】
次いで、燃焼期間θｐの初期暫定値が算出される（ステップ１１４）。具体的には、まず、上記ステップ１０４で算出された総発熱量ｋＱtotalと、上記ステップ１０８で取得された実熱発生率最大値βと、上記ステップ１１２で算出された形状パラメータｍの初期暫定値とを上記（７）式に代入することにより、ａ/θｐm+1が算出される。その後、この算出されたａ/θｐm+1の値と、ａ＝６．９と、形状パラメータｍの初期暫定値とを用いることにより、燃焼期間θｐの初期暫定値が算出される。
【００８１】
最後に、熱発生開始点ズレ量θｂの初期暫定値が算出される（ステップ１１６）。具体的には、ECU５０が設定している点火時期SAと、上記ステップ１０８で取得された実熱発生率が最大となるクランク角度θATDC＊と、上記ステップ１１２で算出された形状パラメータｍの初期暫定値と、上記ステップ１１４で算出されたａ/θｐm+1とを上記（８）および（９）式に代入することにより、熱発生開始点ズレ量θｂの初期暫定値が算出される。
【００８２】
以上説明した図４のルーチンにより、Wiebe関数パラメータであるｍ、ｋ、θｐ、θｂの初期暫定値が求められたら、次に、図６に示すルーチンの処理により、ｍ、ｋ、θｐ、θｂの最終的な適合値（確定値）を求める処理が実行される。図６は、Wiebe関数パラメータの確定値を求めるために、本実施形態においてECU５０が実行するルーチンのフローチャートである。
【００８３】
図６に示すルーチンによれば、まず、上記（１）式のWiebe関数により、燃焼期間θｐ中のクランク角毎のモデル熱発生率dＱ/dθが算出される（ステップ１２０）。この計算に際して、Wiebe関数パラメータのｍ、ｋ、θｐ、θｂには、それぞれ、暫定値が代入される。以下では、一組のWiebe関数パラメータの暫定値［ｍ ｋ θｐ θｂ］を記号ｘで表すことにする。つまり、ｘ＝［ｍ ｋ θｐ θｂ］である。ステップ１２０の処理が最初に実行される場合には、この暫定値ｘとして、図４に示すルーチンの処理によって算出された初期暫定値が用いられる。
【００８４】
次いで、上記ステップ１２０で算出されたクランク角毎のモデル熱発生率dＱ/dθと、上記（１１）式とに基づいて、図５中の破線で示すようなクランク角毎のモデル筒内圧Ｐmodelが算出される（ステップ１２２）。具体的には、上記ステップ１２０で算出されたクランク角毎のモデル熱発生率dＱ/dθの値と、上記（１１）式の右辺とを等しいと置くことにより、筒内圧Ｐの微分方程式が生成され、この微分方程式をECU５０が数値計算によって近似的に解いた値がモデル筒内圧Ｐmodelとされる。この計算処理において、κ、Ｖ、dＶ/dθは前述したように既知の値であるので、ECU５０に予め記憶されている値を用いればよい。
【００８５】
また、上記ステップ１２２の処理において、筒内圧Ｐの微分方程式を解いてモデル筒内圧Ｐmodelを算出する際には、初期値（初期条件）を定める必要がある。すなわち、Wiebe関数の開始点である熱発生開始点（燃焼開始点）におけるモデル筒内圧Ｐmodelの値を定める必要がある。本実施形態では、図５に示すように、このモデル筒内圧Ｐmodelの初期値としては、図４のルーチンのステップ１００で計測された実筒内圧Ｐdataの、対応するクランク角度での値を読み込んで、その値を用いることとした。なお、熱発生開始点は、熱発生開始点ズレ量θｂによって変化するので、モデル筒内圧Ｐmodelの初期値として読み込むべき値は、熱発生開始点ズレ量θｂの暫定値に応じて変化する。このような処理によれば、熱発生開始点において実筒内圧Ｐdataとモデル筒内圧Ｐmodelとを常に一致させることができるので、モデル筒内圧Ｐmodelを実筒内圧Ｐdataに適合させる上で合理的である。
【００８６】
上記ステップ１２２の処理によってモデル筒内圧Ｐmodelが算出されたら、次に、図４のルーチンのステップ１００で計測された実筒内圧Ｐdataと、モデル筒内圧Ｐmodelとの誤差の大きさを表す指標となる誤差評価値Error（ｘ）が算出される。この誤差評価値Error（ｘ）は、本実施形態では、燃焼期間θｐ中のクランク角度毎の実筒内圧Ｐdataとモデル筒内圧Ｐmodelとの誤差（偏差）の二乗和として算出されるものとするが、他の演算手法により算出するようにしてもよい。
【００８７】
誤差評価値Error（ｘ）は、上記ステップ１２０でWiebe関数に代入された一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］の関数として考えることができる。このため、多変数関数の極小値を求める場合と同様の考え方により、誤差評価値Error（ｘ）が最小となるような一組の暫定値ｘを探し出すことが可能である。つまり、一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］の各々について、その値を初期暫定値の前後で少しずつ変化させながら上記ステップ１２０〜１２４の処理を繰り返し行って誤差評価値Error（ｘ）を算出し、誤差評価値Error（ｘ）が極小となるような暫定値を探し出すことにより、誤差評価値Error（ｘ）を最小とするような一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］を探し出すことが可能である。
【００８８】
ステップ１２６では、上記のような考え方に基づき、誤差評価値Error（ｘ）が最小となるような一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］が発見されたか否かが判別される。そして、そのような一組の暫定値ｘが発見されていないと判別された場合には、一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］を更に変化させ（ステップ１２８）、その変化させたｘ＝［ｍ ｋ θｐ θｂ］を用いて、上記ステップ１２０〜１２４の処理が再度実行される。
【００８９】
このようにして上記ステップ１２８を通るループの処理が繰り返し行われることにより、誤差評価値Error（ｘ）を最小とするような一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］を発見することができる。そのような暫定値ｘ＝［ｍ ｋ θｐ θｂ］が発見されて、上記ステップ１２６の判断が肯定されると、次に、その発見された一組の暫定値ｘ＝［ｍ ｋ θｐ θｂ］、つまり誤差評価値Error（ｘ）を最小とするような一組の暫定値ｘ＝［ｍ ｋ θｐθｂ］が、Wiebe関数パラメータｍ、ｋ、θｐ、θｂの確定値として決定される（ステップ１３０）。
【００９０】
上述の通り、図６のルーチンの処理によれば、実筒内圧Ｐdataとモデル筒内圧Ｐmodelとの誤差が最小となるように、Wiebe関数パラメータｍ、ｋ、θｐ、θｂの値を決定することができる。よって、図２（Ｂ）のように実筒内圧とモデル筒内圧と間に大きな誤差が生ずることを確実に防止することができる。このため、本発明によれば、内燃機関１０の実際の筒内圧やトルクをWiebe関数モデルによって精度良く予測することができるように、Wiebe関数パラメータを決定することができる。
【００９１】
また、本実施形態では、図４のルーチンの処理によってWiebe関数パラメータの初期暫定値を実熱発生率に基づいて算出することにより、最終的な適合値（確定値）に近い、精度の良い初期暫定値を得ることができる。このため、図６のルーチンの処理によってWiebe関数パラメータの適合値を探索する上で、その精度の良い初期暫定値を出発点とすることができるので、適合値を容易に見つけ出すことができ、計算負荷を軽減することができる。
【００９２】
ここで、本発明には、以下に述べるような用途がある。図４および図６のルーチンの処理により決定されるWiebe関数パラメータの値は、内燃機関１０の運転条件、すなわち、点火時期SA、機関回転数NE、空燃比A/F、負荷率KL等に応じて変化する。本発明によれば、例えば内燃機関１０の制御ロジックの開発段階において、図４および図６に示すルーチンに従って各Wiebe関数パラメータを決定する操作を、複数の異なる運転条件に対して行うことにより、各Wiebe関数パラメータと、点火時期SA、機関回転数NE、空燃比A/F、負荷率KL等の運転条件パラメータとの関係を正確に把握することができる。そして、そのようにして把握した関係を利用すれば、あらゆる運転条件に対してWiebe関数を精度良く決定することができ、そのWiebe関数を利用して内燃機関１０の筒内圧やトルクなどを正確に予測することができる。
【００９３】
内燃機関１０の開発段階において制御ロジックの検討などを行う場合、従来は様々な運転条件の下で筒内圧やトルクなどを実測する作業が必要である。これに対し、本発明を利用して、上述したようにしてWiebe関数パラメータの値と運転条件との関係を正確に把握すれば、筒内圧やトルクの実測作業をWiebe関数に基づく予測で置き換えることができる。このため、実機データの計測点数を大幅に削減することができ、開発工数、開発コストを大幅に削減することができる。
【００９４】
また、上記のようにして把握されたWiebe関数パラメータの値と運転条件との関係は、マップや演算式の形で、Wiebe関数モデルとして車載ECUに記憶させておくことができる。そのようにすれば、車両に搭載された内燃機関１０について、種々の運転条件下における筒内圧やトルクを、Wiebe関数を用いたシミュレーションによって精度良く予測することができる。そして、その予測結果を利用すれば、例えばトルクデマンド制御やドライバビリティ改善制御などの各種制御を車載ECUが精度良く行うことが可能となる。
【００９５】
また、図４および図６のようなルーチンを、車載ECUに実行させることも可能である。図４および図６に示すようなルーチンを車載ECUによって実行することによってWiebe関数パラメータを算出（決定）し、その算出されたWiebe関数パラメータの値を現在の運転状態と関係付けて学習することとすれば、車載ECUに記憶されたWiebe関数モデルを現在の内燃機関１０の状態に応じて修整することができる。つまり、内燃機関１０の特性の経時変化や個体差を反映させて、車載ECUに記憶されたWiebe関数モデルをより高精度に修整することができる。
【００９６】
以上、本発明の実施の形態について説明したが、本発明は、上述した実施の形態に限定されるものではない。例えば、上述した実施の形態では、Wiebe関数パラメータの初期暫定値を実熱発生率に基づいて算出することとしているが、予め設定された所定の初期暫定値を出発点として各パラメータの適合値を探索するようにしてもよい。
【００９７】
また、本実施形態では、前記のような各種の数式を利用したが、本発明では、それらの数式は前記のものと完全に同じである必要はなく、それらと等価な数式を利用してもよい。また、本発明では、複数のWiebe関数パラメータのうちの少なくとも一つを図６に示すような処理で決定すればよく、残りのパラメータについては他の方法で決定するようにしてもよい。
【００９８】
また、上述した実施の形態１においては、ECU５０が、上記ステップ１００の処理を実行することにより前記第５の発明における「筒内圧計測手段」が、上記ステップ１２０および１２２の処理を実行することにより前記第５の発明における「モデル筒内圧算出手段」が、上記ステップ１２４の処理を実行することにより前記第５の発明における「誤差評価値算出手段」が、上記ステップ１２８，１２０，１２２，１２４および１２６の処理を繰り返し実行した後、上記ステップ１３０の処理を実行することにより前記第５の発明における「確定値取得手段」が、それぞれ実現されている。
【図面の簡単な説明】
【００９９】
【図１】本発明の実施の形態１のシステム構成を説明するための図である。
【図２】（Ａ）は横軸にクランク角、縦軸に熱発生率をとったグラフであり、（Ｂ）は横軸にクランク角、縦軸に筒内圧をとったグラフである。
【図３】Wiebe関数により算出されるモデル熱発生率のグラフである。
【図４】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【図５】実熱発生率と実筒内圧との関係を示す図である。
【図６】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【符号の説明】
【０１００】
１０  内燃機関
１２  クランク角センサ
１４  吸気弁
１６  排気弁
１８  筒内圧センサ
２１  吸気圧センサ
３０  点火プラグ
３３  空燃比センサ
３４  触媒
５０  ECU(Electronic Control Unit)
【図１】
0004577239.tif 000017
【図２】
0004577239.tif 000018
【図３】
0004577239.tif 000019
【図４】
0004577239.tif 000020
【図５】
0004577239.tif 000021
【図６】
0004577239.tif 000022