平成 12 年度卒業論文

4 サイクルガソリンエンジンの動特性について

平成 13 年 2 月 28 日

高知工科大学

工学部

知能機械システム工学科

下 司
佐々木

真 也
兵 太

指導教員

教授 磯村

修一

目次
―――――――――――――――――――――――――――――――――― 1

緒言
第１章

エンジンシミュレータの開発の背景、必要性、目的 ―――――― 1

第２章

エンジン・シミュレータ開発上の問題点 ――――――――――― 2

第３章

シミュレータ・デザイン ―――――――――――――――――― 4

３．１

シミュレーション・プログラム開発手順 ――――――――― 4

３．２

数式モデル

３．３

クランク軸の運動方程式 ―――――――――――――――― 6―――――――――――――――――――――― 4

３．３．１

エネルギの式 ――――――――――――――――――― 6

３．３．２

Euler−Lagrange による運動方程式 ――――――――

３．３．３

トルクについて ―――――――――――――――――― 10

３．３．４

クランク軸系の運動方程式 ――――――――――――― 119――――――――――――――――――――― 11

３．４

ピストン変位

３．５

コンロッドの座標 ――――――――――――――――――― 12

３．６

クランクの座標 ―――――――――――――――――――― 14

第４章

開発環境

――――――――――――――――――――――――― 15―――――――――――――――――――――― 15

４．１

開発ツール

４．２

プログラム開発、実行手順 ――――――――――――――― 15
シミュレーション結果 ――――――――――――――――――― 16

第５章
結言―――――――――――――――――――――――――――――――――― 17

謝辞―――――――――――――――――――――――――――――――――― 17―――――――――――――――――――――――――――――――― 18

参考文献
付録１

シミュレーション結果の連続画像

付録２

連続合成画像

付録３

運動方程式について

付録４

シミュレーション結果のファイルについて

緒言
エンジン・シミュレータにはさまざまな目的・用途があるが、本研究では、加速性と制御の問題を解決するために、エンジンの運動系のダイナミックス（動特性）に注目した。
エンジンのダイナミックシステムは（１） 燃焼系統（２） エンジンの本体の運動系（３） 過給機の動特性からなる一種のフィードバック（特にポジティブフィードバック）システムである点に特徴がある。このシステムの動特性研究では、これらの機器を設計、製作あるいは購入して、解析と試験を実施したいところであるが本学の設備、予算の面から困難であった。
そこで動特性の同定を目的としたエンジンシ動特性ミュレータ開発に取り組もうとしたのである。
必然的にエンジンの燃焼効率や摩擦損失など静的性能、静特性に関する詳細な検討は範囲外となり、エンジンの過渡的な性能に関する部分に開発努力を集中した。
対象としてはガスタービン、またはターボ過給機付ディーゼルエンジンが適していると思われるが、動特性モデル作成用入力データが 4 サイクル‑ガソリンエンジン用しか入手できなかったので、テーマを「4 サイクル‑ガソリンエンジンの動特性」とした。

第１章

エンジンシミュレータ開発の背景、必要性、目的

精密かつ高速のシミュレータ開発は動く機械の製造業にとって最重要課題の一つと言える。コンピュータ内に作られたリアルタイム・シミュレータによって工場や現場でのテス
トを省くことができ、コスト削減や工期の短縮といったメリットが生まれる。また、実機による試験中の事故や試験要員の確保といった問題も解消されるため、今後ますますこのようなシミュレーション技術は必要になると思われる。
一方、産業界での IT 化とダウンサイジング化が急速に進んでいる現状を考え、従来メイ
ンフレーム上の FORTRAN で行われていた科学技術計算を、PC 上のインターネット言語 Javaで実行することは、コンピュータ・グラフィクス（ＣＧ）が脚光を浴びている現状では価値があると思われる。Java は情報家電や携帯端末の開発用言語として将来性があり、数値計算のライブラリも充実してきているので、将来機械製造業にも浸透していくものと考えられる。
本研究の対象であるエンジンについては、レシプロエンジンやガスタービンエンジンは輸送手段である自動車、船舶、航空機、列車の駆動用として採用され、高性能化の要求、
環境基準などを次々にクリアし 20 世紀にはほぼ完成の域に達し、もはや研究テーマは残さ

1

れていないのではないかと思われていた。ところが、最近の電力自由化でビルや家庭向けなどの分散型発電機の需要が世界的に高まる中、米国のガスタービンメーカーと日本のエンジンメーカーが互いの得意分野を補完し合う業界再編も加速される可能性がでてきた。
このような小型エンジンが輸送機関ではない発電機に採用された場合の解決すべき技術
問題点は何か？
効率の追求や環境対策のための技術は輸送機関の場合と特に変わっていないであろう。
今後クローズアップされるのは最新鋭火力発電設備と同様、加速性能と制御であると思われる。そのためには不可欠のエンジン動特性解明の研究に製造業各社は取り組み始めるであろう。
加速性能の優れたマイクロガスタービンが発電用として有力と言われる中、低価格燃料、高効率で優位に立つレシプロエンジンの最大のネックはその加速性にある。現在のレシプ
ロエンジンは効率向上のために必ず排気ガスターボチャージャー（過給機）を装備しているが、負荷上昇時限り過給機の追従遅れのため、出力上昇が追いつかず、発電機の回転数が低下してなかなか回復しないという問題が予測される。そのためビルや家庭の電源周波数が 60Hz（又は 50Hz）から降下し電力系統におけるさまざまな悪影響が予想されるのである。
エンジンのシミュレータと言っても、その範囲は広く、次のような用途がある。
1. 「部品同士が当たって壊れないか？」「スパナが入るか？」などの組み立て時の初歩
的なミスチェックのための CAD/CAM に属するもの
2. 「動くか？」「安全に回転するか？」「爆発しないか？」など異常な圧力や応力が発生しないかを確かめるもの
3. 「訓練・保守に実物の代わりとして使えるか？」4. 「性能が出るか？」「効率は計画通りか？」などを確かめるもの
5. 「材料強度は大丈夫か？」「局部的に応力が高くないか？」を確かめるもの
6. 「振動による騒音はないか？」さらには「破壊しないか？」を確かめるもの
7. 「負荷変動の要求に対しすばやく応答するか？」加速性能を確かめるもの
8. 「疲労破壊の心配はないか？」を確かめるもの
9. 「摩耗で壊れないか？」を確かめるもの

第２章

エンジン・シミュレータ開発上の問題点

このようにエンジン・シミュレータの目的・用途は多岐にわたっており、全ての要求に応じられるシミュレータを開発することは、容易ではない。
先ず、エンジンシミュレーションは Statics と Dynamics に分けられる。Dynamics はさらに Snapshot(Batch processing)と Real‑Time processing(Animation)に分けられる。

2

Statics
Simulation

Snap shot
Dynamics
Real-Time

Statics を計算するシミュレータはエンジンの性能を求めることから Heat balance
program と呼ばれている。このプログラムに対する入力はエンジンの負荷、Fuel、Air など
であり、出力はトルク、排気ガス温度（これは過給機への入力エネルギでもある）である。
このプログラムの目的はエンジンの効率をいかに上げるかを計算することなので、性能にかかわる因子例えば熱力学的性質、燃焼室の形状寸法などのデータが理論あるいは実験式の形で与えられねばならない。
一方、Dynamics の場合、全ての状態変数は時間の関数として扱われる。エンジンの中にはゆっくり変化するものから、非常に速く変化するものまで含まれている。これは固有値の値が広範囲にわたっていることを意味する。例えば、物理学で電磁力、弱い力、強い力、重力の 4 つの力があるようなものである。
ゆっくり変化するものには負荷変動に対するエンジン出力の追従性があり、速く変化す
るものにはクランクシャフトの捩り振動や排気弁（動弁とも言う）の上下振動現象がある。
これらはスティッフな問題とも言われ高い固有振動数を持つ。そのため微分方程式の積分には時間スケールを引き伸ばさないと莫大な計算時間がかかってしまう。
この動弁系の問題はバネを分布系として、つまり無限の固有値とモードを持つ偏微分方
程式として扱わねばならず、5 年と 1 億円をかけて計算、実験を繰り返した例もある。
物理学の 4 つの力を１つの式で計算することが困難であるように、エンジンのあらゆる動的な現象を一挙にシミュレートすることは不可能に近いといえる。
そこで用途別に、つまり目的向き（Object Oriented）にシミュレータを作ることになる。
エンジンシステムのように要素が複雑に構成されているシステムでは、１本のパイプの中の流れの計算のような分けにはいかない。部分的、機能毎に分割して数式モデルを作り、1つ 1 つ積み重ねるしかないのである。部分では意味が無いと言われればそれまでであるが、大学学部の卒業研究では陣営、設備、予算上の制約がありとてもエンジン全体まではできないのである。極論すれば吸気過程の計算すらも完全にはできないであろう。
そのためには前述した目的に応じたシミュレータを作り、使い分けることになる。
エンジンのシミュレータを作成するにあたっては、全システムの微分方程式を一気に書き下ろして計算してしまわけにはいかないので、次のようにシステムを分ける。
1. ピストン、クランク、コンロッドなどの運動部分
2. 燃焼系を吸気、圧縮、爆発、膨張、排気の各過程に分けそれぞれを解析していく部分
3. 過給機（ガスタービン、コンプレッサー、回転部分を含む）
4. 制御システム（演算部、E/H 変換部）

3

これらを一歩一歩進めて最後に統合する。
境界条件の扱いについては、回転部分のシミュレータにピストン圧力が必要となるが、シリンダガス圧を求めて加えるだけでは不十分である。これはピストンの動きがガス圧にも影響するからである。このことは出力インピーダンスが 0 ではないということであり、全体システムのシミュレータで考慮されなければならない。

第３章

シミュレータ・デザイン

３．１

シミュレーション・プログラム開発手順

運動系シミュレーション・プログラムの設計は以下のように行った。
・ 各々の要素を剛体とする正確な数式モデルを作成する。
・ クランク軸系の主要部品、すなわちピストン・クランク・コンロッド及び静止したシリ
ンダーとクランクケースについて、CG 用のインプットデータ１）２）を作成し、個々に CG用に作画する。
・ ピストン・クランク・コンロッドの各要素について集中モデルのスケルトンに関する運動方程式を立て、時間領域で解くプログラムを作る。
・ 直線運動のピストン、複雑な直線・回転運動のコンロッドの分布質量を回転運動のクランク軸に集中させる。
・ スケルトンモデルの運動に、作画したピストン・クランク・コンロッドの CG を貼り付
け、クランクケースに収めて運動系のシミュレータを完成させる。

３．２

数式モデル３）４）

数式モデルの作成に先立って以下のように仮定した。
・ エンジン駆動部分をピストン・クランク・コンロッドの各要素に分け剛体として扱う。
・ ピストンは直線運動、コンロッドは直線運動と同時に回転運動も行うがすべてクランク
軸回りの回転運動に換算する。
・ ピストン‑クランク軸系の運動エネルギを求める。
・ ピストンにかかるガス圧からクランク軸にかかるエンジン駆動トルクを本来求めなけ
ればならないが、今回はプログラムで仮の爆発トルクを与える。
・ 起動時はセルモータのトルクで駆動する。
・ 負荷トルクは自動車の走行速度の二乗に比例するとし、クランク軸の角速度の二乗に比例するトルクとして与える。
・ 駆動系への入力は仮想の爆発によるトルクとセルモーターの駆動トルクとする。
・ 駆動系への出力はピストンとコンロッドの変位、クランク軸の回転角とする。
ピストン、コンロッド、クランク軸の関係の断面図を図 1 に示す。

4

ω

O
ｒ
θ＋φ

θ
D

C

図中の記号説明O 点：クランク軸中心
A 点：ピストンの上死点
B 点：ピストンの下死点
P 点：コンロッドのピストンピン

B

L

C 点：コンロッドのクランクピン
Ｌ ：コンロッド長さ
ｒ ：クランク半径φθ ：クランク角度
φ ：コンロッド傾き角

S＝２ｒ

P
X
A

図 1. ピストン-コンロッド-クランク軸の関係図

5

３．３

クランク軸の運動方程式

３．３．１

エネルギの式

クランク軸の回転運動エネルギー Tk

1
Tk＝ Ikθ& 2
2

[1-1]

ピストンの運動エネルギー Tp

1
Tp＝ mpX& 2
2

[1-2]

コンロッドの回転及び直線運動のエネルギー Tc

T c＝

(

1 &2
Icφ + mcX& 2
2

)

[1-3]

ここで、

Ik ：クランク軸の慣性モーメント Ic ：コンロッドの慣性モーメント
mp ：ピストンの質量
mc ：コンロッドの質量である。
全運動のエネルギーは

(

1
1
1
T＝ Ikθ& 2 + mpX& 2 + Icφ& 2 + mcX& 2
2
2
2

)

1
1
1
＝ Ikθ& 2 + (mp + mc )X& 2 + Icφ& 2
2
2
2

[1-4]

図 1 から X 、 θ 、 φ の関係が次のように求められる。

L sin φ = r sin θ[1-5]

X = (L + r ) − (L cos φ + r cosθ )[1-6]

式[1-5]から

sin φ =

r
sin θ
L

[1-7]

cos 2 φ = 1 − sin 2 φ
2r
cos φ = 1 −   sin 2 θ
L
2

2

r
cosφ = 1 −   sin 2 θ
L

[1-8]

6

式[1-7]を時間 t で微分して

r
cos φφ& = cosθθ&Lr cosθ &φ& =θ
L cos φ
=

r
L

cosθ
2

r
1 −   sin 2 θLrφ& 2 =  
L

2

θ&

cos 2 θ2r
1 −   sin 2 θ
L

[1-9]

θ& 2

[1-10]

が得られる。
また式[1-6]を時間 t で微分して次式を得る。

X& = L sin φφ& + r sin θθ&rr
X& = L sinθ
L
L

cosθ
2

r
1 −   sin2 θ Lθ& + r sinθθ&



r


cosθ

&
L
X& = r sin θ 1 +
θ
2


r
1 −   sin 2 θ 

L






2
2
2
X& = r sin θ 1 +


[1-11]

2


r

cos θ
 &2
L
 θ
2

r 
1 −   sin 2 θ 
L


7

[1-12]

結局、全運動のエネルギー T = T (θ ,θ ′) は
2



r


cosθ
1 &2 1 &2
2
2 L
T＝ Ikθ + (mp + mc )r sin θ 1 +
θ
2
2
2
r


 
1 −   sin 2 θ 

L


+

1 r
Ic 
2  L

2

cos 2 θ2r
1 −   sin 2 θ
L

θ& 2

2






r


2
cosθ
2
 2
r
θ
1
cos


 
L
T＝  Ik + (mp + mc )r 2 sin 2 θ 1 ++ Ic 
θ&

2
2
2
L
r


r
2
1 −   sin 2 θ 
1
−
sin
θ






L
 L





[1-13]となる。ここで
2



r


2
cosθ
cos 2 θ

r
2
2 
L
+
f (θ ) = Ik + (mp + mc )r sin θ 1 +
I
c



2
2
L
r


r
2
1 −   sin 2 θ
1 −   sin θ 

L
L


[1-14]とおくと、運動エネルギー[1-13]はT=

1f (θ )θ& 2
2[1-15]となる。

8

３．３．２Euler-Lagrange による運動方程式５）Lagrangian は

L = T −Uと表されるが、重力を無視すれば U = 0 より

L =T
となる。
Euler-Lagrange の方程式d ∂L ∂L−=Ｑ
dt ∂θ& ∂θ

[1-16]を求める。
式[1-15]から∂L= f (θ )θ&∂θ&d ∂L ∂f (θ ) dθ &=
θ + f (θ )θ&&dt ∂θ&∂θ dt
d ∂L ∂f (θ ) & 2=
θ + f (θ )θ&&dt ∂θ&∂θ∂L 1 ∂f (θ ) & 2=
θ∂θ 2 ∂θd ∂L ∂L ∂f (θ ) & 2
1 ∂f (θ ) & 2−=
θ + f (θ )θ&& −θ&
dt ∂θ ∂θ∂θ2 ∂θとなり、式[1-16]は1 ∂f (θ ) & 2
f (θ )θ&& +θ = Ｑ [1-17]
2 ∂θ
となる。
ここに、 f (θ ) は式[1-14]で与えられる。
また、クランク軸にかかる全トルク Ｑ は

Ｑ = Ｑe + Ｑc − ＱL [1-18]とする。
ここで Ｑe ：エンジントルク、 Ｑc ：セルモータートルク、 ＱL ：負荷トルクである。

9

３．３．３

トルクについてエンジントルクＱe についてだが、ピストンに加えられる力は PA で与えられる。ここに、P は圧力、 A はピストン断面積である。これをクランク軸にかかるトルクに換算する。
コンロッドに働く力は

PA cos φ
クランクピンの接線方向にかかる力は

PA cos φ sin (θ + φ )したがって、クランク軸にかかるトルク Ｑe は

Ｑe = PAr cos φ sin (θ + φ )[1-19]となる。式[1-7]と式[1-8]から

sin φ =

r
sin θ
L
2

r
cosφ = 1 −   sin 2 θ
L
であるから

Ｑe = PAr cos φ (sin θ cos φ + sin φ cosθ )  r  2 2  r
Ｑe = PAr 1 −   sin θ  + cosθ  sin θ L  L 

[1-20]

負荷トルクは車のスピード、すなわちクランク軸の回転数に比例するので

ＱL = ζθ& θ&[1-21]と表される。ここに ζ は比例定数である。

10

３．３．４クランク軸系の運動方程式

クランク軸系の運動方程式は
2






r


2
cosθ
2


r
θ
cos
1 ∂f (θ ) & 2



2
2 
L
+ Ic θ + ζθ& θ& Ik + (mp + mc )r sin θ 1 +
θ&& +

2
2
L
2
∂
θ


r



r
1 −   sin 2 θ 
1 −   sin 2 θ 



L
L





  r  2 2  r
= PAr 1 −   sin θ  + cosθ  sin θ + Ｑc
 L  L 

[1-22]となる。

３．４

ピストン変位

図 1 において、ピストンの上死点からのピストン変位を X とすると

X = OA − OP
= OA − (DP + DO )である。ストローク（行程）を S 、クランク半径を r = S

2、コンロッドの長さを L 、シリ

ンダ中心線からのクランク角度を θ 、コンロッドの傾き角度を φ とし、

λ = Lrとおくと、

OA = L + r
DP = L cos φ
DO = r cosθであるから、

X = (L + r ) − (L cos φ + r cosθ )X = r (1 − cosθ ) + λr (1 − cos φ )である。式[2-1]の φ をクランク角 θ で表すにはr sin θ = L sin φ[2-2]の関係を用いてcosφ = 1 − sin 2 φ

[2-3]

11

[2-1]

sin φ =
=

r
sin θ
L

1
sin θ
λ

[2-4]
2

r
cosφ = 1 −   sin 2 θ
L

[2-5]

式[2-5]を式[2-1]に代入




1
X = r 1 − cosθ + λ 1 − 1 − 2 sin 2 θ 
λ




[2-6]

が得られる。
式[2-6]から、クランク軸系の運動方程式[1-22]で得られた θ を用いてピストン変位 X を計算することができる。

３．５

コンロッドの座標

図 1 においてクランク軸中心 O は座標（0,0）に固定されている。クランクピン中心座標を（ xc 、 yc ）とすると、

xc = r sin θ
yc = r cosθ

[3-1]
[3-2]

となる。
ピストンピン中心の座標（ x p 、 y p ）は

xp = 0[3-3]

y p = yc + L cos φ[3-4]である。式[2-4]から

sin φ =

r
sin θ
L

[3-5]
2

r
cosφ = 1 −   sin 2 θ
L

[3-6]

の関係があり、 sin θ から sin φ と cos φ を求めることができる。

12

コンロッドの運動は P 点を中心に回転すると同時に、 P 点が上下に運動するので、コンロッドの形状をあらわすデータ（座標）は回転中心 P を原点として求める。

P の座標を（ x p 、 y p ）、コンロッドの任意の相対座標を（ ∆x1 、 ∆y1 ）とすると、 P 点が上下し、C 点が P 点の回りに φ だけ回転するコンロッドの座標（ ∆x1 、∆y1 ）は（ ∆x2 、∆y2 ）に移されるとする。
ｙ（ ∆x2 、 ∆y2 ）

Ｌ

Ｌ

φ

（ ∆x1 、 ∆y1 ）αｘ
図 2．コンロッド内の座標変換
図 2 から∆x1 = L cosα∆y1 = L sin α∆x2 = L cos(φ + α )∆y2 = L sin (φ + α )

[3-7]
[3-8]
[3-9]
[3-10]となる。
ここに α はコンロッドの図形のノードをベクトルで表したときの方向を示している。
加法定理

cos(φ + α ) = cosφ cosα − sin φ sin α
sin (φ + α ) = sin φ cosα − sin α cos φ

を用いて
座標変換は

∆x2  cos φ∆y  =  sin φ 2 

− sin φ   L cosα 
cos φ   L sin α 13cos φ= sin φ− sin φ   ∆x1 cos φ  ∆y1 [3-11]となる。
コンロッド上の原点（ x p 、 y p ）はピストンの上下運動座標と等しいので、すなわちxp = 0y p = yc + L cos φx2 = x p + ∆x2 = ∆x2y2 = y p + ∆y2 = yc + L cos φ + ∆y2であるから、結局、求めるコンロッドの座標（ x2 、 y2 ）は0
 x2   cos φ y  = r cosθ + L cos φ  +  sin φ  2 

− sin φ  ∆x1 cos φ  ∆y1 [3-12]

が得られる。

３．６

クランクの座標

クランク軸回りのクランク各点の座標は点 O を中心とする回転のみなので、 θ の回転によって点（ ∆x1 、 ∆y1 ）は点（ ∆x2 、 ∆y2 ）に移される

∆x2  cosθ∆y  =  sin θ 2 

− sin θ  ∆x1 cosθ  ∆y1 [3-13]

14

第４章

開発環境

４．１

開発ツール

ノート型パソコン（TOSHIBA DynaBook Satellite2060）、OS（Windows98 Second
Edition）を使用、これにインターネット経由で Sun Microsystems（http://java.sun.com/）
から Java

2 SDK Standard Edition1.2.2(JDK1.2.2)をダウンロードし、インストールしたものが主な開発ツールである。

４．２

プログラム開発、実行手順

以下に簡単な開発、実行手順を示す。
・ テキストエディタで Java のプログラムソース６）７）を作成する。
・ JDK1.2.2 のコンパイラでソースをコンパイルする。
・ コンパイル正常終了後に Class ファイルが作成される。
・ アプレットを実行するための HTML ファイルを作成する。
appletviewer、NetscapeNavigator、InternetExplorer などで HTML ファイルをブラウズ。

15

第５章

シミュレーション結果

ピストン-クランク-コンロッドを質点とみなした運動系、すなわち集中モデルのスケルトンに関する微分方程式を時間領域で解いたリアルタイム・シミュレーション結果の 1 部を
図 3 に示す。また、このスケルトンモデルの運動に、作画したピストン-クランク-コンロッ
ドのＣＧ画像を貼り付けた運動系のリアルタイム・シミュレーション結果の 1 部を図 4 に示す。また、付録Ａにスケルトンモデルとリアルタイムシミュレータ、付録Ｂにリアルタ
イムシミュレータを連続的に画像として取り込んだものを付属する。

図 3．スケルトンモデル図

図 4．エンジン動的状態図

16結言
ピストン、クランク、コンロッドを含むエンジン運動系動特性シミュレーションを行ってコンピュータ・グラフィックス（CG）表示した。トルクは本来、ピストン面にかかる動的ガス圧を燃焼計算から求めるべきであるが、今回はセルモーターと仮想の爆発によるトルクを与えた。
今後の課題としては、燃焼系のシミュレーション、CG の 3 次元化、多気筒エンジンのシミュレータ等を作成し、最終的には車体を含めた自動車シミュレータの作成を行いたい。
来年以降の卒業研究に期待したい。

謝辞
磯村研究室で開発中の TURBO C 言語によるエンジン動特性プログラムを基に、本研究ではその一環として、Java 言語を用いてインターネット上で動く運動系リアルタイムシミュレーションプログラムを開発した。
また、本研究を進めるにあたり、指導していただいた磯村先生、及び協力を惜しまなかった都築弘孝君にも感謝の意を表したい。

17

参考文献
1）

卒業論文「４サイクルエンジン燃焼系動特性シミュレーション」
今垣真介，谷口貴史
高松工業高等専門学校

2）

機械工学科

1997,2,28

卒業論文「４サイクルガソリンエンジンの動特性シミュレーション」
朝日一平，中林宏之
高松工業高等専門学校

3）

機械工学科

1997,2,28

S,Isomura et al, The Development of Automation for Mitubishi UE Diesel Engine"Studies of Transient Characteristics of Two-Stroke Cycle Full Turbocharged
Diesel Engine"
CIMAC Permanent Committee. Bruxelles, 1968, p5-p84）S,Isomura et al, Dynamic Characteristics of Two-Cycle Full Turbocharged Diesel
Engines,
Technical Review, Mitubishi Heavy Industries, Ltd. Vol.3 No.5, 1966. 9

5）

M.A.LAMPSA; Maximizing Distance of the Golf Drive: An Optimal Control
Study, Transactions of the ASME Dec. 1975, p362-p367

6） 「JAVA による図形処理入門」
山本芳人
工学図書株式会社
7）「やさしい Java プログラミング 改定版」
萩本順三
AI 出版

18

付録１

シミュレーション結果の連続画像

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

付録２

連続合成画像34

付録３

運動方程式についてNewton の運動方程式では慣性力以外の、例えば Coriolis の力などすべての力を式に正しく組み込むのは難しい。運動及び位置エネルギーから Euler-Lagrange の方程式を求める方が間違えない。
例えば惑星の運動について考えると次のようになる。

x = r cosθy = r sin θ

Kinetic energy： T =

ｒ
θ

ｍ

M

(

1
m x& 2 + y& 2
2

Potential energy： U = −)m：惑星を質点とみなした質量γ ⋅m⋅Mr

M：太陽の質量γ：万有引力定数∞ U = − ∞ γ ⋅ m ⋅ M dξ =  γ ⋅ m ⋅ M   ξ ∫r ξ 2r Lagrangian

L： L = T − U

L = T −U=γ ⋅m⋅M1
m x& 2 + y& 2 +
2
r

(

)

x& = r& cosθ − rθ& sin θ
y& = r& sin θ − rθ& cosθ

遠心力： mrθ& 2

x& 2 + y& 2 = r& 2 + r 2θ& 2
T=
L=

(

1
m r& 2 + r 2θ& 2
2

)

引力： −

γ ⋅m⋅M1
m r& 2 + r 2θ& 2 +
2
r

(

)

γ ⋅m⋅M
r2
θ

35

ｍ

x

θに関してｒに関して∂L= mr 2θ&
&
∂θ

∂L
= mr&
∂r&

d ∂L
= 2mrr&θ& + mr 2θ&&
&
dt ∂θ

d ∂L
= m&r&
dt ∂r&

∂L
=0
∂θ

∂L
γ ⋅m⋅M= mrθ& 2 −∂rr2d ∂L ∂L−=0
dt ∂θ& ∂θd ∂L ∂L−=0
dt ∂r& ∂r

mr 2θ&& + 2mrr&θ& = 0
ここでは 2mrr&θ& Coriolis の力

γ ⋅m⋅M
m&r& − mrθ& 2 +=0
r2

mr 2θ&& + 2mrr&θ& = 0 は
d 2&mr θ =0dt2と書けるから r θ& = const が得られる( )これは Kepler の法則と呼ばれる。

r = const のときはどうか？
m = const であるから、慣性能率

I = mr 2 = constとなる。
これは円板のようなものが垂直な中心軸の回りに回転しているようなもの。

mr 2θ&& + 2mrr&θ& = 0に代入して

mr 2θ&& = 0Iθ&& = 0
θ& = const 、 θ& = ω とおくと
γ･m･M
m&r& − mrθ& 2 +
=0
r2
より

36

mrω 2 =

ω=

γ･m･M
r2

γ･M
r3

γ･m･M
がつりあった円軌道上の等速運動になる。
r2
r ≠ const. 、例えば r = r (t ) のとき Coriolis の力 2mrr&θ& が働き、楕円軌道になる。

すなわち遠心力 mrθ& = mrω と引力
22

このように簡単な質点運動力学でも奥が深い。レシプロエンジンにおける回転むらの最大の要因は勿論ある回転角の範囲で働くガス圧の変動であるが、慣性能率の回転角に対する変動も影響していることが以上の議論でわかる。
蒸気機関では、車輪にカウンターバランス質量を付加しているし、自動車エンジンにもフライホイールをつけたり、多気筒にして回転むらを平滑化している。

エンジンの慣性能率を I 、フライホイールのを I F 、全体のを I T とするとIT = I F + I ≒ I F

IF ≫ Iとすることにより、滑らかに回る円板に近づく。
ところが、 I F を大きくすると加速性能が悪くなる。発電機では回転むら（これは電源の周波数変動を引き起こす）を少なくし、加速性能を改善するアクティブな制御システムの開発が期待される所以である。
一般の場合

( )

慣性能率が I = I θ ,θ& と θ 、θ& の関数となる場合について考える。ゴルフスイング、フィギュアスケートのスピン、宇宙ステーション、エンジン回転系などがその例である。
(1) Newton の運動方程式
慣性能率×角加速度＝トルク
とすると

( )

I θ ,θ& ･θ&& = Ｑとなる。
(2) Momentum の変化から運動方程式を求める

( )

d &
Iθ = Ｑdt
dI &θ + Iθ&& = Ｑ
dtIθ&& + I&θ& = Ｑ

37

I = mr 2 のとき
I&θ& = 2mrr&θ&となるから I&θ& は Coriolis の力である。
(3) Euler‐Lagrange Equation

T=

( )

1
I θ ,θ& ･θ& 2
2

L = T −U =（簡単化のため U = 0 とする）

( )1
I θ ,θ& ･θ& 2
2( ) ( )

∂L 1 ∂=I θ ,θ& + I θ ,θ& ･θ&∂θ& 2 ∂θ&

( )

( )

( )

( )

( )

( )

( )

( )

d ∂L 1 d ∂∂=I θ ,θ& ･θ& 2 +I θ ,θ& ･θ&θ&& + I θ ,θ& ･θ& + I θ ,θ& ･θ&&&&dt ∂θ 2 dt ∂θ
∂θ&

( )

∂L 1 ∂=I θ ,θ& ･θ& 2∂θ 2 ∂θをd ∂L ∂L−=Ｑ
dt ∂θ& ∂θ
に代入して

( )

( )

( )

1 d ∂
& ･θ& 2 + ∂ I θ ,θ& ･θ&θ&& + I θ ,θ& ･θ& + I θ ,θ& ･θ&& − 1 ∂ I θ ,θ& ･θ& 2 = ＱIθ,θ
2 dt ∂θ&
∂θ&
2 ∂θ
よって

( )

( )

( )

∂
∂
1d ∂



I θ ,θ& −I θ ,θ& θ& 2 + I θ ,θ& ･θ& = Ｑ I θ ,θ& + & I θ ,θ& θ& θ&& + &∂θ
2  dt ∂θ
∂θ




( )
( )
( )
∂
∂
1d ∂

I (θ ,θ& )･θ&θ&& 、 I (θ ,θ& ) −I (θ ,θ& )θ& が加わっていることがわかる。
&&∂θ2 dt ∂θ∂θI θ ,θ& が Coriolis の力である。Newton の慣性項 I θ ,θ& ･θ&& 、Coriolis の力 I θ ,θ& ･θ& 以外に
2

慣性能率がスイングアングル θ の関数となるゴルフスイングの dynamics では、このような理由から、Euler‐Lagrange Equation から出発している。

38

付録４

今回作成したスケルトンモデル、リアルタイムシミュレーションを
CD-R で論文に添付する。
ファイル名は以下の通りである。
・スケルトンモデル：skeleton_model.html
・リアルタイムシミュレーション：Engine.html

39