(57)【要約】【課題】この発明は、内燃機関の熱発生率のシミュレーション方法、内燃機関のトルクモデル作成方法および内燃機関のトルク推定方法に関し、Wiebe関数を利用して内燃機関の熱発生率、筒内圧、トルクなどを推定する場合の推定精度を向上することを目的とする。【解決手段】内燃機関の熱発生率を、複数のWiebe関数によってシミュレートする。１番目のWiebe関数は、実熱発生率の最大値およびその最大値を示すクランク角度と、１番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、決定される。ｉ番目（ｉは２以上の整数）のWiebe関数は、（ｉ−１）番目までのWiebe関数を実熱発生率から差し引いた波形の最大値およびその最大値を示すクランク角度と、ｉ番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、決定される。ｉ番目までのWiebe関数と実熱発生率との誤差を評価するための誤差評価値が基準値以下となるまで、前記ｉの値を増加させる。【選択図】図３
2008215204.tif 000002
【特許請求の範囲】
【請求項１】
内燃機関の熱発生率を、複数のWiebe関数に基づいてシミュレートする方法であって、
実熱発生率を求めるステップと、
実熱発生率の最大値およびその最大値を示すクランク角度と、１番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、１番目のWiebe関数を定めるステップと、
ｉ番目（ｉは２以上の整数）のWiebe関数について、（ｉ−１）番目までのWiebe関数を実熱発生率から差し引いた波形の最大値およびその最大値を示すクランク角度と、ｉ番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、ｉ番目のWiebe関数を定めるステップと、
ｉ番目までのWiebe関数と実熱発生率との誤差を評価するための誤差評価値が基準値以下であるか否かを判別するステップとを備え、
前記誤差評価値が前記基準値以下となるまで、前記ｉの値を増加させることを特徴とする内燃機関の熱発生率のシミュレーション方法。
【請求項２】
内燃機関の熱発生率をWiebe関数によってシミュレートする方法であって、
実熱発生率を求めるステップと、
クランク角度毎の実熱発生率とWiebe関数との偏差に対し、各クランク角度における筒内容積の逆数を重みとして乗じた値に基いて、実熱発生率と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関数について、（ｉ−１）番目までのWiebe関数を実熱発生率から差し引いた波形の最大値およびその最大値を示すクランク角度と、ｉ番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、ｉ番目のWiebe関数を定めるステップと、
ｉ番目までのWiebe関数と実熱発生率との誤差を評価するための誤差評価値が基準値以下であるか否かを判別するステップとを備え、
前記誤差評価値が前記基準値以下となるまで、前記ｉの値を増加させることを特徴とする。
【０００９】
また、第２の発明は、内燃機関の熱発生率をWiebe関数によってシミュレートする方法であって、
実熱発生率を求めるステップと、
クランク角度毎の実熱発生率とWiebe関数との偏差に対し、各クランク角度における筒内容積の逆数を重みとして乗じた値に基いて、実熱発生率とWiebe関数との誤差を評価するための誤差評価値を算出するステップと、
前記誤差評価値が最小となるようにWiebe関数を定めるステップと、
を備えることを特徴とする。
【００１０】
また、第３の発明は、内燃機関のトルクを推定するトルクモデルを作成する方法であって、
前記内燃機関の実筒内圧を運転状態毎に計測するステップと、
前記実筒内圧に基づいて、Wiebe関数によって熱発生率をシミュレートするWiebe関数モデルを作成するステップと、
前記実筒内圧に基いて実トルクを算出するステップと、
前記Wiebe関数モデルによってシミュレートされる熱発生率に基いて、モデル筒内圧を算出するステップと、
前記モデル筒内圧に基いて、モデルトルクを算出するステップと、
前記実トルクと前記モデルトルクとの差であるトルク誤差と、運転状態との関係に基いて、誤差補正モデルを作成するステップと、
を備えることを特徴とする。
【００１１】
また、第４の発明は、内燃機関のトルクを推定する方法であって、
Wiebe関数モデルを用いてモデルトルクを算出するステップと、
前記モデルトルクに与えるべき補正量を、誤差補正モデルを用いて算出するステップと、
前記モデルトルクに前記補正量を与えることにより推定トルクを算出するステップと、
を備え、
前記誤差補正モデルは、前記Wiebe関数モデルを作成する際に計測された前記内燃機関の実筒内圧から算出された実トルクと、前記モデルトルクとの差に基いて作成されたものであることを特徴とする。
【発明の効果】
【００１２】
第１の発明によれば、内燃機関の熱発生率を、複数のWiebe関数に基づいてシミュレートすることができる。その際、１番目のWiebe関数は、実熱発生率の最大値およびその最大値を示すクランク角度と、１番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、決定される。また、ｉ番目（ｉは２以上の整数）のWiebe関数は、（ｉ−１）番目までのWiebe関数を実熱発生率から差し引いた波形の最大値およびその最大値を示すクランク角度と、ｉ番目のWiebe関数の最大値およびその最大値を示すクランク角度とが一致するように、決定される。これにより、第１の発明では、複数のWiebe関数を適切な方法で組み合わせて、熱発生率を精度良くシミュレートすることができる。特に、第１の発明によれば、燃焼開始直後の熱発生率や、熱発生率の最大値付近において、誤差を極めて小さくすることができる。その領域での熱発生率は、筒内圧やトルクの推定値に大きく影響する。すなわち、第１の発明によれば、筒内圧やトルクへの影響が大きい領域の熱発生率を、特に高い精度でシミュレートすることができる。このため、筒内圧やトルクの推定精度を十分に向上することができる。更に、第１の発明によれば、ｉ番目までのWiebe関数と実熱発生率との誤差を評価するための誤差評価値が基準値以下となるまで、上記ｉの値（組み合わせるWiebe関数の数）を増やすことができる。これにより、要求される推定精度に応じて、組み合わせるWiebe関数の数を適切に選択することができる。よって、演算負荷を大きくしすぎることなく、十分なシミュレーション精度を得ることができる。
【００１３】
第２の発明によれば、内燃機関のクランク角度毎の実熱発生率とWiebe関数との偏差に対し、各クランク角度における筒内容積の逆数を重みとして乗じた値に基いて、実熱発生率とWiebe関数との誤差を評価するための誤差評価値を算出し、その誤差評価値が最小となるように、Wiebe関数を決定することができる。熱力学に基づく関係式によれば、クランク角度をθとすると、熱発生率dＱ/dθに筒内容積Ｖの逆数１/Ｖを乗じた値が、筒内圧Ｐの変化率dＰ/dθに影響する。つまり、筒内容積Ｖが小さいクランク角度領域における熱発生率が、筒内圧に大きく影響する。第２の発明によれば、そのような、筒内圧に及ぼす影響の大きいクランク角度領域における誤差が特に小さくなるように、熱発生率をシミュレートすることができる。その結果、筒内圧や、筒内圧から算出されるトルク（図示トルク）を推定する場合の推定精度を向上することができる。
【００１４】
第３の発明によれば、Wiebe関数モデルに基いて算出されるモデルトルクに対し、誤差補正モデルを用いて、補正を加えることができる。このため、Wiebe関数を利用して内燃機関のトルクを推定するトルクモデルにおいて、その推定精度を向上することができる。また、上記誤差補正モデルは、Wiebe関数モデルを作成する際に計測された実筒内圧のデータを利用して、作成することができる。つまり、誤差補正モデルの作成のために新たな計測作業が必要とされることはない。このため、開発の負担を増大することなく、トルクモデルの精度を向上することができる。
【００１５】
第４の発明によれば、Wiebe関数モデルに基いて算出されるモデルトルクに対し、誤差補正モデルを用いて、補正を加えることができる。このため、Wiebe関数を利用して内燃機関のトルクを推定するトルクモデルにおいて、その推定精度を向上することができる。また、上記誤差補正モデルは、Wiebe関数モデルを作成する際に計測された実筒内圧のデータを利用して、作成することができる。つまり、誤差補正モデルの作成のために新たな計測作業が必要とされることはない。このため、開発の負担を増大することなく、トルクモデルの精度を向上することができる。
【発明を実施するための最良の形態】
【００１６】
実施の形態１．
［システム構成の説明］
図１は、本発明の実施の形態１で用いられるシステム構成を説明するための図である。図１に示すように、本実施形態のシステムは、火花点火式の内燃機関１０を備えている。内燃機関１０の各気筒には、吸気弁１４、排気弁１６および点火プラグ３０が設けられている。
【００１７】
また、内燃機関１０には、クランク軸の回転位置（クランク角度）を検出するクランク角センサ１２が設けられている。クランク角センサ１２は、クランク軸が所定回転角だけ回転する毎に、Hi出力とLo出力を反転させるセンサである。クランク角センサ１２の出力によれば、機関回転数NEを検知することもできる。
【００１８】
また、内燃機関１０には、筒内圧センサ１８が組み込まれている。筒内圧センサ１８によれば、気筒内（燃焼室内）に生ずる圧力を検出することができる。
【００１９】
内燃機関１０の吸気通路１９の途中には、サージタンク２０が設けられている。サージタンク２０には、その内部の圧力、すなわち吸気管圧力を検出する吸気圧センサ２１が設置されている。また、吸気通路１９には、その内部を流通する吸入空気量GAを検出するエアフローメータ２２が配置されている。吸気圧センサ２１の出力によれば、内燃機関１０の負荷率KL［％］などを取得することができる。なお、負荷率KLは、吸入空気量GAと機関回転数NEとに基づいて算出してもよい。
【００２０】
エアフローメータ２２の下流には、スロットル弁２４が配置されている。スロットル弁２４の近傍には、スロットル開度TAを検出するスロットル開度センサ２６が組み付けられている。
【００２１】
内燃機関１０の吸気ポートには、ガソリン等の燃料を噴射するためのインジェクタ２８が配設されている。なお、本発明における内燃機関は、このようなポート噴射式の内燃機関に限定されるものではなく、筒内直接噴射式の内燃機関や、ポート噴射と筒内噴射とを併用する内燃機関であってもよい。
【００２２】
内燃機関１０の排気通路３２には、排気ガスの空燃比を検出する空燃比センサ３３が設置されている。また、排気通路３２には、排気ガスを浄化するための触媒３４が組み込まれている。
【００２３】
本実施形態のシステムは、制御装置として、ＥＣＵ(Electronic Control Unit)５０を備えている。ＥＣＵ５０には、上述した各種のセンサからセンサ信号が供給されている。ＥＣＵ５０は、それらのセンサ信号に基づいて、スロットル弁２４や、インジェクタ２８、点火プラグ３０などの各種アクチュエータを制御することができる。更に、ＥＣＵ５０は、クランク角センサ１２、筒内圧センサ１８等のセンサ信号に基づいて、後述するWiebe関数パラメータの決定処理を行うことができる。
【００２４】
（Wiebe関数）
本実施形態のWiebe関数は、次式で表される。
【数１】
2008215204.tif 000003
【００２５】
上記（１）式中の各記号の意味は、次の通りである。
Ｑ：内燃機関１０の気筒内の発熱量[J]
Ｑtotal：気筒内に供給された燃料の熱量[J]
θ：クランク角度[degCA]
θp：燃焼期間[degCA]
ｍ：形状パラメータ
ｋ：効率
【００２６】
上記（１）式は、内燃機関１０の気筒内の熱発生率dＱ/dθをシミュレート（模擬）するためのWiebe関数である。以下、本実施形態のWiebe関数の各パラメータについて説明する。
【００２７】
（クランク角度θ）
上記（１）式におけるクランク角度θは、気筒内の熱発生が開始した点（燃焼開始点）以降の経過クランク角度を意味する。なお、θ＝０の点の決め方については後述する。
【００２８】
（燃焼期間θp）
燃焼期間θpは、燃焼による熱の発生が継続する期間、つまり燃焼の開始から終了までの期間をクランク角度で表したものとしての物理的意味を有している。Wiebe関数は、燃焼期間θpの範囲、つまり０≦θ≦θpの範囲において定義される。
【００２９】
（形状パラメータｍ）
形状パラメータｍは、グラフ上でのWiebe関数の形状に大きく関与するパラメータである。
【００３０】
（効率ｋ）
気筒内に供給された燃料の熱量Ｑtotalは、気筒内に供給された燃料量に、その燃料の低位発熱量を乗ずることにより算出することができる。なお、低位発熱量は、真発熱量とも呼ばれる物性値である。低位発熱量とは、単位量の燃料が完全燃焼したときに発生する熱量から、燃料中に含まれる水分および燃焼によって生じる水分を蒸発させるのに必要な熱量（潜熱）を差し引いた残りの熱量を意味する。図１に示すシステムにおいては、Ｑtotalの値は、燃料の低位発熱量が既知であるものとして、インジェクタ２８からの燃料噴射量に基づいて算出することができる。あるいは、空燃比A/Fおよび筒内空気量（負荷率KL）からＱtotalを算出することもできる。
【００３１】
内燃機関１０における燃焼では、冷却損失や燃料の燃え残りなどに起因して、何らかの熱損失を伴うのが普通である。すなわち、気筒内に供給された燃料の熱量Ｑtotalが全部そのまま気筒内の発熱量Ｑに変換されることは実際上はあり得ない。本実施形態では、このことをWiebe関数に反映させるパラメータとして、効率ｋを導入している。すなわち、効率ｋは、気筒内に供給された燃料の熱量Ｑtotalが発熱量Ｑに変換される効率としての物理的意味を有しており、０＜ｋ＜１なる範囲の数である。
【００３２】
効率ｋを用いると、気筒内での総発熱量は、Ｑtotalにｋを乗じたｋＱtotalとして表すことができる。したがって、あるクランク角度θにおける燃焼割合は、そのクランク角度θまでの発熱量Ｑと、総発熱量ｋＱtotalとの比Ｑ/ｋＱtotalとして表すことができる。この燃焼割合Ｑ/ｋＱtotalは、次式で表される。
【数２】
2008215204.tif 000004
【００３３】
上記（２）式をθで微分すると、上記（１）式が得られる。つまり、上記（２）式は、上記（１）式と等価な式である。
【００３４】
（パラメータａ）
Wiebe関数における燃焼終了時の燃焼割合Ｑ/ｋＱtotalの値は、上記（２）式にθ＝θpを代入することにより、次式で表すことができる。
Ｑ/ｋＱtotal（燃焼終了時）＝１−exp（−ａ）    ・・・（３）
【００３５】
燃焼割合の定義からすると、燃焼割合の値は、燃焼終了時に１（１００％）となるべきである。しかしながら、上記（３）式の右辺で表されるように、Wiebe関数においては、燃焼終了時の燃焼割合の値が１とならない。上記（３）式の右辺をできるだけ１に近い値（ここでは０．９９９とする）にするには、ａ≒６．９とすればよい。そこで、本実施形態では、パラメータａを定数６．９として扱うものとする。ただし、本発明は、パラメータａを定数として扱う場合に限定されるものではなく、ａの値を他の方法で決定するようにしてもよい。
【００３６】
図２は、Wiebe関数および実熱発生率の一例を示すグラフである。図２中の破線で示す実熱発生率は、次のようにして算出することができる。まず、内燃機関１０の運転中に、クランク角センサ１２および筒内圧センサ１８の出力に基づいて、クランク角度毎の実筒内圧Ｐdataを計測する。実筒内圧Ｐdataを計測するクランク角度の間隔は、特に限定されず、例えば１degCA間隔とすることができる。また、実筒内圧Ｐdataの計測は、全クランク角度範囲に渡って行わなくてもよく、少なくとも実際の燃焼期間が含まれる範囲で行えばよい。
【００３７】
次に、その計測されたクランク角度毎の実筒内圧Ｐdataを下記式中のＰに代入することにより、クランク角度毎の実熱発生率dＱ/dθを算出することができる。
【数３】
2008215204.tif 000005
【００３８】
なお、上記（４）式は、熱力学上のエネルギー保存則である。同式中、Ｐ[Pa]は筒内圧、Ｖ[m3]は筒内容積、κは比熱比である。また、同式において、クランク角度θは、上死点基準の絶対的なクランク角度であるものとする。比熱比κは、燃焼ガスの組成などに基づいて定めることのできる値であり、既知である。また、筒内容積Ｖおよびその変化率dＶ/dθは、クランク角度θの関数であり、それらの関数は、内燃機関１０のクランク半径やコンロッド長などによって決定される。
【００３９】
Wiebe関数を用いて内燃機関１０の熱発生率をシミュレートする際には、前述した形状パラメータｍ、効率ｋなどのパラメータ（Wiebe関数パラメータ）の値を、Wiebe関数と実熱発生率との波形がなるべく一致するように、決定することが行われる。しかしながら、単一のWiebe関数を用いて熱発生率をシミュレートしようとした場合には、図２の破線の楕円内に示すように、Wiebe関数パラメータをどのように定めても、実熱発生率とWiebe関数との誤差を小さくすることができない場合も多い。特に、実熱発生率とWiebe関数との全体的な形状を一致させようとした場合、図２に示すように、燃焼開始直後（熱発生開始直後）や、熱発生率の最大値付近において、誤差が大きくなり易い。その結果、Wiebe関数によってシミュレートした熱発生率に基いて、内燃機関１０の筒内圧やトルクを推定した場合の推定誤差が大きくなり易い。
【００４０】
本発明者は、このような問題を解決するべく、鋭意研究を重ねた結果、複数のWiebe関数を以下に説明するような方法で組み合わせて熱発生率をシミュレートすることにより、実熱発生率との誤差を小さくすることができることを見出した。
【００４１】
図３は、本実施形態における熱発生率のシミュレーション手法を説明するための図である。図３中、グラフＡは実熱発生率であり、黒い星印は実際の点火時期SAである。図３に示す例では、３つのWiebe関数を組み合わせて熱発生率をシミュレートする場合を示している。以下、iを自然数（１，２，３，…）として、i番目のWiebe関数を、Wiebe(i)と記す。Wiebe(i)は、以下のようにして、順次定められる。
【００４２】
（Wiebe(1)の定め方）
Wiebe(1)は、以下の各条件を満足するように定める。
条件１：Wiebe(1)の起点は、点火時期SAに一致させる。すなわち、Wiebe(1)のパラメータとしてのクランク角度をθ(1)としたとき、点火時期においてθ(1)＝０となるようにする。
条件２：Wiebe(1)の最大値と、実熱発生率（すなわちグラフＡ）の最大値とを一致させる。
条件３：Wiebe(1)が最大値をとるクランク角度と、実熱発生率が最大値をとるクランク角度とを一致させる。なお、以下では、Wiebe(1)が最大値をとるクランク角度θ(1)をθ*(1)とする。
条件４：条件１〜３を満足することを前提とした上で、０≦θ(1)≦θ*(1)の範囲、すなわち、Wiebe(1)の起点からWiebe(1)および実熱発生率が最大値をとるクランク角度までの範囲における、Wiebe(1)と実熱発生率との誤差評価値（例えば、クランク角度毎の偏差の２乗和）を最小化する。
【００４３】
以上のような条件を満足するようにWiebe(1)を定めることにより、図３に示すように、点火時期SAから、Wiebe(1)および実熱発生率が最大値をとるクランク角度までの範囲において、Wiebe(1)を実熱発生率に精度良く一致させることができる。
【００４４】
その一方で、実熱発生率およびWiebe(1)が最大値をとるクランク角度以降の範囲については、実熱発生率とWiebe(1)との間に乖離が存在する。図３中のグラフＢは、その乖離の大きさを示している。つまり、このグラフＢは、実熱発生率からWiebe(1)を差し引いた波形を示している。
【００４５】
（Wiebe(2)の定め方）
Wiebe(2)は、以下の各条件を満足するように定める。
条件１：Wiebe(2)の起点は、Wiebe(1)および実熱発生率が最大値をとるクランク角度に一致させる。すなわち、Wiebe(2)のパラメータとしてのクランク角度をθ(2)としたとき、Wiebe(1)および実熱発生率が最大値をとるクランク角度においてθ(2)＝０となるようにする。
条件２：Wiebe(2)の最大値と、グラフＢ（すなわち、実熱発生率−Wiebe(1)）の最大値とを一致させる。
条件３：Wiebe(2)が最大値をとるクランク角度と、グラフＢが最大値をとるクランク角度とを一致させる。なお、以下では、Wiebe(2)が最大値をとるクランク角度θ(2)をθ*(2)とする。
条件４：条件１〜３を満足することを前提とした上で、０≦θ(2)≦θ*(2)の範囲、すなわち、Wiebe(2)の起点からWiebe(2)およびグラフＢが最大値をとるクランク角度までの範囲における、Wiebe(2)とグラフＢとの誤差評価値を最小化する。
【００４６】
図３中のグラフＣは、グラフＢと、Wiebe(2)との乖離の大きさを示している。すなわち、このグラフＣは、実熱発生率からWiebe(1)およびWiebe(2)を差し引いた波形を示している。
【００４７】
（Wiebe(3)の定め方）
Wiebe(3)は、以下の各条件を満足するように定める。
条件１：Wiebe(3)の起点は、Wiebe(2)およびグラフＢが最大値をとるクランク角度に一致させる。すなわち、Wiebe(3)のパラメータとしてのクランク角度をθ(3)としたとき、Wiebe(2)およびグラフＢが最大値をとるクランク角度においてθ(3)＝０となるようにする。
条件２：Wiebe(3)の最大値と、グラフＣ（すなわち、実熱発生率−Wiebe(1)−Wiebe(2)）の最大値とを一致させる。
条件３：Wiebe(3)が最大値をとるクランク角度と、グラフＣが最大値をとるクランク角度とを一致させる。なお、以下では、Wiebe(3)が最大値をとるクランク角度θ(3)をθ*(3)とする。
条件４：条件１〜３を満足することを前提とした上で、０≦θ(3)≦θ*(3)の範囲、すなわち、Wiebe(3)の起点からWiebe(3)およびグラフＣが最大値をとるクランク角度までの範囲における、Wiebe(3)とグラフＣとの誤差評価値を最小化する。
【００４８】
本実施形態では、上述したようにして定めたWiebe(1)〜Wiebe(3)の足し合わせによって、熱発生率をシミュレートする。
【００４９】
なお、上記の説明では、３つのWiebe関数を組み合わせる場合について述べたが、本発明において組み合わせるWiebe関数の数は、２つでも、４つ以上でもよい。４つ以上のWiebe関数を組み合わせる場合には、上記の手法と同様にして、４つ目以降のWiebe関数を順次定めればよい。
【００５０】
図４は、Wiebe(1)のみのグラフと、Wiebe(1)とWiebe(2)との足し合わせのグラフと、Wiebe(1)からWiebe(3)までの足し合わせのグラフとを、実熱発生率と比較した図である。図４から分かるように、組み合わせるWiebe関数の数を多くしていくほど、実熱発生率のグラフに近づいていく。よって、本実施形態では、必要とされるシミュレーション精度に応じて、組み合わせるWiebe関数の数を決めればよい。
【００５１】
以下、上記の手法を実現する際に必要となる数式の導出方法について説明する。まず、上記（１）式のWiebe関数が最大値をとるクランク角度θを求める。具体的には、上記（１）式をθで微分し、d2Ｑ/dθ2＝０とおいた上で、下記の式展開を行うことにより、Wiebe関数が最大値をとるクランク角度θ*を表す式が得られる。
【数４】
2008215204.tif 000006
【００５２】
上記（５）式で表されるクランク角度θ*において、上記（１）式のWiebe関数は最大値をとる。そこで、このクランク角度θを上記（１）式に代入することにより、上記（１）式のWiebe関数の最大値（dＱ/dθ）maxは、次式のように求められる。
【数５】
2008215204.tif 000007
【００５３】
また、上記（５）式より、次式が成り立つ。
【数６】
2008215204.tif 000008
【００５４】
上記（６）式中のθpに、上記（７）式の右辺を代入すると、次式が得られる。
【数７】
2008215204.tif 000009
【００５５】
なお、上記（８）式中、パラメータａは前述した通り６．９である。また、同式中、Ｑtotalは、既述した手法に従って算出可能である。
【００５６】
［実施の形態１における具体的処理］
以下、本実施形態の具体的処理について説明する。図５は、上述した手法に従ってWiebe(i)を順次定めるべく、本実施形態においてＥＣＵ５０が実行するルーチンのフローチャートである。
【００５７】
Wiebe(i)がシミュレートすべき波形をdＱ/dθ(i)とすると、図３から分かるように、dＱ/dθ(1)は、実熱発生率そのものであり、dＱ/dθ(2)は、実熱発生率からWiebe(1)を引いたものである。また、i≧３の場合には、dＱ/dθ(i)は、Wiebe(i)からWiebe(i-1)までの総和を実熱発生率から引いたものである。以上をまとめると、dＱ/dθ(i)は、次式のように表すことができる。
【数８】
2008215204.tif 000010
【００５８】
図５に示すルーチンの処理によれば、まず、上記（９）式に基いて、dＱ/dθ(i)が算出される。なお、本ルーチンの処理において、iの初期値は１とされる。また、実熱発生率の算出方法は、既述した通りである。そして、その算出されたdＱ/dθ(i)の最大値(dＱ/dθ)max(i)と、その最大値(dＱ/dθ)max(i)を示すクランク角度θ*(i)とを求める処理が実行される（ステップ１００）。
【００５９】
なお、上記クランク角度θ*(i)は、図３から分かるように、i＝１の場合には、点火時期SAから、dＱ/dθ(1)すなわち実熱発生率が最大値(dＱ/dθ)max(1)を示すクランク角度までのクランク角度幅である。また、i≧２の場合には、クランク角度θ*(i)は、dＱ/dθ(i-1)が最大値(dＱ/dθ)max(i-1)を示すクランク角度から、dＱ/dθ(i)が最大値(dＱ/dθ)max(i)を示すクランク角度までのクランク角度幅である。
【００６０】
上記ステップ１００の処理に続いて、Wiebe(i)の形状パラメータｍ（以下、ｍ(i)と記す）および効率ｋ（以下、ｋ(i)と記す）を決定する処理が実行される（ステップ１０２）。以下、この処理について、具体的に説明する。まず、Wiebe(i)のパラメータとしてのクランク角度θ（以下、θ(i)と記す）の起点、すなわちθ(i)＝０の点は、i＝１の場合には点火時期SAとされ、i≧２の場合には、dＱ/dθ(i-1)が最大値(dＱ/dθ)max(i-1)を示すクランク角度とされる。次いで、０≦θ(i)≦θ*(i)の範囲において、dＱ/dθ(i)に対するWiebe(i)の誤差が最小となるように、所定の最適化手法により、最適なｍ(i)およびｋ(i)の値が算出される。本実施形態では、上記最適化手法として、最小２乗法を用いるものとする。具体的には、各クランク角度におけるdＱ/dθ(i)とWiebe(i)との偏差の２乗を、０≦θ(i)≦θ*(i)の範囲に渡って足し合わせることにより算出される誤差評価値が最小となるように、ｍ(i)およびｋ(i)の値が決定される。また、この誤差評価値を算出するに際しては、上記（８）式の左辺に上記ステップ１００で算出された(dＱ/dθ)max(i)を代入し、同式の右辺のθ*に上記ステップ１００で算出されたθ*(i)を代入した式が満足されるように、ｍ(i)およびｋ(i)の組み合わせが選択されるものとする。これにより、dＱ/dθ(i)の最大値およびその最大値を示すクランク角度と、Wiebe(i)の最大値およびその最大値を示すクランク角度とが一致するように、ｍ(i)およびｋ(i)を決定することができる。
【００６１】
なお、Wiebe(i)の燃焼期間θp（以下、θp(i)と記す。）は、上記（７）式に、上記ステップ１００において算出されたθ*(i)と、上記ステップ１０２において決定されたｍ(i)と、ａ＝６．９とを代入することにより、算出することができる。よって、本実施形態では、各Wiebe(i)を特定するために必要なパラメータは、θ*(i)、ｍ(i)およびｋ(i)の３つとなる。
【００６２】
上記ステップ１０２の処理に続いて、実熱発生率の波形と、ΣWiebe(i)の波形との誤差が、基準以下になっているか否かが判別される（ステップ１０４）。ここで、ΣWiebe(i)は、１番目からi番目までのWiebe関数の足し合わせを意味する。上記ステップ１０４においては、具体的には、例えば上述した最小２乗法等の所定の方法に従って、実熱発生率の波形の全体に渡り、実熱発生率とΣWiebe(i)との誤差評価値が算出され、その誤差評価値が基準値以内に収まっているか否かが判別される。その結果、実熱発生率とΣWiebe(i)との誤差が基準以下になっていないと判別された場合には、iをi＝i+1と更新した上で、上記ステップ１００以下の処理が再度実行される。つまり、追加して足し合わせるためのWiebe関数を求める処理が実行される。
【００６３】
一方、上記ステップ１０４において、実熱発生率とΣWiebe(i)との誤差が基準以下になっていると判別された場合には、本ルーチンの処理が終了される。
【００６４】
以上説明したように、本実施形態によれば、複数のWiebe関数を適切な方法で組み合わせることにより、実熱発生率を精度良くシミュレートすることができる。特に、図４から分かるように、燃焼開始直後の熱発生率や、熱発生率の最大値付近において、誤差を極めて小さくすることができる。熱発生率のシミュレーション結果を利用して筒内圧やトルクを推定する場合には、燃焼開始直後の熱発生率や、熱発生率の最大値付近が、推定値に大きく影響する。このため、本実施形態によれば、筒内圧やトルクの推定精度を十分に向上することができる。
【００６５】
また、本実施形態では、組み合わせるWiebe関数の数を増やすほど、シミュレーション精度を高めることができる。一方、本実施形態では、上述したように、組み合わせるWiebe関数の各々につき、θ*、ｍおよびｋの３つのパラメータを定める必要がある。このため、組み合わせるWiebe関数の数を増やすと、演算負荷が大きくなる。本実施形態では、要求精度に応じて、組み合わせるWiebe関数の数を設定することにより、演算負荷が大きくなりすぎることなく、十分なシミュレーション精度を得ることができる。
【００６６】
以上説明したような本実施形態の熱発生率のシミュレーション方法は、例えば、内燃機関１０の制御ロジックを開発する段階などにおいて、利用することができる。従来、内燃機関１０の制御ロジックの開発に際しては、様々な運転状態の下で、筒内圧やトルクなどを実測する作業が行われている。本実施形態の熱発生率のシミュレーション方法を次のように利用すれば、このような作業を簡略化することができる。まず、内燃機関１０のいくつかの運転状態の下で、組み合わせる各Wiebe関数のパラメータθ*、ｍ、ｋの値を上述した手法に従って決定することにより、それらの値と運転状態との関係を把握する。その関係を利用すれば、各運転状態毎に、熱発生率をシミュレーションによって高精度に求めることができる。そして、その求められた熱発生率と、上記（４）式の関係とを利用すれば、クランク角度毎の筒内圧を算出することができる。そして、クランク角度毎の筒内圧が求まれば、内燃機関１０のトルク（図示トルク）を算出することができる。以上のようにして、筒内圧やトルクを実測する作業を、熱発生率シミュレーションに基づく予測で置き換えることができる。本実施形態の熱発生率のシミュレーション方法を利用すれば、実機データの計測点数を大幅に削減することができる。このため、開発工数、開発コスト、開発期間を削減することができる。
【００６７】
また、上記のようにして把握された熱発生率のシミュレーションモデルを、車載ＥＣＵに記憶させるようにしてもよい。これにより、内燃機関１０の筒内圧やトルクを車両上で精度良く予測することができ、その予測結果を利用して、例えばトルクデマンド制御やドライバビリティ改善制御などの各種制御を精度良く行うことが可能となる。
【００６８】
実施の形態２．
次に、本発明の実施の形態２について説明するが、上述した実施の形態１との相違点を中心に説明し、同様の事項については、その説明を簡略化または省略する。
【００６９】
本実施形態のWiebe関数は、次式で表される。
【数９】
2008215204.tif 000011
【００７０】
本実施形態では、上記（１０）式のWiebe関数を単独で利用して、熱発生率をシミュレートする。上記（１０）式中のクランク角度θは、点火時期SAからのクランク角度を表す。また、θbは、着火遅れである。着火遅れθbは、点火時期SAと、熱発生が開始するクランク角度との差である。本実施形態では、θ*、ｍ、ｋに加えて、着火遅れθbが、決定すべきパラメータとなる。
【００７１】
実施の形態１で説明した（４）式を変形すると、次式が得られる。
【数１０】
2008215204.tif 000012
【００７２】
上記（１１）式からは、熱発生率dＱ/dθに筒内容積Ｖの逆数１/Ｖを乗じた値が、筒内圧Ｐの変化率dＰ/dθに影響することが分かる。すなわち、筒内容積Ｖが小さいクランク角度領域における熱発生率が、筒内圧に大きく影響する。従って、筒内圧の推定を目的として熱発生率シミュレーションを行う場合であれば、筒内容積Ｖが小さいクランク角度領域において、Wiebe関数と実熱発生率とが精度良く一致していることが望ましいと言える。
【００７３】
そこで、本実施形態では、上記（１０）式のWiebe関数のパラメータを決定する際に、Wiebe関数と実熱発生率との偏差に、１/Ｖを重みとして乗じた値に基いて誤差評価値を算出することとした。具体的には、本実施形態では、誤差評価値を次式により算出する。
【数１１】
2008215204.tif 000013
【００７４】
なお、上記（１２）式中、dＱ/dθは実熱発生率であり、Wiebeは上記（１０）式のWiebe関数で算出される熱発生率である。筒内容積Ｖは、前述したようにクランク角度θの関数である。Σは、クランク角度毎の値の総和をとることを意味する。
【００７５】
本実施形態では、決定すべきWiebe関数パラメータθ*、ｍ、ｋおよびθbの値を変化させながら、上記（１２）式の誤差評価値を繰り返し算出し、誤差評価値が最小となるように、それらのWiebe関数パラメータを決定する。これにより、筒内容積Ｖの小さいクランク角度領域、すなわち筒内圧に及ぼす影響の大きいクランク角度領域における誤差が特に小さくなるように、熱発生率をシミュレートすることができる。よって、単一のWiebe関数を用いてシミュレーションを行う場合であっても、内燃機関１０の筒内圧や、筒内圧から算出されるトルク（図示トルク）を高精度に推定することができる。
【００７６】
なお、上述した実施形態では、単一のWiebe関数を用いる場合について説明したが、本発明はこれに限定されるものではない。すなわち、実施の形態１のように複数のWiebe関数を組み合わせてシミュレーションを行う場合に、本実施形態の手法を適用するようにしてもよい。
【００７７】
実施の形態３．
次に、図６を参照して、本発明の実施の形態３について説明するが、上述した実施の形態との相違点を中心に説明し、同様の事項については、その説明を簡略化または省略する。
【００７８】
本実施形態では、Wiebe関数を用いて内燃機関１０のトルクを推定する場合に、統計的な手法によって誤差を補正して、推定精度を向上させる。
【００７９】
図６は、本実施形態におけるトルクモデル作成手法を説明するための図である。図６に示すように、本実施形態では、まず、Wiebe関数を用いて内燃機関１０の熱発生率を推定するWiebe関数モデル６０を作成する。このWiebe関数モデルの作成方法は、実施の形態１で述べた手法を用いても、実施の形態２で述べた手法を用いても、あるいは他の公知の手法を用いてもよい。
【００８０】
Wiebe関数モデル６０により算出された熱発生率は、シリンダモデル６２を用いて、トルクに換算することができる。シリンダモデル６２は、まず、上記（４）式を利用して熱発生率を筒内圧（以下「モデル筒内圧」という）に換算する。次いで、そのモデル筒内圧を、次式を利用してトルクTorqueに換算する。
【数１２】
2008215204.tif 000014
【００８１】
本実施形態では、上述したようにして、Wiebe関数モデル６０およびシリンダモデル６２を用いて算出されるトルク（図示トルク）を「モデルトルク」と称する。
【００８２】
ところで、Wiebe関数モデル６０を作成する際には、前述したように、いくつかの運転状態の下で、クランク角度毎の実筒内圧を計測することが必要である。つまり、Wiebe関数モデル６０の作成がなされた場合には、筒内圧の実測データが存在していることになる。その実筒内圧を上記（１３）式に代入すれば、内燃機関１０の実際のトルク（図示トルク）を算出することができる。このようにして算出されるトルクを以下「実トルク」と称する。
【００８３】
本実施形態では、上記実トルクと上記モデルトルクとの差（トルク誤差）と、運転条件（空燃比A/F、機関回転数NE、負荷率KL、点火時期SA等）との関係に基いて、誤差補正モデル（誤差補正関数）を作成する。このような誤差補正モデルによれば、内燃機関１０の運転状態に基づいて、モデルトルクと実トルクとの間に存在する誤差を精度良く推定することができる。そこで、本実施形態では、内燃機関１０のトルクを推定する際、Wiebe関数モデル６０およびシリンダモデル６２を用いて算出されるモデルトルクを、上記誤差補正モデルによって算出されるトルク誤差で補正することによって、最終的な推定トルクを算出することとした。これにより、モデルトルクに内在する誤差を精度良くする補正することができ、正確な推定トルクを算出することができる。
【００８４】
また、本実施形態では、Wiebe関数モデルを作成する際に計測された実筒内圧のデータを利用して誤差補正モデルを作成することができる。よって、誤差補正モデルの作成のために新たな計測作業が必要とされることはない。このため、開発の負担を増大することなく、トルクモデルの精度を向上することができる。
【図面の簡単な説明】
【００８５】
【図１】本発明の実施の形態１で用いられるシステム構成を説明するための図である。
【図２】Wiebe関数および実熱発生率の一例を示すグラフである。
【図３】本発明の実施の形態１における熱発生率のシミュレーション手法を説明するための図である。
【図４】複数のWiebe関数の足し合わせのグラフと、実熱発生率とを比較した図である。
【図５】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【図６】本発明の実施の形態３におけるトルクモデル作成手法を説明するための図である。
【符号の説明】
【００８６】
１０  内燃機関
１２  クランク角センサ
１４  吸気弁
１６  排気弁
１８  筒内圧センサ
２１  吸気圧センサ
３０  点火プラグ
３３  空燃比センサ
３４  触媒
５０  ＥＣＵ(Electronic Control Unit)
【図１】
2008215204.tif 000015
【図２】
2008215204.tif 000016
【図３】
2008215204.tif 000017
【図４】
2008215204.tif 000018
【図５】
2008215204.tif 000019
【図６】
2008215204.tif 000020