(57)【要約】【課題】火炎面がピストン頂面あるいはボア壁面に到達する時期以降においても燃焼割合経過を精度良く予測することが可能な内燃機関の燃焼割合予測システムを提供する。【解決手段】１次元ＣＦＤ演算により得られた燃焼速度経過Sf,s(θ)と、ピストン速度経過Sp(θ)と、最短距離経過Lw(θ)とを用いて、火炎面到達時期θwを算出する（ステップ１０６）。点火時期θiから火炎面到達時期θwの燃焼割合経過Xb,s(θ)に対して、Wiebe関数Xb,wb(θ)を最小二乗法を用いてフィッティングして、Wiebeパラメータa,m,θb,θpを算出する（ステップ１０８）。このWiebeパラメータa,m,θb,θpをWiebe関数Xb,wb(θ)に代入することで、燃焼割合経過Xb,wb(θ)を予測する（ステップ１１０）。【選択図】図５
2008121439.tif 000002
【特許請求の範囲】
【請求項１】
内燃機関の燃焼割合の経過を予測する燃焼割合経過予測システムであって、
１次元ＣＦＤ演算により燃焼速度経過を算出する燃焼速度経過算出手段と、
１次元ＣＦＤ演算により燃焼割合経過を算出する燃焼割合経過算出手段と、
ピストンの移動速度の経過であるピストン速度経過を算出するピストン速度経過算出手段と、
点火プラグからピストン頂面までの距離と、該点火プラグからボア壁面までの距離とのうちの最短距離の経過である最短距離経過を算出する最短距離経過算出手段と、
前記燃焼速度経過と、前記ピストン速度経過と、前記最短距離経過とに基づいて、火炎面がピストン頂面又はボア壁面に到達する火炎面到達時期を算出する火炎面到達時期算出手段と、
前記火炎面到達時期以前の前記燃焼割合経過に対してＷｉｅｂｅ関数をフィッティングさせることで、Ｗｉｅｂｅパラメータを算出するＷｉｅｂｅパラメータ算出手段と、
算出されたＷｉｅｂｅパラメータをＷｉｅｂｅ関数に代入することで、燃焼割合経過を予測する燃焼割合経過予測手段とを備えたことを特徴とする内燃機関の燃焼割合経過予測システム。
【請求項２】
請求項１に記載の内燃機関の燃焼割合経過予測システムにおいて、
内燃機関の各サイクルの筒内圧経過を取得する筒内圧経過取得手段と、
前記各サイクルの筒内圧経過に基づいて、燃焼割合のサイクル間ばらつきを算出する燃焼割合ばらつき算出手段とを更に備え、
前記パラメータ算出手段は、前記燃焼割合のサイクル間ばらつきを考慮して、Ｗｉｅｂｅ関数をフィッティングさせることを特徴とする内燃機関の燃焼割合経過予測システム。
【請求項３】
内燃機関の燃焼割合の経過を予測する燃焼割合経過予測システムであって、
１次元ＣＦＤ演算により燃焼速度経過を算出する燃焼速度経過算出手段と、
１次元ＣＦＤ演算により燃焼割合経過を算出する燃焼割合経過算出手段と、
ピストンの移動速度の経過であるピストン速度経過を算出するピストン速度経過算出手段と、
点火プラグからピストン頂面までの距離と、該点火プラグからボア壁面までの距離とのうちの最短距離の経過である最短距離経過を算出する最短距離経過算出手段と、
前記燃焼速度経過と、前記ピストン速度経過と、前記最短距離経過とに基づいて、火炎面がピストン頂面又はボア壁面に到達する火炎面到達時期を算出する火炎面到達時期算出手段と、
内燃機関の各サイクルの筒内圧経過を取得する筒内圧経過取得手段と、
前記各サイクルの筒内圧経過に基づいて、燃焼割合のサイクル間ばらつきを算出する燃焼割合ばらつき算出手段と、
前記燃焼割合のサイクル間ばらつきを考慮し、かつ、前記火炎面到達時期以前よりも前記火炎面到達時期以降の重みを軽くして、点火時期から排気弁閉弁時期までの前記燃焼割合経過に対してＷｉｅｂｅ関数をフィッティングさせることで、Ｗｉｅｂｅパラメータを算出するＷｉｅｂｅパラメータ算出手段と、
算出されたＷｉｅｂｅパラメータをＷｉｅｂｅ関数に代入することで、燃焼割合経過を予測する燃焼割合経過予測手段とを備えたことを特徴とする内燃機関の燃焼割合経過予測システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、内燃機関の燃焼割合の経過を予測する燃焼割合経過予測システムに関する。
【背景技術】
【０００２】
内燃機関のシミュレーションにおいて、吸気及び排気行程での吸排気の流れを１次元及び３次元のＣＦＤ演算により予測するシステムが知られている（例えば、特許文献１参照。）。このシステムによれば、圧縮及び膨張行程では、燃焼室の作動ガスを化学反応式演算により予測し、１次元ＣＦＤ演算が併用されている。
【０００３】
【特許文献１】特開２００４−２３９１３０号公報
【特許文献２】特開２００４−２３９１２８号公報
【特許文献３】特開２００４−２３９１２９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかしながら、燃焼割合経過を予測する際に圧縮及び膨張行程で１次元ＣＦＤ演算を用いると、火炎面がピストン頂面あるいはボア壁面に到達する時期以降において、燃焼割合経過の実測値と予測値とが大きく乖離する可能性がある。
【０００５】
本発明は、上述のような課題を解決するためになされたもので、火炎面がピストン頂面あるいはボア壁面に到達する時期以降においても燃焼割合経過を精度良く予測することが可能な内燃機関の燃焼割合予測システムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
第１の発明は、上記の目的を達成するため、内燃機関の燃焼割合の経過を予測する燃焼割合経過予測システムであって、
１次元ＣＦＤ演算により燃焼速度経過を算出する燃焼速度経過算出手段と、
１次元ＣＦＤ演算により燃焼割合経過を算出する燃焼割合経過算出手段と、
ピストンの移動速度の経過であるピストン速度経過を算出するピストン速度経過算出手段と、
点火プラグからピストン頂面までの距離と、該点火プラグからボア壁面までの距離とのうちの最短距離の経過である最短距離経過を算出する最短距離経過算出手段と、
前記燃焼速度経過と、前記ピストン速度経過と、前記最短距離経過とに基づいて、火炎面がピストン頂面又はボア壁面に到達する火炎面到達時期を算出する火炎面到達時期算出手段と、
前記火炎面到達時期以前の前記燃焼割合経過に対してＷｉｅｂｅ関数をフィッティングさせることで、Ｗｉｅｂｅパラメータを算出するＷｉｅｂｅパラメータ算出手段と、
算出されたＷｉｅｂｅパラメータをＷｉｅｂｅ関数に代入することで、燃焼割合経過を予測する燃焼割合経過予測手段とを備えたことを特徴とする。
【０００７】
また、第２の発明は、第１の発明において、
内燃機関の各サイクルの筒内圧経過を取得する筒内圧経過取得手段と、
前記各サイクルの筒内圧経過に基づいて、燃焼割合のサイクル間ばらつきを算出する燃焼割合ばらつき算出手段とを更に備え、
前記パラメータ算出手段は、前記燃焼割合のサイクル間ばらつきを考慮して、Ｗｉｅｂｅ関数をフィッティングさせることを特徴とする。
【０００８】
また、第３の発明は、内燃機関の燃焼割合の経過を予測する燃焼割合経過予測システムであって、
１次元ＣＦＤ演算により燃焼速度経過を算出する燃焼速度経過算出手段と、
１次元ＣＦＤ演算により燃焼割合経過を算出する燃焼割合経過算出手段と、
ピストンの移動速度の経過であるピストン速度経過を算出するピストン速度経過算出手段と、
点火プラグからピストン頂面までの距離と、該点火プラグからボア壁面までの距離とのうちの最短距離の経過である最短距離経過を算出する最短距離経過算出手段と、
前記燃焼速度経過と、前記ピストン速度経過と、前記最短距離経過とに基づいて、火炎面がピストン頂面又はボア壁面に到達する火炎面到達時期を算出する火炎面到達時期算出手段と、
内燃機関の各サイクルの筒内圧経過を取得する筒内圧経過取得手段と、
前記各サイクルの筒内圧経過に基づいて、燃焼割合のサイクル間ばらつきを算出する燃焼割合ばらつき算出手段と、
前記燃焼割合のサイクル間ばらつきを考慮し、かつ、前記火炎面到達時期以前よりも前記火炎面到達時期以降の重みを軽くして、点火時期から排気弁閉弁時期までの前記燃焼割合経過に対してＷｉｅｂｅ関数をフィッティングさせることで、Ｗｉｅｂｅパラメータを算出するＷｉｅｂｅパラメータ算出手段と、
算出されたＷｉｅｂｅパラメータをＷｉｅｂｅ関数に代入することで、燃焼割合経過を予測する燃焼割合経過予測手段とを備えたことを特徴とする。
【発明の効果】
【０００９】
第１の発明によれば、１次元ＣＦＤ演算により算出された燃焼割合経過のうちの火炎面到達時期以前の燃焼割合経過に対して、Ｗｉｅｂｅ関数をフィッティングさせることで、Ｗｉｅｂｅパラメータが算出される。そして、この算出されたＷｉｅｂｅパラメータをＷｉｅｂｅ関数に代入することで、燃焼割合経過が予測される。このように、１次元ＣＦＤ演算により得られた燃焼割合経過と実測値との乖離が小さい火炎面到達時期以前において、Ｗｉｅｂｅ関数のフィッティングが実行される。よって、火炎面到達時期以降においても燃焼割合経過を精度良く予測することが可能となる。
【００１０】
第２の発明によれば、各サイクルの筒内圧経過に基づいて算出された燃焼割合のサイクル間ばらつきを考慮して、Ｗｉｅｂｅ関数のフィッティングが実行される。これにより、火炎面到達時期以前の燃焼割合経過に対して、Ｗｉｅｂｅ関数を精度良くフィッティングさせることができる。
【００１１】
第３の発明によれば、１次元ＣＦＤ演算により算出された燃焼割合経過に対してＷｉｅｂｅ関数のフィッティングを実行する際、各サイクルの筒内圧経過に基づいて算出された燃焼割合のサイクル間ばらつきが考慮されると共に、火炎面到達時期以前よりも火炎面到達時期以降の重みが軽くされる。これにより、火炎面到達時期以降の１次元ＣＦＤ演算結果を活用しつつ、Ｗｉｅｂｅ関数のフィッティングを精度良く実行することができる。従って、火炎面到達時期以降においても燃焼割合経過を精度良く予測することが可能となる。
【発明を実施するための最良の形態】
【００１２】
以下、図面を参照して本発明の実施の形態について説明する。尚、各図において共通する要素には、同一の符号を付して重複する説明を省略する。
【００１３】
実施の形態１．
［システム構成の説明］
図１は、本発明の実施の形態１のシステム構成を説明するための図である。図１に示すシステムは、内燃機関１として、例えば、４サイクル４気筒ガソリンエンジンを備えている。図１には、複数の気筒２のうちの１気筒のみを示している。
【００１４】
内燃機関１は、内部にピストン３を有するシリンダブロック４を備えている。ピストン３は、クランク機構を介してクランクシャフト５と接続されている。クランクシャフト５の近傍には、クランク角センサ６が設けられている。クランク角センサ６は、クランクシャフト５の回転角度（クランク角θ）を検出するように構成されている。
【００１５】
シリンダブロック４の上部にはシリンダヘッド８が組み付けられている。ピストン３上面からシリンダヘッド８までの空間は燃焼室１０を形成している。シリンダヘッド８には、燃焼室１０内の混合気に点火する点火プラグ１２が設けられている。また、シリンダヘッド８には、筒内圧を検出する筒内圧センサ１３が設けられている。
【００１６】
シリンダヘッド８は、燃焼室１０と連通する吸気ポート１４を備えている。吸気ポート１４と燃焼室１０との接続部には吸気バルブ１６が設けられている。吸気バルブ１６は、吸気バルブ１６の開弁特性を変更可能な可変動弁機構１８に接続されている。
【００１７】
吸気ポート１４には、吸気通路２０が接続されている。吸気通路２０における吸気ポート１４近傍には、該近傍に燃料を噴射するインジェクタ２２が設けられている。吸気通路２０の途中にはサージタンク２４が設けられている。サージタンク２４の上流にはスロットルバルブ２６が設けられている。スロットルバルブ２６は、スロットルモータ２８により駆動される電子制御式のバルブである。スロットルバルブ２６は、アクセル開度センサ３２により検出されるアクセル開度AAに基づいて駆動されるものである。スロットルバルブ２６の近傍には、スロットル開度TAを検出するスロットル開度センサ３０が設けられている。スロットルバルブ２６の上流には、エアフロメータ３４が設けられている。エアフロメータ３４は吸入空気量Gaを検出するように構成されている。エアフロメータ３４の上流にはエアクリーナ３６が設けられている。
【００１８】
また、シリンダヘッド８は、燃焼室１０と連通する排気ポート３８を備えている。排気ポート３８と燃焼室１０との接続部には排気バルブ４０が設けられている。排気バルブ４０は、排気バルブ４０の開弁特性を変更可能な可変動弁機構４２に接続されている。排気ポート３８には排気通路４４が接続されている。排気通路４４には、排気ガスを浄化する排気浄化触媒４６が設けられている。排気浄化触媒４６の上流には、排気空燃比を検出する空燃比センサ４８が設けられている。
【００１９】
また、本実施の形態のシステムは、制御装置としてのＥＣＵ（Electronic Control Unit）６０を備えている。ＥＣＵ６０の出力側には、点火プラグ１２、可変動弁機構１８，４２、インジェクタ２２、スロットルモータ２８等が接続されている。ＥＣＵ６０の入力側には、クランク角センサ６、筒内圧センサ１３、スロットル開度センサ３０、アクセル開度センサ３２、エアフロメータ３４、空燃比センサ４８等が接続されている。
ＥＣＵ６０は、各センサの出力に基づいて、燃料噴射制御や点火時期制御のような内燃機関全体の制御を実行する。
更に、ＥＣＵ６０は、詳細は後述するが、内燃機関１における燃焼割合の経過を予測する。
【００２０】
［実施の形態１の特徴］
従来より、内燃機関についての各種シミュレーションが実行されている。例えば、既述した特許文献１のシミュレーションによれば、圧縮及び膨張行程で、燃焼室の作動ガスを化学反応式演算により予測され、１次元ＣＦＤ演算も併用されている。
【００２１】
図２は、燃焼割合経過の実測値と、１次元ＣＦＤ演算により得られた燃焼割合経過との乖離を示す図である。図２に示すように、火炎面がピストン頂面３Ａもしくはボア壁面２Ａに到達する時期（以下「火炎面到達時期」という。）θw以降において、実測値と１次元ＣＦＤ演算結果とが大きく乖離してしまう。
【００２２】
この乖離の理由について、図３及び図４を参照して説明する。図３は、実際の火炎面の挙動を示す図である。図４は、１次元ＣＦＤ演算における火炎面の挙動を示す図である。図３及び図４において、太実線は火炎面を表している。
【００２３】
図３において矢印で示すように、実際の火炎面はピストン頂面３Ａもしくはボア壁面２Ａに到達した後、かかるピストン頂面３Ａもしくはボア壁面２Ａに沿って拡がっていく。
これに対して、１次元ＣＦＤ演算では、点火プラグ１２からピストン頂面もしくはボア壁面へ向かう軸方向のみが考慮されている。このため、図４において矢印Ａ１，Ａ２で示すように、火炎面がピストン頂面もしくはボア壁面に沿って拡がることが考慮されていない。かかる火炎面の挙動の相違により、上記のような実測値と１次元ＣＦＤ演算結果との乖離が起こる。
【００２４】
そこで、本実施の形態１では、実測値と１次元ＣＦＤ演算結果の乖離が大きい火炎面到達時期θw以降の燃焼割合経過が、火炎面到達時期θw以前の燃焼割合経過Xb,s(θ)を外挿することによって予測される。以下、図５〜図９を参照しつつ、具体的な予測手段について説明する。
【００２５】
図５は、本実施の形態１において、ＥＣＵ６０が実行するルーチンを示すフローチャートである。
図５に示すルーチンによれば、先ず、１次元ＣＦＤ演算により、燃焼速度経過Sf,s(θ)と燃焼割合経過Xb,s(θ)とをそれぞれ算出する（ステップ１００）。このステップ１００により得られた燃焼速度経過Sf,s(θ)及び燃焼割合経過Xb,s(θ)の一例を、図６に示す。図６は、１次元ＣＦＤ演算により算出された燃焼速度経過Sf,s(θ)と燃焼割合経過Xb,s(θ)を示す図である。
【００２６】
次に、ピストン３の移動速度の経過であるピストン速度経過Sp(θ)を算出する（ステップ１０２）。
ここで、このステップ１０２での処理を分かりやすくするため、クランク角θがθt（例えば６０°）からθt＋Δθ（例えば、６１°）までの区間（θt≦θ≦θt＋Δθ）において、ピストン速度経過Sp(θt)を算出する場合について説明する。
先ず、θ＝θtにおけるピストン位置Xp(θt)と、θ＝θt＋Δθにおけるピストン位置Xp(θt＋Δθ)とをそれぞれ算出する。ここで、ピストン位置Xp(θ)は、ＴＤＣ時のピストン頂面位置を基準としたピストン頂面位置であり、ハードウェア設計諸元から算出することができる。
【００２７】
次に、θ＝θtからθ＝θt＋Δθとなるまでに要する時間Δｔを算出する。そして、上記算出されたピストン位置Xp(θt),Xp(θt＋Δθ)及び時間Δｔを用いて、次式(1)によりピストン速度経過Sp(θt)を算出する。
Sp(θt)={Xp(θt)-Xp(θt+Δθ)}/Δt・・・(1)
【００２８】
続いて、点火プラグ１２からピストン頂面３Ａもしくはボア壁面２Ａまでの最短距離経過Lw(θ)を算出する（ステップ１０４）。
このステップ１０４では、先ず、図７に示される距離経過Lp(θ)、すなわち、点火プラグ１２からピストン頂面３Ａまでの距離経過Lp(θ)を算出する。図７は、点火プラグ１２からピストン頂面３Ａまでの距離経過Lp(θ)と、点火プラグ１２からボア壁面２Ａまでの最短距離Lb（後述）とを示す図である。この距離経過Lp(θ)を算出する際、ＴＤＣ時のピストン頂面３Ａと点火プラグ１２との距離Lp,0を算出する。この距離Lp,0は、ハードウェア設計諸元から算出することができる。この距離Lp,0を用いて、次式(2)により距離経過Lp(θ)を算出する。次式(2)における「Xp(θ)」は、上述したように、ＴＤＣ時のピストン頂面位置を基準としたピストン頂面位置である。
Lp(θ)=Lp,0+Xp(θ)・・・(2)
【００２９】
続いて、図７に示される最短距離Lb、すなわち、点火プラグ１２からボア壁面２Ａまでの最短距離Lbを算出する。この最短距離Lbは、ピストン位置に関わらず一定であり、ハードウェア設計諸元から算出することができる。
そして、次式(3)により最短距離経過Lw(θ)を算出する。次式(3)で表されるように、距離経過Lp(θ)と最短距離Lbのうち小さい方が、最短距離経過Lw(θ)とされる。
Lw(θ)=min{Lp(θ),Lb}・・・(3)
【００３０】
図８は、最短距離経過Lw(θ)を示す図である。図８に示すように、最短距離経過Lw(θ)は、ＴＤＣ時を含むθa≦θ≦θbの区間では距離経過Lp(θ)であり、それ以外の区間（θ<θa,θ>θb）では最短距離Lbである。
【００３１】
次に、上記ステップ１００で１次元ＣＦＤ演算により求められた燃焼速度経過Sf,s(θ)と、上記ステップ１０４で算出された最短距離経過Lw(θ)とを用いて、次式(4)により、火炎面がピストン頂面３Ａもしくはボア壁面２Ａに到達するクランク角θwを算出する（ステップ１０６）。次式(4)によれば、点火時期θi以降に火炎面が進む距離と、点火プラグ１２からの最短距離経過Lwとが等しくなる時期が求められ、その時期が火炎面到達時期θwとされる。
2008121439.tif 000003
【００３２】
次に、図９に示すように、θi≦θ≦θwの区間において、上記ステップ１００で１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)に、次式(5)で表されるWiebe関数Xb,wb(θ)を最小二乗法を用いてフィッティングすることで、Wiebeパラメータa,m,θb,θpを算出する（ステップ１０８）。図９は、Wiebe関数のフィッティングを説明するための図である。すなわち、このステップ１０８では、次式(6)により得られる値が最小となるように、Wiebeパラメータa,m,θb,θpが算出される。次式(6)における「N」は、点火時期θiから火炎面到達時期θiまでに取得されるデータ数である。
2008121439.tif 000004
2008121439.tif 000005
【００３３】
そして、上記ステップ１０８で算出されたWiebeパラメータa,m,θb,θpを上式(5)に代入することで、燃焼割合経過Xb,wb(θ)を予測する（ステップ１１０）。
【００３４】
以上説明したように、本実施の形態１によれば、１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)のうちθi≦θ≦θwの部分に対して、Wiebe関数がフィッティングされる。すなわち、実測値と１次元ＣＦＤ演算結果の乖離が大きい火炎面到達時期θw以降の燃焼割合経過が、火炎面到達時期θw以前の燃焼割合経過Xb(θ)を外挿することによって予測される。
従って、予測された燃焼割合経過Xb,wb(θ)は、θi≦θ≦θwの場合だけでなく、θ≧θwの場合においても実測値との乖離が小さい。よって、θ≧θwにおいて、１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)を用いずに、精度良く燃焼割合経過を予測することができる。その結果、筒内圧経過やトルク等を精度良く予測することが可能となる。
【００３５】
ところで、本実施の形態１においては、車両に搭載されたＥＣＵ６０により燃焼割合経過を予測しているが、該ＥＣＵ６０と通信可能なコンピュータにより予測することも可能である（後述する実施の形態２，３についても同様）。
【００３６】
尚、本実施の形態１においては、内燃機関１が第１及び第３の発明における「内燃機関」に、ピストン３が第１及び第３の発明における「ピストン」に、点火プラグ１２が第１及び第３の発明における「点火プラグ」に、筒内圧センサ１３が第２の発明における「筒内圧経過取得手段」に、それぞれ相当する。
【００３７】
また、本実施の形態１においては、ＥＣＵ６０が、ステップ１００の処理を実行することにより第１の発明における「燃焼速度経過算出手段」及び「燃焼割合経過算出手段」が、ステップ１０２の処理を実行することにより第１の発明における「ピストン速度経過算出手段」が、ステップ１０４の処理を実行することにより第１の発明における「最短距離経過算出手段」が、それぞれ実現されている。
【００３８】
また、本実施の形態１においては、ＥＣＵ６０が、ステップ１０６の処理を実行することにより第１の発明における「火炎面到達時期算出手段」が、ステップ１０８の処理を実行することにより第１の発明における「Ｗｉｅｂｅパラメータ算出手段」が、ステップ１１０の処理を実行することにより第１の発明における「燃焼割合経過予測手段」が、それぞれ実現されている。
【００３９】
実施の形態２．
次に、図１０及び図１１を参照して、本発明の実施の形態２について説明する。
本実施の形態２のシステムは、図１に示すハードウェア構成を用いて、ＥＣＵ６０に、後述する図１１に示すルーチンを実行させることにより実現することができる。
【００４０】
［実施の形態２の特徴］
上記実施の形態１では、１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)に対してWiebe関数Wb,wb(θ)をフィッティングする際、Ｎ個のデータ（すなわち、Ｎ個の「Wb,wb(θ)−Xb,s(θ)」）が等価に扱われている。
【００４１】
ところで、火炎面到達時期θw以前であっても、燃焼割合経過Xb,s(θ)のサイクル間ばらつきが生じる場合がある。このサイクル間ばらつきにより、燃焼割合結果Xb,wb(θ)の精度が低下する可能性がある。
【００４２】
そこで、本実施の形態２では、かかるサイクル間ばらつきを考慮して、Wiebe関数のフィッティングを行うようにする。
より詳細には、先ず、サイクル毎の筒内圧経過Pm,i(θ)に基づいてサイクル毎の燃焼割合経過Xb,m,i(θ)を算出し、各クランク角θにおける燃焼割合のサイクル間ばらつき（標準偏差）σxb,m(θ)を算出する。図１０は、各サイクルの筒内圧経過Pm,i(θ)と各サイクルの燃焼割合経過Xb,m,i(θ)の一例を示す図である。
そして、この燃焼割合経過Xb,m,i(θ)の標準偏差、すなわち、燃焼割合経過のサイクル間ばらつきσxb,m(θ)を求める。この求められたサイクル間ばらつきσxb,m(θ)を重み付けとして、Wiebe関数のフィッティングを行う。以下、具体的に説明する。
【００４３】
図１１は、本実施の形態２において、ＥＣＵ６０が実行するルーチンを示すフローチャートである。
図１１に示すルーチンによれば、図５に示すルーチンのステップ１００〜１０６の処理を順次実行する。
【００４４】
その後、各サイクルの筒内圧経過Pm,i(θ)から、各サイクルの燃焼割合経過Xb,m,i(θ)を算出する（ステップ１１２）。この筒内圧経過Pm,i(θ)は、筒内圧センサ１３（図１）により検出することができる。
【００４５】
このステップ１１２では、先ず、図１０に示されるような各サイクルの筒内圧経過Pm,i(θ)から、次式(7)により各サイクルの熱発生率経過dQ(θ)/dθを算出する。次式(7)は、熱力学第１法則から導出可能である。次式(7)中の「Cv」は定積比熱を、「R」は気体定数を、「V(θ)」は筒内体積経過を、それぞれ表している。筒内体積経過V(θ)は、ハードウェア設計諸元により、あるいは、上記ピストン位置Xp(θ)から算出することができる。
2008121439.tif 000006
【００４６】
次に、上式(7)により算出された各サイクルの熱発生率経過dQ(θ)/dθを用いて、次式(8)により、各サイクルの総熱発生量Qfを算出する。次式(8)中の「θEVO」は、排気バルブ４０の開弁時期（クランク角）を表している。
2008121439.tif 000007
【００４７】
最後に、上式(7)により算出された各サイクルの熱発生率経過dQ(θ)/dθと、上式(8)により算出された各サイクルの総熱発生量Qfとを用いて、次式(9)により、各サイクルの燃焼割合経過Xb,m,i(θ)を算出する。このステップ１１２で算出された各サイクルの燃焼割合経過Xb,m,i(θ)の一例を、図１０に示す。
2008121439.tif 000008
【００４８】
次に、上記ステップ１１２で算出された各サイクルの燃焼割合経過Xb,m,i(θ)を用いて、各クランク角θにおける燃焼割合のばらつきである標準偏差σxb,m(θ)を、公知の手法により算出する（ステップ１１４）。
【００４９】
次に、θi≦θ≦θwの区間において、上記ステップ１１４で算出された標準偏差σxb,m(θ)を重み付けとする最小二乗法を用いて、上記ステップ１００で１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)に、上式(5)で表されるWiebe関数Xb,wb(θ)をフィッティングすることで、Wiebeパラメータa,m,θb,θpを算出する（ステップ１１６）。すなわち、このステップ１１６では、次式(10)により得られる値が最小となるように、Wiebeパラメータa,m,θb,θpが算出される。
2008121439.tif 000009
【００５０】
ここで、上式(10)によれば、あるクランク角において燃焼割合のサイクル間ばらつきが大きい場合（すなわち、標準偏差σxb,m(θ)が大きい場合）には、そのクランク角におけるデータが過小評価されることとなる。一方、燃焼割合のサイクル間ばらつきが小さい場合（すなわち、標準偏差σxb,m(θ)が小さい場合）には、そのクランク角におけるデータが過大評価されることとなる。従って、燃焼割合のサイクル間ばらつきを考慮することで、Wiebe関数のフィッティングを精度良く行うことができる。
【００５１】
次に、上記ステップ１１６で算出されたWiebeパラメータa,m,θb,θpを上式(5)に代入することで、燃焼割合経過Xb,wb(θ)を予測する（ステップ１１０）。
【００５２】
以上説明したように、本実施の形態２によれば、各サイクルの筒内圧経過Pm,i(θ)に基づいて各サイクルの燃焼割合経過Xb,m,i(θ)を算出し、この燃焼割合経過Xb,m,i(θ)から各クランク角θにおける燃焼割合のサイクル間ばらつき（標準偏差）σxb,m(θ)が算出される。そして、この標準偏差σxb,m(θ)を重み付けとする最小二乗法を用いて、１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)に、Wiebe関数Xb,wb(θ)がフィッティングされる。すなわち、燃焼割合のサイクル間ばらつきを考慮して、Wiebe関数のフィッティングが実行される。従って、火炎面到達時期θw以前の燃焼割合経過Xb,s(θ)に対するWiebe関数Xb,wb(θ)のフィッティングを、上記実施の形態１よりも更に精度良く実行することができる。
【００５３】
ところで、本実施の形態２においては、各サイクルの燃焼割合経過Xb,m,i(θ)を算出する際に筒内圧センサ１３により検出された各サイクルの筒内圧経過Pm,i(θ)を用いているが、公知の手法により推定された各サイクルの筒内圧経過を用いてもよい。例えば、クランク角θに基づいて推定された各サイクルの筒内圧経過を用いることができる。この場合も、上記実施の形態２と同様の効果を得ることができる。後述する実施の形態３についても同様である。
【００５４】
尚、本実施の形態２においては、ＥＣＵ６０が、ステップ１１２及び１１４の処理を実行することにより第２の発明における「燃焼割合ばらつき算出手段」が、ステップ１１６の処理を実行することにより第２の発明における「パラメータ算出手段」が、それぞれ実現されている。
【００５５】
実施の形態３．
次に、図１２及び図１３を参照して、本発明の実施の形態３について説明する。
本実施の形態３のシステムは、図１に示すハードウェア構成を用いて、ＥＣＵ６０に、後述する図１３に示すルーチンを実行させることにより実現することができる。
【００５６】
［実施の形態３の特徴］
上記実施の形態１，２では、点火時期θiから火炎面到達時期θwまでの区間の燃焼割合経過Xb,s(θ)が、Wiebe関数Wb,wb(θ)のフィッティングに用いられている。
【００５７】
ところで、１次元ＣＦＤ演算により得られる燃焼割合経過Xb,s(θ)は、火炎面到達時期θwの直後から急激に精度が悪くなるわけではなく、徐々に悪くなっていく。
従って、火炎面到達時期θw以降の１次元ＣＦＤ演算結果Xb,s(θ)を有効に活用できる余地がある。
【００５８】
そこで、本実施の形態３では、Wiebe関数フィッティング時に、上記実施の形態２による標準偏差σxb,m(θ)の重み付けの概念に加えて、火炎面到達時期θw以降の重み付けがθw以前に比して軽くなるようにするように、重み付け係数k(θ)を設定する。
【００５９】
図１２は、本実施の形態３において算出される重み付け係数k(θ)の一例を示す図である。図１２に示すように、点火時期θiから火炎面到達時期θwまでは、重み付け係数k(θ)は１に算出される。さらに、火炎面到達時期θwから排気バルブ開弁時期θEVOまでは、１よりも小さな値とされ、かつ、クランク角θの増加に伴い徐々に小さくなるように算出される。
【００６０】
このように算出された重み付け係数k(θ)によれば、Wiebe関数のフィッティング時に、火炎面到達時期θw以降のデータが、火炎面到達時期θw以前のデータよりも過小評価される。従って、火炎面到達時期θw以降の１次元ＣＦＤ演算結果を活用しつつ、燃焼割合経過Xb,wb(θ)を精度良く算出することができる。
【００６１】
図１３は、本実施の形態３において、ＥＣＵ６０が実行するルーチンを示すフローチャートである。
図１３に示すルーチンによれば、図１１に示すルーチンのステップ１００〜１１４の処理を順次実行する。
【００６２】
その後、図１１に示すルーチンとは異なり、排気バルブ開弁時期θEVOまでのデータを取得する（ステップ１１８）。よって、図１１に示すルーチンでは点火時期θiから火炎面到達時期θwまでのＮ個のデータが取得されるのに対し、図１３に示すルーチンでは点火時期θiから排気バルブ開弁時期θEVOまでのＭ個のデータが取得される（図１２参照）。
【００６３】
次に、火炎面到達時期θw以前の値に比して火炎面到達時期θw以降の重み付けが軽くなるように、重み付け係数k(θ)を算出する（ステップ１２０）。このステップ１２０では、例えば、図１２に示すように、重み付け係数k(θ)が算出される。
【００６４】
次に、θi≦θ≦θEVOの区間において、上記ステップ１１４で算出された標準偏差σxb,m(θ)と、上記ステップ１２０で算出された重み付け係数k(θ)と重み付けとする最小二乗法を用いて、上記ステップ１００で１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)に、上式(5)で表されるWiebe関数Xb,wb(θ)をフィッティングすることで、Wiebeパラメータa,m,θb,θpを算出する（ステップ１２２）。すなわち、このステップ１２２では、次式(11)により得られる値が最小となるように、Wiebeパラメータa,m,θb,θpが算出される。
2008121439.tif 000010
【００６５】
ここで、上式(11)によれば、上記実施の形態２の式(10)と同様に、標準偏差σxb,m(θ)が大きいクランク角におけるデータは過小評価されると共に、標準偏差σxb,m(θ)が小さいクランク角におけるデータは過大評価される。
【００６６】
さらに、火炎面到達時期θw以前の重み付け係数k(θ)は１とされ、火炎面到達時期θw以降の重み付け係数k(θ)は１よりも小さくされる。このため、上式(11)の分母における逆数（=1/k(θ)）は、火炎面到達時期θw以前に比して火炎面到達時期θwの方が大きな値とされる。よって、火炎面到達時期θw以降のデータは、火炎面到達時期θw以前のデータよりも過小評価される。
【００６７】
また、火炎面到達時期θw以降において、クランク角θが増加するにつれて重み付け係数k(θ)は徐々に小さくされる。よって、火炎面到達時期θw以降において、クランク角θが増加するほどデータがより過小評価される。
【００６８】
次に、上記ステップ１２２で算出されたWiebeパラメータa,m,θb,θpを上式(5)に代入することで、燃焼割合経過Xb,wb(θ)を予測する（ステップ１１０）。
【００６９】
以上説明したように、本実施の形態３によれば、１次元ＣＦＤ演算により得られた燃焼割合経過Xb,s(θ)に対してWiebe関数をフィッティングする際に、各クランク角θにおける燃焼割合のサイクル間ばらつき（標準偏差）σxb,m(θ)が考慮されると共に、重み付け係数k(θ)により火炎面到達時期以前よりも火炎面到達時期以降の重みが軽くされる。これにより、火炎面到達時期以降の１次元ＣＦＤ演算結果を活用しつつ、Wiebe関数のフィッティングを精度良く実行することができる。従って、火炎面到達時期以降においても燃焼割合経過を精度良く予測することが可能となる。その結果、筒内圧経過やトルク等を精度良く予測することが可能となる。
【００７０】
尚、本実施の形態３においては、ＥＣＵ６０が、ステップ１００の処理を実行することにより第３の発明における「燃焼速度経過算出手段」及び「燃焼割合経過算出手段」が、ステップ１０２の処理を実行することにより第３の発明における「ピストン速度経過算出手段」が、ステップ１０４の処理を実行することにより第３の発明における「最短距離経過算出手段」が、ステップ１０６の処理を実行することにより第３の発明における「火炎面到達時期算出手段」が、それぞれ実現されている。
【００７１】
また、本実施の形態３においては、ＥＣＵ６０が、ステップ１１２及び１１４の処理を実行することにより第３の発明における「燃焼割合ばらつき算出手段」が、ステップ１２０及び１２２の処理を実行することにより第３の発明における「Ｗｉｅｂｅパラメータ算出手段」が、ステップ１１０の処理を実行することにより第３の発明における「燃焼割合経過予測手段」が、それぞれ実現されている。
【図面の簡単な説明】
【００７２】
【図１】本発明の実施の形態１のシステム構成を説明するための図である。
【図２】燃焼割合経過の実測値と、１次元ＣＦＤ演算により得られた燃焼割合経過との乖離を示す図である。
【図３】実際の火炎面の挙動を示す図である。
【図４】１次元ＣＦＤ演算における火炎面の挙動を示す図である。
【図５】本発明の実施の形態１において、ＥＣＵ６０が実行するルーチンを示すフローチャートである。
【図６】１次元ＣＦＤ演算により算出された燃焼速度経過Sf,s(θ)と燃焼割合経過Xb,s(θ)を示す図である。
【図７】点火プラグ１２からピストン頂面３Ａまでの距離経過Lp(θ)と、点火プラグ１２からボア壁面２Ａまでの最短距離Lbとを示す図である。
【図８】最短距離経過Lw(θ)を示す図である。
【図９】Wiebe関数のフィッティングを説明するための図である。
【図１０】各サイクルの筒内圧経過Pm,i(θ)と各サイクルの燃焼割合経過Xb,m,i(θ)の一例を示す図である。
【図１１】本発明の実施の形態２において、ＥＣＵ６０が実行するルーチンを示すフローチャートである。
【図１２】本実施の形態３において算出される重み付け係数k(θ)の一例を示す図である。
【図１３】本発明の実施の形態３において、ＥＣＵ６０が実行するルーチンを示すフローチャートである。
【符号の説明】
【００７３】
１    内燃機関
２Ａ  ボア壁面
３Ａ  ピストン頂面
６    クランク角センサ
１２  点火プラグ
１３  筒内圧センサ
４０  排気バルブ
６０  ＥＣＵ
【図１】
2008121439.tif 000011
【図２】
2008121439.tif 000012
【図３】
2008121439.tif 000013
【図４】
2008121439.tif 000014
【図５】
2008121439.tif 000015
【図６】
2008121439.tif 000016
【図７】
2008121439.tif 000017
【図８】
2008121439.tif 000018
【図９】
2008121439.tif 000019
【図１０】
2008121439.tif 000020
【図１１】
2008121439.tif 000021
【図１２】
2008121439.tif 000022
【図１３】
2008121439.tif 000023