4ストロークエンジンの物理モデル

4ストロークのガソリンエンジンの物理モデリングをしてみました．
シミュレーションはMATLAB/Simulink，可視化ツールはProcessingです．

参考元はこのあたり．
双方ともに高知工科大学の磯村先生の文献です．

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

並列計算機によるエンジンのシミュレータ開発

---

以前にも倒立振子とかパーティクルフィルタでのグリッドマップマッチングとか2リンク平面マニピュレータとかのモデリングをやってきたのですが，そのどれとも色が違うなあ，ということで．

例によって，ラグランジュ法で運動方程式を立てます．
このあたりは参考文献の上の方が詳しいかな．

---

幾何学的関係

定義から，
$$\sin\phi = \frac{r}{L}\sin\theta$$
$$\cos^2\phi = 1-\sin^2\theta$$
代入して，
$$\cos^2\phi = 1-\left(\frac{r}{L}\right)^2\sin^2\theta$$
$$\cos\phi = \sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta}$$

---

最初の式を$t$で微分して，
$$\dot{\phi}\cos\phi = \frac{r}{L}\dot{\theta}\cos{\theta}$$

$$\dot{\phi} = \frac{r}{L}\frac{\cos{\theta}}{\cos\phi}\dot{\theta}$$

$$\dot{\phi} = \frac{r}{L}\frac{\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\dot{\theta}$$

$\cos\phi$に関する式を代入すると，以下の関係が得られます．

$$\dot{\phi}^2 = \left(\frac{r}{L}\right)^2\frac{\cos^2{\theta}}{1-\left(\frac{r}{L}\right)^2\sin^2\theta }\dot{\theta}^2$$

---

定義より，

$$x = (L+r)-(L\cos\phi+r\cos\theta)$$

$t$で微分すると，

$$\dot{x} = L\dot{\phi}\sin\phi+r\dot{\theta}\sin\theta$$

$\dot{\phi}$，$\sin\phi$を代入して，

$$\dot{x} = L \frac{r}{L}\frac{\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\frac{r}{L}\dot{\theta}\sin\theta+r\dot{\theta}\sin\theta$$

$$\dot{x} = r\dot{\theta}\sin\theta\left[1+\frac{\frac{r}{L}\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\right]$$

$$\dot{x}^2 = r^2\dot{\theta}^2\sin^2\theta\left[1+\frac{\frac{r}{L}\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\right]^2$$

---

運動エネルギ

クランク軸の運動エネルギ
$$T_{cr} = \frac{1}{2}I_{cr} \dot{\theta}^2$$
ただし，$I_{cr}$はクランク軸の慣性モーメントで，クランク軸は並進運動をしないため，回転運動の運動エネルギのみとなります．

---

ピストンの運動エネルギ
$$T_{p} = \frac{1}{2}M_{p} \dot{x}^2$$
ただし，$M_{p}$はピストンの質量で，ピストンは回転運動をしないため，並進運動のみの運動エネルギとなります．

---

コネクションロッドの運動エネルギ
$$T_{cn} = \frac{1}{2}\left(I_{cn}\dot{\phi}^2+M_{cn}\dot{x}^2\right)$$
ただし，$I_{cn}$，$M_{cn}$はそれぞれコネクションロッドの慣性モーメントおよび質量です．

---

全運動エネルギ
$$T = T_{cr}+T_{p}+T_{cn}$$

$$T = \frac{1}{2}I_{cr} \dot{\theta}^2+\frac{1}{2}M_{p} \dot{x}^2+\frac{1}{2}\left(I_{cn}\dot{\phi}^2+M_{cn}\dot{x}^2\right)$$

$$T = \frac{1}{2}\left[I_{cr} \dot{\theta}^2+I_{cn}\dot{\phi}^2+\left(M_{p} +M_{cn}\right)\dot{x}^2 \right]$$

$\dot{\phi}^2$と$\dot{x}^2$を代入すると，

$$T = \frac{1}{2}\left[I_{cr} \dot{\theta}^2+I_{cn}\left(\frac{r}{L}\right)^2\frac{\cos^2{\theta}}{1-\left(\frac{r}{L}\right)^2\sin^2\theta }\dot{\theta}^2+\left(M_{p} +M_{cn}\right)r^2\dot{\theta}^2\sin^2\theta\left(1+\frac{\frac{r}{L}\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\right)^2 \right]$$

$\dot{\theta}$でまとめると，以下の関係が得られます．

$$T = \frac{1}{2}\left[I_{cr}+I_{cn}\left(\frac{r}{L}\right)^2\frac{\cos^2{\theta}}{1-\left(\frac{r}{L}\right)^2\sin^2\theta }+\left(M_{p} +M_{cn}\right)r^2\sin^2\theta\left(1+\frac{\frac{r}{L}\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\right)^2 \right]\dot{\theta}^2$$

---

ここで，

$$f(\theta) = I_{cr}+I_{cn}\left(\frac{r}{L}\right)^2\frac{\cos^2{\theta}}{1-\left(\frac{r}{L}\right)^2\sin^2\theta }+\left(M_{p} +M_{cn}\right)r^2\sin^2\theta\left(1+\frac{\frac{r}{L}\cos{\theta}}{\sqrt{1-\left(\frac{r}{L}\right)^2\sin^2\theta} }\right)^2$$

とおくと，全運動エネルギ$T$は以下に書き換えることができます．

$$T = \frac{1}{2}f(\theta) \dot{\theta}^2$$

---

ラグランジュ法

ラグランジアン

$$L = T-U$$

ただし，位置エネルギ$U$は無視できるため実際は，

$$L=T$$

となります．
ラグランジュ法では系全体の運動方程式は，

$$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}}-\frac{\partial L}{\partial \theta} = Q$$

から求めることができます．

---

左辺第1項
$$\frac{\partial L}{\partial \dot{\theta}} = f(\theta) \dot{\theta}$$

$$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}} = \frac{d f(\theta)}{dt}\dot{\theta}+f(\theta) \ddot{\theta}$$

---

左辺第2項
$$\frac{\partial L}{\partial \theta} = \frac{1}{2}\frac{\partial f(\theta)}{\partial \theta} \dot{\theta}$$

ここで注意したいのですが，MATLAB/SimulinkではTime-domainでシミュレーションを実行するため，$\theta$での微分項は少し使いづらいです．
（前サンプリング時刻の$\theta$と現在のサンプリング時刻の$\theta$の差で除算する必要があり，クランク角速度ゼロの状態ではゼロ割になってしまう）
このため，$\dot{\theta}$を$d\theta/dt$として，少し変形しています．

$$\frac{\partial L}{\partial \theta} = \frac{1}{2}\frac{d f(\theta)}{dt}$$

これで左辺は求まりました．

---

右辺
右辺は，爆発により生じるトルク（外力）です．
クランク軸に作用するトルク$Q$は，幾何学的関係から以下となります．

$$Q = PA_{p}r\sin\theta$$

ただし，$P$はシリンダ内圧力，$A_{p}$はピストンの断面積です．

文献ではもう少し複雑な式になっていますが，上の式で合っているような気がします．

これで運動方程式は求まりました．

---

燃焼系

燃焼系はまともにモデリングするのは骨が折れるので，簡単化します．
具体的には，吸気，圧縮，膨張，排気にそれぞれの行程で以下の仮定を置いています．

• 吸気
  シリンダ内圧力は大気圧（1気圧）で一定．
  管内摩擦も無視し，吸気行程が始まった瞬間からエンジンの慣性力に任せるままの運動となる．

• 圧縮
  シリンダ内圧力は断熱圧縮．
  気体の状態方程式により圧力が決まる．
  圧縮開始時の圧力は大気圧とする．（つまり吸気，排気によってシリンダ内部の圧力が初期状態に戻る）

• 膨張
  シリンダ内圧力は断熱膨張．
  圧縮行程の逆の動きをする．
  膨張開始時の圧力は定数とする．
  圧縮から膨張への切り替わり時（爆発時）には等積加熱の動作を行う．（微小時間で圧力が最大値まで上昇する）

• 排気
  シリンダ内圧力は大気圧で一定．
  管内摩擦も無視する．
  エンジンの慣性力のみで動作する．

これらの仮定を基に，行程をモデル化しています．
クランク軸が2回転でサイクルのため，クランク軸角度$\pi$(rad)ごとに区分しています．

---

吸気 $(0\leq\theta<\pi)$

$$P = P_{Air}$$

ただし，$P_{Air}$は大気圧です．

---

圧縮 $(\pi\leq\theta<2\pi)$

$$P = P_{Air}\left(\frac{V_{bdc}}{A_{p}x+V_{clr}}\right)^{\kappa}$$

$$x = (L+r)-(L\cos\phi+r\cos\theta)$$

$$x = (L+r)-\left(\sqrt{L^2-r^2\sin^2\theta}+r\cos\theta\right)$$

ただし，$V_{bdc}$は下死点でのシリンダ容積，$V_{clr}$は上死点での隙間容積，$\kappa$は比熱比です．

---

膨張 $(2\pi\leq\theta<3\pi)$

$$P = P_{Exp}\left(\frac{V_{clr}}{A_{p}x+V_{clr}}\right)^{\kappa}$$
$$x = (L+r)-\left(\sqrt{L^2-r^2\sin^2\theta}+r\cos\theta\right)$$

ただし，$P_{Exp}$は爆発時の圧力です．

---

排気 $(3\pi\leq\theta<4\pi)$

$$P = P_{Air}$$

---

これで全ての式が導出できました．
これらをSimulinkで表現すると，次のような感じです．

全体

---

$f(\theta)$の中身

---

シリンダ内圧力

---

圧縮行程

---

膨張行程

---

$\theta$から$x$への変換

このモデルでシミュレーションした結果が，冒頭の動画です．
リンクによって拘束される系のダイナミクスは解いたことがなかったのでなかなか面白かったように思います．

通常だと下死点から始まるため，初期位置を爆発直後の位置にずらしています．
スタータとかを考慮する場合には燃焼部分を少しいじる感じでしょうか．