(57)【特許請求の範囲】
【請求項１】
エンジン回転数とエンジン冷却水温度、または前記エンジン回転数とエンジン潤滑油の動粘度に基づいて、内燃機関のフリクションを推定するフリクション推定手段と、
クランク角度とクランク軸に入力する入力トルクとの関係を規定する所定の関係式に従って、前記エンジン回転数およびまたは前記クランク角度の実測値または推定値に基づいて前記入力トルクを算出する入力トルク算出手段と、
前記入力トルク算出手段により算出される前記入力トルクと、前記フリクション推定手段により算出されるフリクショントルクと、前記内燃機関の負荷トルクとに基づいて、前記内燃機関の発生トルクを算出する発生トルク算出手段と、
前記発生トルク算出手段により算出される前記発生トルクと筒内圧力との関係を規定する所定の関係式に従って、前記発生トルクに基づいて前記筒内圧力を算出する筒内圧力算出手段と、
前記筒内圧力算出手段により算出される前記筒内圧力に基づいて、前記内燃機関の各気筒の燃焼状態を表す燃焼状態指標を算出する燃焼状態指標算出手段と、
当該各気筒の前記燃焼状態指標に対して、所定サイクル数に対する燃焼状態指標平均値をそれぞれ算出する燃焼状態指標平均値算出手段と、
全気筒の前記燃焼状態指標平均値が所定の範囲内にあるか否かを判定する判定手段とを備え、
全気筒の前記燃焼状態指標平均値が前記所定の範囲内にないと判定された場合に、前記フリクション推定手段による前記フリクションの推定にずれがあると判断することを特徴とする内燃機関の制御装置。
【請求項２】
前記フリクション推定手段による前記フリクションの推定にずれがあると判断された場合に、エンジン回転数に応じて前記フリクションを学習することを特徴とする請求項１記載の内燃機関の制御装置。
【請求項３】
クランク軸の回転位置を検出するクランク角センサと、
前記筒内圧力算出手段により算出される前記筒内圧力とクランク軸周りの運動方程式とに基づいて、エンジン回転数およびまたはクランク角度の推定値を算出する演算部と、
前記演算部により算出される前記推定値を前記クランク角センサにより検出される実測値に一致させるためのフィードバック補正量を算出し、当該フィードバック補正量を前記演算部の入力値に反映させるフィードバック手段と、
を更に備えることを特徴とする請求項１または２記載の内燃機関の制御装置。
【請求項４】
前記燃焼状態指標は、筒内圧力の最大圧力値であることを特徴とする請求項１乃至３の何れか１項記載の内燃機関の制御装置。
【請求項５】
前記燃焼状態指標算出手段は、Weibe関数による燃焼圧力推定手段を含み、
前記燃焼状態指標は、当該Weibe関数の同定により算出されるWeibe関数のパラメータであることを特徴とする請求項１乃至３の何れか１項記載の内燃機関の制御装置。
【請求項６】
エンジン回転数とエンジン冷却水温度、または前記エンジン回転数とエンジン潤滑油の動粘度に基づいて、内燃機関のフリクションを推定するフリクション推定手段と、
筒内圧力と前記内燃機関の発生トルクとの関係を規定する所定の関係式に従って、前記筒内圧力に基づいて前記発生トルクを算出する発生トルク算出手段と、
前記発生トルク算出手段により算出される前記発生トルクと、前記フリクション推定手段により算出されるフリクショントルクと、前記内燃機関の負荷トルクとに基づいて、クランク軸に入力する入力トルクを算出する入力トルク算出手段と、
前記入力トルク算出手段により算出される前記入力トルクとクランク角度と前記エンジン回転数との関係を規定する所定の関係式に従って、前記エンジン回転数およびまたは前記クランク角度の推定値を算出するクランク推定値算出手段と、
クランク軸の回転位置を検出するクランク角センサと、
前記クランク推定値算出手段により算出される前記推定値を前記クランク角センサにより検出される前記エンジン回転数およびまたは前記クランク角度の実測値に一致させるためのフィードバック補正量を、前記内燃機関の各気筒の燃焼状態を表す燃焼状態指標として算出する燃焼状態指標算出手段と、
当該各気筒の前記燃焼状態指標に対して、所定サイクル数に対する燃焼状態指標平均値をそれぞれ算出する燃焼状態指標平均値算出手段と、
全気筒の前記燃焼状態指標平均値が所定の範囲内にあるか否かを判定する判定手段とを備え、
全気筒の前記燃焼状態指標平均値が前記所定の範囲内にないと判定された場合に、前記フリクション推定手段による前記フリクションの推定にずれがあると判断することを特徴とする内燃機関の制御装置。
【請求項７】
前記フリクション推定手段による前記フリクションの推定にずれがあると判断された場合に、エンジン回転数に応じて前記フリクションを学習することを特徴とする請求項６記載の内燃機関の制御装置。
【請求項８】
前記フィードバック補正量を前記筒内圧力に反映させるフィードバック手段を更に備えることを特徴とする請求項６または７記載の内燃機関の制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、内燃機関の制御装置に関する。
【背景技術】
【０００２】
従来、例えば特許文献１には、内燃機関の失火検出装置が開示されている。この従来の失火検出装置は、気筒間の発生トルクの差および気筒間のフリクショントルクの差を学習する手段を備え、それらの学習値に応じて、クランク軸の角速度（クランク角回転速度）の計測値を補正するというものである。
【０００３】
【特許文献１】特開平４−１０１０７１号公報
【特許文献２】特公平３−２０５９０号公報
【特許文献３】特開２００２−７１０００号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
内燃機関を高精度に制御するためには、燃焼状態を正確に推定する手段が必要とされる。そして、燃焼状態を正確に推定するうえで、フリクションがもたらす影響は大きなものである。上記従来の装置によれば、気筒間のフリクションの差は学習されるが、環境変化や経年変化に起因して生ずるフリクションのずれが考慮されていない。このため、上記従来の装置は、燃焼状態を正確に推定する制御ロジックを構築するうえで、未だ改善の余地を有するものであった。
【０００５】
この発明は、上述のような課題を解決するためになされたもので、経年変化等に依らずに正確なフリクション情報を得ることができ、その結果として、高精度なエンジン制御を実現し得る内燃機関の制御装置を提供することを目的とする。
【課題を解決するための手段】
【０００６】
第１の発明は、エンジン回転数とエンジン冷却水温度、または前記エンジン回転数とエンジン潤滑油の動粘度に基づいて、内燃機関のフリクションを推定するフリクション推定手段と、
クランク角度とクランク軸に入力する入力トルクとの関係を規定する所定の関係式に従って、前記エンジン回転数およびまたは前記クランク角度の実測値または推定値に基づいて前記入力トルクを算出する入力トルク算出手段と、
前記入力トルク算出手段により算出される前記入力トルクと、前記フリクション推定手段により算出されるフリクショントルクと、前記内燃機関の負荷トルクとに基づいて、前記内燃機関の発生トルクを算出する発生トルク算出手段と、
前記発生トルク算出手段により算出される前記発生トルクと筒内圧力との関係を規定する所定の関係式に従って、前記発生トルクに基づいて前記筒内圧力を算出する筒内圧力算出手段と、
前記筒内圧力算出手段により算出される前記筒内圧力に基づいて、前記内燃機関の各気筒の燃焼状態を表す燃焼状態指標を算出する燃焼状態指標算出手段と、
当該各気筒の前記燃焼状態指標に対して、所定サイクル数に対する燃焼状態指標平均値をそれぞれ算出する燃焼状態指標平均値算出手段と、
全気筒の前記燃焼状態指標平均値が所定の範囲内にあるか否かを判定する判定手段とを備え、
全気筒の前記燃焼状態指標平均値が前記所定の範囲内にないと判定された場合に、前記フリクション推定手段による前記フリクションの推定にずれがあると判断することを特徴とする。
【０００７】
また、第２の発明は、第１の発明において、前記フリクション推定手段による前記フリクションの推定にずれがあると判断された場合に、エンジン回転数に応じて前記フリクションを学習することを特徴とする。
【０００８】
また、第３の発明は、第１または第２の発明において、クランク軸の回転位置を検出するクランク角センサと、
前記筒内圧力算出手段により算出される前記筒内圧力とクランク軸周りの運動方程式とに基づいて、エンジン回転数およびまたはクランク角度の推定値を算出する演算部と、
前記演算部により算出される前記推定値を前記クランク角センサにより検出される実測値に一致させるためのフィードバック補正量を算出し、当該フィードバック補正量を前記演算部の入力値に反映させるフィードバック手段と、
を更に備えることを特徴とする。
【０００９】
また、第４の発明は、第１乃至第３の発明の何れかにおいて、前記燃焼状態指標は、筒内圧力の最大圧力値であることを特徴とする。
【００１０】
また、第５の発明は、第１乃至第３の発明の何れかにおいて、前記燃焼状態指標算出手段は、Weibe関数による燃焼圧力推定手段を含み、
前記燃焼状態指標は、当該Weibe関数の同定により算出されるWeibe関数のパラメータであることを特徴とする。
【００１１】
また、第６の発明は、エンジン回転数とエンジン冷却水温度、または前記エンジン回転数とエンジン潤滑油の動粘度に基づいて、内燃機関のフリクションを推定するフリクション推定手段と、
筒内圧力と前記内燃機関の発生トルクとの関係を規定する所定の関係式に従って、前記筒内圧力に基づいて前記発生トルクを算出する発生トルク算出手段と、
前記発生トルク算出手段により算出される前記発生トルクと、前記フリクション推定手段により算出されるフリクショントルクと、前記内燃機関の負荷トルクとに基づいて、クランク軸に入力する入力トルクを算出する入力トルク算出手段と、
前記入力トルク算出手段により算出される前記入力トルクとクランク角度と前記エンジン回転数との関係を規定する所定の関係式に従って、前記エンジン回転数およびまたは前記クランク角度の推定値を算出するクランク推定値算出手段と、
クランク軸の回転位置を検出するクランク角センサと、
前記クランク推定値算出手段により算出される前記推定値を前記クランク角センサにより検出される前記エンジン回転数およびまたは前記クランク角度の実測値に一致させるためのフィードバック補正量を、前記内燃機関の各気筒の燃焼状態を表す燃焼状態指標として算出する燃焼状態指標算出手段と、
当該各気筒の前記燃焼状態指標に対して、所定サイクル数に対する燃焼状態指標平均値をそれぞれ算出する燃焼状態指標平均値算出手段と、
全気筒の前記燃焼状態指標平均値が所定の範囲内にあるか否かを判定する判定手段とを備え、
全気筒の前記燃焼状態指標平均値が前記所定の範囲内にないと判定された場合に、前記フリクション推定手段による前記フリクションの推定にずれがあると判断することを特徴とする。
また、第７の発明は、第６の発明において、前記フリクション推定手段による前記フリクションの推定にずれがあると判断された場合に、エンジン回転数に応じて前記フリクションを学習することを特徴とする。
また、第８の発明は、第６または第７の発明において、前記フィードバック補正量を前記筒内圧力に反映させるフィードバック手段を更に備えることを特徴とする。
【発明の効果】
【００１２】
第１の発明によれば、燃焼状態指標算出手段によって、内燃機関の燃焼状態を正確に取得しているので、フリクション推定手段によるフリクションの推定にずれが生じているか否かを的確に判断することができる。
【００１３】
第２の発明によれば、フリクションの推定にずれがあると認められる場合に、フリクションのずれが学習されるので、経年変化等によるフリクションのずれを補正することができる。このため、より正確なエンジン制御を行うことが可能となる。
【００１４】
第３の発明によれば、内燃機関の出力軸であるクランク軸周りの運動方程式に基づいて、軸トルク情報を用いてクランク軸に作用するフリクションを正確に見積もったうえで、そのフリクションのずれを検証することができる。このため、各気筒の燃焼状態のばらつきに占めるフリクションの影響をより明確にすることができる。
【００１５】
第４の発明によれば、筒内圧力の最大圧力値を燃焼状態指標として、燃焼状態を正確に取得することができる。
【００１６】
第５の発明によれば、Weibe関数のパラメータを燃焼状態指標として、燃焼状態を正確に取得することができる。
【００１７】
第６の発明によれば、フィードバック補正量を燃焼状態指標として用いる燃焼状態指標算出手段によって、内燃機関の燃焼状態を正確に取得しているので、フリクション推定手段によるフリクションの推定にずれが生じているか否かを的確に判断することができる。
第７の発明によれば、フリクションの推定にずれがあると認められる場合に、フリクションのずれが学習されるので、経年変化等によるフリクションのずれを補正することができる。このため、より正確なエンジン制御を行うことが可能となる。
第８の発明によれば、軸トルク情報を用いてクランク軸に作用するフリクションを正確に見積もったうえで、そのフリクションのずれを検証することができる。このため、各気筒の燃焼状態のばらつきに占めるフリクションの影響をより明確にすることができる。
【発明を実施するための最良の形態】
【００１８】
実施の形態１．
［実施の形態１の装置の構成］
図１は、本発明の実施の形態１の内燃機関の停止位置制御装置が適用された内燃機関１０の構成を説明するための図である。本実施形態のシステムは、内燃機関１０を備えている。ここでは、内燃機関１０は、直列４気筒型エンジンであるものとする。内燃機関１０の筒内には、ピストン１２が設けられている。ピストン１２は、コンロッド１４を介してクランク軸１６と連結されている。また、内燃機関１０の筒内には、ピストン１２の頂部側に燃焼室１８が形成されている。燃焼室１８には、吸気通路２０および排気通路２２が連通している。
【００１９】
吸気通路２０には、スロットルバルブ２４が設けられている。スロットルバルブ２４は、アクセル開度と独立してスロットル開度を制御することのできる電子制御式スロットルバルブである。スロットルバルブ２４の近傍には、スロットル開度TAを検出するスロットルポジションセンサ２６が配置されている。スロットルバルブ２４の下流には、内燃機関１０の吸気ポートに燃料を噴射するための燃料噴射弁２８が配置されている。また、内燃機関が備えるシリンダヘッドには、気筒毎に、燃焼室１８の頂部から燃焼室１８内に突出するように点火プラグ３０がそれぞれ取り付けられている。吸気ポートおよび排気ポートには、それぞれ、燃焼室１８と吸気通路２０、或いは燃焼室１８と排気通路２２を導通状態または遮断状態とするための吸気弁３２および排気弁３４が設けられている。
【００２０】
吸気弁３２および排気弁３４は、それぞれ吸気可変動弁（VVT）機構３６および排気可変動弁（VVT）機構３８により駆動される。可変動弁機構３６、３８は、それぞれ、クランク軸の回転と同期して吸気弁３２および排気弁３４を開閉させると共に、それらの開弁特性（開弁時期、作用角、リフト量など）を変更することができる。
【００２１】
内燃機関１０は、クランク軸の近傍にクランク角センサ４０を備えている。クランク角センサ４０は、クランク軸が所定回転角だけ回転する毎に、Hi出力とLo出力を反転させるセンサである。クランク角センサ４０の出力によれば、クランク軸の回転位置やその回転速度（エンジン回転数Ne）を検知することができる。また、内燃機関１０は、吸気カム軸の近傍にカム角センサ４２を備えている。カム角センサ４２は、クランク角センサ４０と同様の構成を有するセンサである。カム角センサ４２の出力によれば、吸気カム軸の回転位置（進角量）などを検知することができる。
【００２２】
図１に示すシステムは、ECU(Electronic Control Unit)５０を備えている。ECU５０には、上述した各種センサに加え、排気通路２２内の排気空燃比を検出するための空燃比センサ５２や内燃機関１０の冷却水温度を検出するための水温センサ５４が接続されている。また、ECU５０には、上述した各種アクチュエータが接続されている。ECU５０は、それらのセンサ出力、およびECU５０内に仮想的に構成されたエンジンモデル６０を用いた演算結果に基づいて、内燃機関１０の運転状態を制御することができる。
【００２３】
［エンジンモデルの概要］
図２は、図１に示すECU５０が備えるエンジンモデル６０の構成を示すブロック図である。図２に示すように、エンジンモデル６０は、クランク軸周りの運動方程式演算部６２と、フリクションモデル６４と、吸気圧力推定モデル６６と、筒内圧推定モデル６８と、燃焼波形算出部７０と、大気圧補正項算出部７２と、大気温補正項算出部７４とを含んでいる。以下、これらの各部の詳細な構成について説明を行う。
【００２４】
（１）クランク軸周りの運動方程式演算部について
クランク軸周りの運動方程式演算部６２は、クランク角度θおよびエンジン回転数Ne（クランク角回転速度dθ/dt）のそれぞれの推定値を求めるためのものである。クランク軸周りの運動方程式演算部６２は、筒内圧推定モデル６８または燃焼波形算出部７０から内燃機関１０の筒内圧力Pの入力を受け、演算開始時には、更に、初期クランク角度θ0および初期エンジン回転数Ne0の入力を受ける。
【００２５】
クランク軸周りの運動方程式演算部６２によって算出される推定クランク角度θおよび推定エンジン回転数Neは、図２に示すPIDコントローラ７６によって、実クランク角度θおよび実エンジン回転数Neとの偏差が無くなるようにフィードバック制御される。また、クランク軸周りの運動方程式演算部６２の演算結果には、フリクションモデル６４によって、内燃機関１０の内部のフリクションに関する影響が反映される。
【００２６】
次に、クランク軸周りの運動方程式演算部６２の内部で実行される具体的な演算内容について説明する。
図３は、クランク軸周りの各要素に付す記号を示す図である。図３に示すように、ここでは、筒内圧力Pを受けるピストン１２の頂部の表面積をAとする。コンロッド１４の長さをL、クランクの回転半径をrとする。そして、コンロッド１４のピストン取り付け点とクランク軸１６の軸中心とを結ぶ仮想線（シリンダの軸線）と、コンロッド１４の軸線とがなす角度をφ（以下、「コンロッド角度φ」と称する）とし、シリンダの軸線とクランクピン１７の軸線とがなす角度をθとする。
【００２７】
４つの気筒を有する内燃機関１０では、気筒間のクランク角度の位相差は180°CAであるため、それらの気筒間のクランク角度の関係は、次の（１ａ）式のように定義することができる。また、各気筒のクランク角回転速度dθ/dtは、それぞれ各気筒のクランク角度θの時間微分となるため、それぞれ次の（１ｂ）式のように表すことができる。
【数１】
0004665788.tif 000002
【００２８】
ただし、上記（１ａ）式および（１ｂ）式において、クランク角度θおよびクランク角回転速度dθ/dtに付された符号1〜4は、内燃機関１０の所定の爆発順序に従って燃焼が到来する気筒の順番に対応しており、また、後述する数式においては、それらの符号1〜4を「i」で代表させることがある。
【００２９】
また、図３に示すピストン・クランク機構においては、クランク角度θiとコンロッド角度φiとは、次の（２）式で表される関係を有することになる。
【数２】
0004665788.tif 000003
ただし、上記（２）式において、dXi/dtはピストン速度である。
【００３０】
また、クランク軸周りの全運動エネルギTは、次の（３）式のように表すことができる。（３）式を展開すると、（３）式中の各項の諸々のパラメータを1/2(dθ/dt)2の係数としてまとめることができる。ここでは、そのようにまとめられた係数を、クランク角度θの関数f(θ)として表現している。
【数３】
0004665788.tif 000004
【００３１】
ただし、上記（３）式において、右辺第１項はクランク軸１６の回転運動に関する運動エネルギに、右辺第２項はピストン１２およびコンロッド１４の直進運動に関する運動エネルギに、右辺第３項はコンロッド１４の回転運動に関する運動エネルギに、それぞれ対応している。また、上記（３）式において、Ikはクランク軸１６の軸周りの慣性モーメントであり、Iflはフライホイールの回転軸周りの慣性モーメントであり、Imiは内燃機関１０と組み合わされる変速機以下の回転部（すなわち、変速機、駆動軸、タイヤ等）の回転軸周りの慣性モーメントであり、Icはコンロッドに関する慣性モーメントである。また、mpはピストン１２の変位であり、mcはコンロッド１４の変位である。
【００３２】
次に、ラグラジアンLを、系の全運動エネルギTと位置エネルギUとの偏差として、次の（４ａ）式のように定義する。そして、クランク軸１６に作用する入力トルクをTRQとすると、ラグランジュの運動方程式を用いて、ラグラジアンLとクランク角度θと入力トルクTRQとの関係を、次の（４ｂ）式のように表すことができる。
【数４】
0004665788.tif 000005
【００３３】
ここで、上記（４ａ）式において、位置エネルギUの影響は運動エネルギTの影響に比して小さく、その影響を無視することができる。従って、上記（４ｂ）式の左辺第１項は、上記（３）式をクランク角回転速度(dθ/dt)で偏微分して得られた値を時間微分することで、クランク角度θの関数として、次の（４ｃ）式のように表すことができる。また、上記（４ｂ）式の左辺第２項は、上記（３）式をクランク角度θで偏微分することで、クランク角度θの関数として、次の（４ｄ）式のように表すことができる。
【００３４】
従って、上記（４ｂ）式は、次の（４ｅ）式のようにして表すことができ、これにより、クランク角度θと入力トルクTRQとの関係を得ることができる。また、ここでは、その入力トルクTRQを、次の（５）式のように、３つのパラメータからなるものと定義する。
【数５】
0004665788.tif 000006
【００３５】
ただし、上記（５）式において、TRQeは、エンジン発生トルクであり、より具体的には、ガス圧力（筒内圧力P）を受けるピストン１２からクランク軸１６に作用するトルクである。TRQLは、負荷トルクであり、内燃機関１０が搭載される車両の特性に応じて異なる既知の値として、ECU５０に記憶されている。TRQfは、フリクショントルク、すなわち、ピストン１２、クランク軸１６等の摺動部分の摩擦損失に対応するトルクである。このフリクショントルクTRQfは、フリクションモデル６４から得られる値である。
【００３６】
次に、エンジン発生トルクTRQeは、次の（６ａ）式〜（６ｃ）式に従って算出することができる。すなわち、先ず、筒内圧力Pに基づいてコンロッド１４に作用する力Fcは、ピストン１２の頂部に作用する力PAのコンロッド１４の軸線方向成分として、（６ａ）式のように表すことができる。そして、図３に示すようにコンロッド１４の軸線とクランクピン１７の軌跡の接線とがなす角度αが{π/2−(φ＋θ)}であるため、筒内圧力Pに基づいてクランクピン１７の軌跡の接線方向に作用する力Fkは、コンロッド１４に作用する力Fcを用いて、（６ｂ）式のように表すことができる。従って、エンジン発生トルクTRQeは、クランクピン１７の軌跡の接線方向に作用する力Fkとクランクの回転半径rとの積であるため、（６ａ）式および（６ｂ）式を用いて、（６ｃ）式のように表すことができる。
【数６】
0004665788.tif 000007
【００３７】
以上説明したクランク軸周りの運動方程式演算部６２の構成によれば、筒内圧推定モデル６８または燃焼波形算出部７０によって筒内圧力Pを取得することにより、（６ｃ）式および（５）式に従って入力トルクTRQを得ることができる。そして、（４ｅ）式を解くことにより、クランク角度θやクランク角回転速度dθ/dtを得ることが可能となる。
【００３８】
（２）フリクションモデルについて
図４は、図２に示すフリクションモデル６４がフリクショントルクTRQfを取得するために備えているフリクションマップの一例を示している。図４に示すフリクションマップでは、フリクショントルクTRQfを、エンジン回転数Neとエンジン冷却水温度との関係で定めている。このようなフリクションマップの特性は、予め実験等により定められたものであり、フリクショントルクTRQfは、エンジン冷却水温度が低くなると大きくなる傾向を有している。
【００３９】
尚、ここでは、ECU５０の計算負荷の低減のため、フリクションモデル６４として、上記のようなフリクションマップを備えるようにしているが、フリクションモデルの構成は、これに限定されるものではなく、以下の（７）式のような関係式を用いるものであってもよい。この（７）式では、フリクショントルクTRQfが、エンジン回転数Neと内燃機関１０の潤滑油の動粘度νとをパラメータとする関数となるように構成されている。
【数７】
0004665788.tif 000008
ただし、上記（７）式において、C1、C2、C3は、それぞれ実験等により適合される係数である。
【００４０】
（３）吸気圧力推定モデルについて
吸気圧力推定モデル６６は、吸気圧力を推定するための吸気圧マップ（図示省略）を備えている。この吸気圧マップは、吸気圧力を、スロットル開度TA、エンジン回転数Ne、および吸排気弁のバルブタイミングVVTとの関係で定めたものである。このような吸気圧力推定モデルの構成によれば、ECU５０の計算負荷を低く抑えつつ、吸気圧力を取得することができる。尚、詳細に吸気圧力を計算する場合には、上記のような吸気圧マップを用いずに、スロットルバルブ２４を通過する空気流量を推定するスロットルモデルと、吸気弁３２の周囲を通過する空気流量（すなわち、筒内吸入空気流量）を推定するバルブモデルとを用いて、吸気圧力推定モデルを構成するようにしてもよい。
【００４１】
（４）筒内圧推定モデル
筒内圧推定モデル６８は、燃焼が行われない状況下で、筒内圧力Pを算出するために用いられるモデルである。この筒内圧推定モデル６８では、内燃機関１０の各行程における筒内圧力Pを、次の（８ａ）式〜（８ｄ）式を用いて算出するようにしている。すなわち、先ず、吸気行程の経過中の筒内圧力Pは、（８ａ）式で示すように、上述した吸気圧力推定モデル６６が有する吸気圧マップから得られる筒内圧力のマップ値Pmapから得るようにしている。
【数８】
0004665788.tif 000009
【００４２】
次に、圧縮行程の経過中の筒内圧力Pは、気体の可逆断熱変化の式に基づいて、（８ｂ）式のように表すことができる。
ただし、上記（８ｂ）式において、VBDCはピストン１２が吸気下死点にあるときの行程容積Vであり、κは比熱比である。
【００４３】
また、膨張行程の経過中の筒内圧力Pについても、圧縮行程の場合と同様にして、（８ｃ）式のように表すことができる。
ただし、上記（８ｃ）式において、VTDCはピストン１２が圧縮上死点にあるときの行程容積Vであり、Pcは圧縮行程の終了時における筒内圧力である。
【００４４】
また、排気行程の経過中の筒内圧力Pは、(８ｄ)式で示すように、排気通路２２内の圧力Pexであるものとしている。この圧力Pexは、ほぼ大気圧力Pairに等しいとみなすことができるものである。従って、ここでは、大気圧力Pairを、排気行程の経過中の筒内圧力Pに使用している。
【００４５】
（５）燃焼波形算出部について
燃焼波形算出部７０は、圧縮行程の途中から膨張行程の途中までの燃焼が行われている期間における筒内圧力（燃焼圧力）Pを算出するために用いられるモデルである。この燃焼波形算出部７０では、Weibe関数を用いた関係式である（９ａ）式と、後述する（１０）式とを用いて、燃焼圧力Pの推定値が算出される。
【数９】
0004665788.tif 000010
【００４６】
より具体的には、燃焼波形算出部７０では、先ず、（９ａ）式を用いて、現在のクランク角度θに対応する熱発生率dQ/dθを算出することとしている。
ただし、上記（９ａ）式において、mは形状係数、kは燃焼効率、θbは着火遅れ期間、aは燃焼速度（ここでは固定値6.9）である。これらの各パラメータは、事前に適合された値が使用される。また、Qは発熱量である。
【００４７】
上記（９ａ）式を用いて熱発生率dQ/dθを算出するには、発熱量Qを算出する必要がある。発熱量Qは、微分方程式である（９ａ）式を解くことにより算出することができる。そのために、先ず、（９ｂ）式では、（９ａ）式におけるWeibe関数に相当する部分をg(θ)と置き換えている。そうすると、（９ａ）式を（９ｃ）式のように表すことが可能となる。次いで、（９ｃ）式の両辺をクランク角度θで積分した後に、当該（９ｃ）式を展開することで、発熱量Qを（９ｄ）式のように表すことができる。次いで、（９ｄ）式に従って算出された発熱量Qを、再度（９ａ）式に代入することで、熱発生率dQ/dθが算出される。
【００４８】
熱発生率dQ/dθと筒内圧力（燃焼圧力）Pとは、エネルギ保存則に基づく関係式を用いて（１０）式のように表すことができる。従って、（９ａ）式に従って算出された熱発生率dQ/dθを代入して当該（１０）式を解くことにより、燃焼圧力Pを算出することができる。
【数１０】
0004665788.tif 000011
【００４９】
以上説明した筒内圧推定モデル６８および燃焼波形算出部７０によれば、筒内圧推定モデル６８を用いて燃焼が行われていない状況下での筒内圧力Pを算出するととともに、燃焼波形算出部７０を用いて燃焼が行われている期間中の筒内圧力Pを算出することにより、燃焼実行の有無に関係なく、内燃機関１０の筒内圧力Pの履歴を取得することができる。
【００５０】
尚、内燃機関１０の筒内圧力Pの履歴を取得する手法は、上記の手法に限定されるものではなく、例えば、以下の図５を参照して示すような手法であってもよい。
図５は、そのような変形例の手法を説明するための図である。この手法では、上記（９ａ）式および（１０）式を用いて、所定のクランク角度θ毎に燃焼圧力Pを計算することを行うのではなく、事前に、上記（９ａ）式および（１０）式を用いて、図５（Ａ）に示すような燃焼パターン、すなわち、燃焼に付されることで変化する筒内圧力Pの波形の変化分（燃焼による圧力増加分）のみを算出しておく。
【００５１】
より具体的には、そのような燃焼パターンを決定する３つのパラメータである着火遅れ期間、燃焼期間、およびΔPmax（燃焼時の最大圧力Pmaxと燃焼無し時の最大圧力Pmax0との偏差）を、エンジン回転数Ne、空気充填率KL、吸排気弁のバルブタイミングVVT、および点火時期のそれぞれとの関係で定めたマップを記憶しておく。そして、燃焼による圧力増加分に対応する波形を、２次関数などの簡易な関数を組み合わせて近似させた波形として算出するために、当該近似波形の各係数を上記のエンジン回転数Neとの関係でマップ化しておく。そして、図５（Ｂ）に示すように、そのようなマップを参照して得られた燃焼による圧力増加分の波形を、筒内圧推定モデル６８で算出される筒内圧力Pの値と足し合わせることで、燃焼圧力Pを取得するようにする。
【００５２】
（６）大気圧補正項算出部について
大気圧補正項算出部７２は、筒内に吸入される筒内充填空気量Mcを推定するモデル（ここでは「エアモデル」と称する）を含んでいる。このエアモデルでは、筒内充填空気量Mcを次の（１１）式に従って算出することとしている。
【数１１】
0004665788.tif 000012
【００５３】
ただし、上記（１１）式において、a、bは、それぞれ運転条件（エンジン回転数NeやバルブタイミングVVTなど）に応じて適合される係数である。尚、Pmは、吸気圧力であり、例えば、上述した吸気圧力推定モデル６６によって算出される値を使用することができる。
【００５４】
また、大気圧補正項算出部７２は、筒内に吸入される燃料量fcを推定するモデル（ここでは「燃料モデル」と称する）を含んでいる。燃料噴射弁２８から噴射された後の燃料の挙動を考慮すると、すなわち、噴射された燃料の一部の吸気ポートの内壁等への付着やその付着燃料の気化という現象を考慮すると、第kサイクルにおける燃料噴射の開始時における壁面付着燃料量がfw(k)であり、第kサイクルにおける実燃料噴射量がfi(k)である場合、第kサイクルの終了後に発生している壁面付着燃料量fw(k+1)、および第kサイクルにおいて筒内に吸入される燃料量fcは、次の（１２a）式および（１２ｂ）式のように表すことができる。
【数１２】
0004665788.tif 000013
【００５５】
ただし、上記（１２）式において、Pは、付着率、より具体的には、噴射燃料量fiのうちの吸気ポートの内壁等に付着する燃料量の割合である。Rは、残留率、より具体的には、吸気行程の実行後に付着燃料量fwが壁面等に付着したままの状態で残る割合である。
上記（１２）式によれば、付着率Pおよび残留率Rをパラメータとして、個々のサイクル毎に上記燃料量fcを算出することができる。
【００５６】
従って、上記のエアモデルおよび燃料モデルの算出結果を用いて、空燃比A/Fの推定値を算出することができる。大気圧補正項算出部７２では、次いで、この推定空燃比A/Fと、噴射された燃料が燃焼に付された後に空燃比センサ５２に到達するまでの輸送遅れを考慮したタイミングで検出する空燃比A/Fの実測値との定常偏差を算出する。そして、この定常偏差が筒内充填空気量Mcの誤差であるため、当該定常偏差が大きい場合には、大気圧がずれているものとして、大気圧補正係数kairpを算出する。具体的には、上記エアモデルより吸気圧力Pmを逆算し、その吸気圧力Pmに基づいて標準大気圧Pa0に対する補正率として大気圧補正係数kairpを算出する。この大気圧補正係数kairpは、上述した吸気圧力推定モデル６６および筒内圧推定モデル６８において、吸気圧力Pmapと排気圧力（大気圧Pair）の補正に用いられる。
【００５７】
（７）大気温補正項算出部について
大気温補正項算出部７４では、排気行程中の行程容積V、残留ガス質量（排気上死点でのすきま容積Vcに基づいて算出）m、残留ガス（既燃ガス）のガス定数R、および大気温度Tairの実測値を理想気体の状態方程式に代入することで、筒内圧力Pthを算出する。当該筒内圧力Pthと、筒内圧推定モデル６８で算出される筒内圧力Pとの偏差を算出する。そして、その偏差が大きい場合には、上記偏差に基づいて補正係数を算出する。この補正係数は、上述した吸気圧力推定モデル６６において、吸気圧力Pmapの補正に用いられる。
【００５８】
［フリクションモデルの学習］
内燃機関１０において、高精度な空燃比制御と良好なドライバビリティの確保を実現するためには、燃焼状態を正確に推定する手段が必要とされる。上述したエンジンモデル６０によれば、燃焼状態を表す指標として、筒内圧力Pの最大圧力値Pmaxを算出することができる。燃焼状態を正確に推定するうえで、フリクションがもたらす影響は大きなものである。
【００５９】
そのようなフリクションの影響を明確にするためには、環境変化や経年変化に起因して生ずるフリクションのずれが考慮されているのが望ましい。そこで、本実施形態のシステムでは、経年変化等に依らずに正確なフリクション情報を得られるようにすべく、ECU５０に以下の図６に示すルーチンを実行させることとした。
【００６０】
図６は、上記の目的を実現するために、本実施の形態１においてECU５０が実行するルーチンのフローチャートである。尚、本ルーチンは、エンジンモデル６０が高い演算精度を発揮する状態で処理が行われるように、フリクションが適合されている状態で実行されるものとする。ここでいう、フリクションが適合されている状態とは、エンジンモデル６０が算出する推定クランク角度θおよび推定エンジン回転数Neと、実クランク角度θおよび実エンジン回転数Neとの偏差が十分に小さい状態をいう。
【００６１】
図６に示すルーチンでは、先ず、実エンジン回転数Neおよび実クランク角度θが、30°CA毎に検出される（ステップ１００）。次いで、内燃機関１０の現在の運転状態が定常状態にあるか否かが判別される（ステップ１０２）。
【００６２】
上記ステップ１０２において、内燃機関１０が定常状態にあると判定された場合には、エンジンモデル６０を用いて、燃焼圧力Pが算出される（ステップ１０４）。より具体的には、本ステップ１０４では、先ず、上記ステップ１００において検出された実エンジン回転数Neおよび実クランク角度θが初期値としてエンジンモデル６０に入力される。次いで、エンジンモデル６０により演算されるエンジン回転数Neおよびクランク角度θの推定値を上記（４ｅ）式に代入することで、入力トルクTRQが算出される。また、フリクショントルクTRQfはフリクションモデル６４より得ることができ、負荷トルクTRQLはECU５０に記憶されている。従って、上記（５）式を用いて、エンジン発生トルクTRQeを算出することができる。その結果、上記（６ｃ）式に従って、燃焼圧力Pを算出することができる。より詳細には、ここでは、燃焼圧力Pは、クランク角度θに対する圧力履歴として算出される。
【００６３】
次に、気筒毎に、所定のサイクル数に対する最大圧力値Pmaxの平均値が算出される（ステップ１０６）。より具体的には、本ステップ１０６では、上記ステップ１０４において算出された燃焼圧力Pの履歴から、クランク角度θの情報に基づいて、各気筒の最大圧力値Pmaxがそれぞれ特定されたうえで、気筒毎に、所定のサイクル数に対する最大圧力値Pmaxの平均値が算出される。
【００６４】
次に、上記ステップ１０６において算出された最大圧力値Pmaxの平均値が所定の範囲内にあるか否かが、気筒毎に判別される（ステップ１０８）。その結果、上記平均値が所定の範囲内にあると判定された場合には、燃焼状態が悪くなっている気筒はないものと判断され、今回のルーチンの処理が速やかに終了される。
【００６５】
一方、上記ステップ１０８において、上記平均値が所定の範囲から外れている気筒があると判定された場合には、該当する気筒において燃焼状態が悪くなっていると判断され、その該当気筒に対して点火時期の補正が実行される（ステップ１１０）。
【００６６】
次に、全気筒において、上記平均値が所定の範囲を外れているか否かが判別される（ステップ１１２）。その結果、全気筒において、上記平均値が所定の範囲を外れていると判定された場合には、フリクションのずれが発生しているものと判断される。そこで、本ルーチンでは、次いで、エンジンモデル６０を用いて、エンジン回転数Neに応じて、フリクションの学習が実行される（ステップ１１４）。
【００６７】
図７は、本ステップ１１４におけるフリクション学習の手法を説明するための図である。本ステップ１１４では、上記ステップ１００において算出される実エンジン回転数Neとエンジンモデル６０により算出される推定エンジン回転数Neとの偏差に所定のフィードバックゲインを乗じた値として、PIDコントローラ７６によって算出されるPID補正量を、フリクションモデル６４が備えるフリクションマップ（図４参照）のマップ値に反映させるようにしている。
【００６８】
図７は、そのようなフリクションマップの補正の仕方を表している。尚、図７において、破線で示す波形は本ステップ１０４の学習がなされる前のマップ値に、実線で示す波形は当該学習のなされた後のマップ値に、それぞれ対応しており、また、図７中の波形上の丸印および三角印は、それぞれ、学習前後の各マップ値に対応している。
【００６９】
図７に示すように、上記のPID補正量は、ノイズ的な挙動を除去すべく、各マップ点に対する所定の領域を考慮して、当該領域の中で算出された補正量の平均値や時間的な積分値として算出されたものである。このようなPID補正量が各マップ値（丸印の値）に反映されることで、フリクションの値が新たなマップ値（三角印の値）に学習更新される。
【００７０】
以上説明した図６に示すルーチンによれば、エンジンモデル６０を用いて、内燃機関１０の燃焼状態を正確に取得しているので、フリクションモデル６４にずれが生じているか否かを的確に判断することができる。そして、そのような的確な判断に基づいて、フリクションモデル６４のずれを学習し、フリクションモデル６４の出力を更新することができる。このため、経年変化や環境変化に起因して生ずるフリクションモデル６４のずれを補正することができる。その結果として、フリクションモデル６４の正確な出力値が反映されたエンジンモデル６０に基づいて、より正確なエンジン制御を行うことが可能となる。
【００７１】
ところで、従来から、筒内圧力を検出する筒内圧センサの出力値に基づいて燃焼状態を推定する手法が知られている。この筒内圧センサを用いる手法によれば、内燃機関１０の図示トルクを得ることはできるが、当該手法からは、出力軸であるクランク軸１６の軸トルクを得ることはできない。フリクションモデル６４が出力するフリクション情報が正確であるか否かは、図示トルクではなく軸トルクに依らなければ的確に判断することはできない。
【００７２】
上記の筒内圧センサを用いる手法に対し、本実施形態のエンジンモデル６０を用いる手法によれば、上記図６に示すルーチンで説明したように、実機のエンジン回転数Neおよびクランク角度θから燃焼圧力Pが逆算されることになる。このような手法によれば、フリクションの影響が反映されたクランク軸１６の入力トルクTRQ、すなわち、軸トルクに基づいて、燃焼圧力Pを算出することができる。
【００７３】
つまり、本実施形態の手法によれば、内燃機関１０の出力軸であるクランク軸１６の周りで運動方程式を立てているので、クランク軸１６に作用するフリクショントルクTRQfを正確に見積もったうえで、そのフリクションのばらつきが検証される。このため、上記の筒内圧センサを用いる手法に比して、フリクションの影響を明確にすることができる。このように、本実施形態の手法によれば、高価なセンサを必要とせず、より正確なエンジン制御を行うことが可能となる。
【００７４】
尚、上述した実施の形態１においては、ECU５０が、上記ステップ１０４の処理を実行することにより前記第１の発明における「燃焼状態指標算出手段」が、上記ステップ１０６の処理を実行することにより前記第１の発明における「燃焼状態指標平均値算出手段」が、上記ステップ１１２の処理を実行することにより前記第１の発明における「判定手段」が、それぞれ実現されている。また、フリクションモデル６４が前記第１の発明における「フリクション推定手段」に相当している。
また、上述した実施の形態１においては、クランク軸周りの運動方程式演算部６２が前記第３の発明における「演算部」に、PID補正量が前記第３の発明における「フィードバック補正量」に、PIDコントローラ７６が前記第３の発明における「フィードバック手段」に、それぞれ相当している。
【００７５】
実施の形態２．
次に、図８を参照して、本発明の実施の形態２について説明する。
本実施形態のシステムは、図１に示すハードウェア構成および図２に示すエンジンモデル６０の構成を用いて、ECU５０に図６に示すルーチンに代えて後述する図８に示すルーチンを実行させることにより実現することができるものである。
【００７６】
上述した実施の形態１では、燃焼状態を表す指標として、最大圧力値Pmaxを使用している。これに対し、本実施形態では、燃焼状態を表す指標として、上述したWeibe関数のパラメータ（上記（９ａ）式参照）を使用する点に特徴を有している。
【００７７】
［実施の形態２における具体的処理］
図８は、各気筒の燃焼状態に基づきフリクション学習を行うために、本実施の形態２においてECU５０が実行するルーチンのフローチャートである。尚、図８において、実施の形態１における図６に示すステップと同一のステップについては、同一の符号を付してその説明を省略または簡略する。また、本ルーチンについても、フリクションが適合されている状態で実行されるものとする。
【００７８】
図８に示すルーチンでは、エンジンモデル６０を用いて燃焼圧力Pが逆算された後（ステップ１０４）、次いで、燃焼波形算出部（燃焼モデル）７０を利用して、気筒毎にWeibe関数のパラメータの平均値が算出される（ステップ２００）。
【００７９】
より具体的には、本ステップ２００では、上記（９ａ）式および（１０）式を用いて燃焼波形算出部７０によって算出される燃焼圧力Pを、上記ステップ１０４の処理により逆算された燃焼圧力Pと一致させるためのWeibe関数のパラメータ（上述した形状係数m、燃焼効率kなど）が、最小二乗法を用いたWeibe関数の同定によって算出される。本ステップ２００では、そのようなWeibe関数のパラメータの同定が、気筒毎に所定のサイクル数に渡って実行される。そして、所定のサイクル数に対して、気筒毎に同定されたパラメータの平均値が算出される。
【００８０】
次に、上記ステップ２００において算出されたWiebe関数のパラメータの平均値が所定の範囲内にあるか否かが、気筒毎に判別される（ステップ２０２）。その結果、上記平均値が所定の範囲内にあると判定された場合には、燃焼状態が悪くなっている気筒はないものと判断され、今回のルーチンの処理が速やかに終了される。
【００８１】
一方、上記ステップ２０２において、上記平均値が所定の範囲から外れている気筒があると判定された場合には、該当する気筒において燃焼状態が悪くなっていると判断され、その該当する気筒に対して点火時期の補正が実行される（ステップ１１０）。
【００８２】
次に、全気筒において、上記平均値が所定の範囲を外れているか否かが判別される（ステップ２０４）。その結果、全気筒において、上記平均値が所定の範囲を外れていると判定された場合には、フリクションのずれが発生しているものと判断される。そこで、本ルーチンでは、次いで、エンジンモデル６０を用いて、エンジン回転数Neに応じて、フリクションの学習が実行される（ステップ１１４）。
【００８３】
図８に示すルーチンを参照して説明したように、燃焼状態を表す指標としてWeibe関数のパラメータを利用する手法によっても、エンジンモデル６０を用いた軸トルク情報に基づいて内燃機関１０の燃焼状態を正確に取得しているので、フリクションモデル６４にずれが生じているか否かを的確に判断することができる。
【００８４】
尚、上述した実施の形態２においては、燃焼波形算出部７０が前記第４の発明における「燃焼圧力推定手段」に相当している。
【００８５】
実施の形態３．
次に、図９を参照して、本発明の実施の形態３について説明する。
本実施形態のシステムは、図１に示すハードウェア構成および図２に示すエンジンモデル６０の構成を用いて、ECU５０に図６に示すルーチンに代えて後述する図９に示すルーチンを実行させることにより実現することができるものである。
【００８６】
上述した実施の形態１では、燃焼状態を表す指標として、最大圧力値Pmaxを使用している。これに対し、本実施形態では、燃焼状態を表す指標として、PIDコントローラ７６によって算出されるPID補正量を使用する点に特徴を有している。
【００８７】
［実施の形態３における具体的処理］
図９は、各気筒の燃焼状態に基づきフリクション学習を行うために、本実施の形態３においてECU５０が実行するルーチンのフローチャートである。尚、図９において、実施の形態１における図６に示すステップと同一のステップについては、同一の符号を付してその説明を省略または簡略する。また、本ルーチンについても、フリクションが適合されている状態で実行されるものとする。
【００８８】
図９に示すルーチンでは、内燃機関１０の現在の運転状態が定常状態にあると判定された場合には（ステップ１０２）、上記ステップ１００において算出された実エンジン回転数Neと、エンジンモデル６０により算出される推定エンジン回転数Neとの偏差に所定の係数を乗じた値として、PID補正量がPIDコントローラ７６によって算出される（ステップ３００）。尚、ここで算出されるPID補正量は、上記のエンジン回転数Neの偏差を筒内圧力Pにフィードバックするための補正量である。
【００８９】
次に、上記ステップ３００において算出されたPID補正量が、気筒毎に所定のサイクル間で積分され、気筒毎にPID補正量の平均値が算出される（ステップ３０２）。
【００９０】
次に、上記ステップ３０２において算出されたPID補正量の平均値が所定の範囲内にあるか否かが、気筒毎に判別される（ステップ３０４）。その結果、上記平均値が所定の範囲内にあると判定された場合には、燃焼状態が悪くなっている気筒はないものと判断され、今回のルーチンの処理が速やかに終了される。
【００９１】
一方、上記ステップ３０４において、上記平均値が所定の範囲から外れている気筒があると判定された場合には、該当する気筒において燃焼状態が悪くなっていると判断され、その該当気筒に対して点火時期の補正が実行される（ステップ１１０）。
【００９２】
次に、全気筒において、上記平均値が所定の範囲を外れているか否かが判別される（ステップ３０６）。その結果、全気筒において、上記平均値が所定の範囲を外れていると判定された場合には、フリクションのずれが発生しているものと判断される。そこで、本ルーチンでは、次いで、エンジンモデル６０を用いて、エンジン回転数Neに応じて、フリクションの学習が実行される（ステップ１１４）。
【００９３】
図９に示すルーチンを参照して説明したように、燃焼状態を表す指標として、クランク軸周りの運動方程式演算部６２の出力に基づくPID補正量を利用する手法によっても、エンジンモデル６０を用いた軸トルク情報に基づいて、内燃機関１０の燃焼状態を正確に取得しているので、フリクションモデル６４にずれが生じているか否かを的確に判断することができる。
【図面の簡単な説明】
【００９４】
【図１】本発明の実施の形態１の内燃機関の停止位置制御装置が適用された内燃機関の構成を説明するための図である。
【図２】図１に示すECUが備えるエンジンモデルの構成を示すブロック図である。
【図３】クランク軸周りの各要素に付す記号を示す図である。
【図４】図２に示すフリクションモデルがフリクショントルクTRQfを取得するために備えているフリクションマップの一例を示している。
【図５】筒内圧力Pの履歴取得の変形例の手法を説明するための図である。
【図６】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【図７】図６に示すルーチンのステップ１１４におけるフリクション学習の手法を説明するための図である。
【図８】本発明の実施の形態２において実行されるルーチンのフローチャートである。
【図９】本発明の実施の形態３において実行されるルーチンのフローチャートである。
【符号の説明】
【００９５】
１０  内燃機関
１２  ピストン
１４  コンロッド
１６  クランク軸
２４  スロットルバルブ
２６  スロットルポジションセンサ
４０  クランク角センサ
４２  カム角センサ
５０  ECU(Electronic Control Unit)
５２  空燃比センサ
５４  水温センサ
６０  エンジンモデル
６２  クランク軸周りの運動方程式演算部
６４  フリクションモデル
６６  吸気圧力推定モデル
６８  筒内圧推定モデル
７０  燃焼波形算出部
７２  大気圧補正項算出部
７４  大気温補正項算出部
７６  PIDコントローラ
dQ/dθ  熱発生率
dθ/dt  クランク角回転速度
【図１】
0004665788.tif 000014
【図２】
0004665788.tif 000015
【図３】
0004665788.tif 000016
【図４】
0004665788.tif 000017
【図５】
0004665788.tif 000018
【図６】
0004665788.tif 000019
【図７】
0004665788.tif 000020
【図８】
0004665788.tif 000021
【図９】
0004665788.tif 000022