発電機の周波数と電気的出力の関係

# はじめに
発電機の周波数は一般的に重負荷になるほど低下する傾向にある。また、その場合の周波数と電気的出力は負の相関があり、線形関係が近似的に成立する場合が多い。この考えをもとにして速度調定率が定義されている。そこで今回は、そのような線形関係が負荷の変動後の安定状態で成立するかどうかを発電機の動揺方程式という微分方程式を数値解析的に解くことで検証する。ただし、以下のグラフがシミュレーション結果である。

![発電機の周波数と電気的出力の関係.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/270d16a3-fc68-a32a-5039-7493a6907490.png)

![発電機の周波数と電気的出力の関係_等高線図.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/cd2ca504-1549-f702-409e-95816e91e87e.png)

# 導入
## 動揺方程式


回転体の回転エネルギーは慣性モーメント（定数）$J$を用いて以下のようにあらわすことができる。

```math
W=\frac{1}{2}mv^2=\frac{1}{2}m(r\omega)^2=\frac{1}{2}(mr^2)\omega^2=\frac{1}{2}J\omega^2
```
ただし、$J=mr^2$として、$J$は慣性モーメントと定義する。

したがって、エネルギーの変化、つまり仕事率$P$は以下のようにあらわすことができる。

```math
P=\frac{dW}{d\omega}=J\omega \frac{d\omega}{dt}
```

これは、単位時間あたりの仕事である仕事率の定義からもいえる。


```math
P=Fv=(ma)v=(mr\frac{d\omega}{dt})(r\omega)=J\omega \frac{d\omega}{dt}
```

ところで、仕事率$P$は機械的入力$P_m$と電気的出力$P_e$を用いて以下の様に表すことができる。ただし今回は、$P_m$は定数とする。


```math
J\omega \frac{d\omega}{dt}=P_m-P_e
```

これが簡易的な発電機の動揺方程式である。

ゆえに、

```math
\frac{d\omega}{dt}=\frac{P_m-P_e}{J\omega }
```

という式を差分法を用いた解析で用いるものとする。

## アルゴリズム

今回のプログラムによる解析では、$P_e$が与えられたとき、シミュレーション時間$t$を十分大きな値にしてそのときの$\omega$を差分法により算出する。

そこで、以下の様なプログラムを作成した。

```python gene_p_f.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

P_m=1.0
J=10

delta_t=0.01
num=100

t_ary=[]
omega_ary=[]
f_ary=[]

P_e_ary=np.linspace(0,2.0,num)

for i in range(num):
  t=0
  f=50
  omega=2*math.pi*f
  while t<=10:
    P_e=P_e_ary[i]
    d_omega= ((P_m-P_e)/J)*1/omega
    omega = omega + d_omega*delta_t
    t=t+delta_t
    f=omega/(2*math.pi)
  f_ary.append(f)

plt.title("発電機の周波数と電気的出力の関係")
plt.xlabel("電気的出力")
plt.ylabel("周波数")
plt.plot(P_e_ary,f_ary)
plt.savefig("発電機の周波数と電気的出力の関係.png")
plt.show()
```

これを実行すると以下のようなグラフが出力される。

![発電機の周波数と電気的出力の関係.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/270d16a3-fc68-a32a-5039-7493a6907490.png)

予想どおり、$P_m,\omega$には傾きが負となる線形関係が近似的に成立するようである。しかし、これだと、過渡的な応答から安定状態に移るまでの間に、どのように変化していったかが分からない。そこで、以下の様なプログラムを作成した。

```python gene_p_f_contourf.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

P_m=1.0
J=10

delta_t=0.01
num=100

t_ary=[]
omega_ary=[]
f_ary=[]
P_e_ary=np.linspace(0,2.0,num)
P_e_ary2=[]
for i in range(num):
  t=0
  f=50
  omega=2*math.pi*f
  while t<=10:
    P_e=P_e_ary[i]
    d_omega= ((P_m-P_e)/J)*1/omega
    omega = omega + d_omega*delta_t
    t=t+delta_t
    f=omega/(2*math.pi)
    t_ary.append(t)
    omega_ary.append(omega)
    f_ary.append(f)
    P_e_ary2.append(P_e)
    # plt.scatter(P_e,omega/(2*math.pi),c=t)
  # f_ary.append(omega/(2*math.pi))

plt.title("発電機の周波数と電気的出力の関係")
plt.xlabel("電気的出力")
plt.ylabel("周波数")
plt.scatter(P_e_ary2,f_ary,c=t_ary)
plt.colorbar(label="経過時間")
plt.savefig("発電機の周波数と電気的出力の関係_等高線図.png")
plt.show()
```

これを実行すると以下のようになる。

![発電機の周波数と電気的出力の関係_等高線図.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/cd2ca504-1549-f702-409e-95816e91e87e.png)

このように、周波数は初期状態から過渡状態、そして安定状態にいたるまでに変化しているということが分かる。

# まとめ
今回は、速度調定率の根拠となる、発電機の周波数と送電電力が傾き負の線形関係が近似的に成立することについて疑問に思ったので調査を行った。するとシミュレーション結果では確かに、傾き負の線形関係が近似的に成立することが分かった。このような原理は電験の施設管理の科目の『系統の周波数変動($\Delta P=-K \Delta f$)』と深く関わってくるので、苦手な人は復習することを強く勧めたい。