鉄道のPythonコード


---

#### 【1】列車の速度

$$
v = 2\pi r \cdot \frac{N}{60}
$$

* $v$：列車の速度（m/s）
* $r$：車輪半径（m） 
* $N$：車輪回転数（rpm）

---

#### 【2】運動方程式（加速度運動）

$$
F_{\text{牽引}} - F_{\text{抵抗}} = M \cdot a
$$

* $F_{\text{牽引}}$：モータによる牽引力（N）
* $F_{\text{抵抗}}$：走行抵抗（N）
* $M$：列車質量（kg）≒ 300,000
* $a$：加速度（m/s²）

---

#### 【3】国鉄式抵抗モデル（近似）

$$
F_{\text{抵抗}} = A + Bv + Cv^2
$$

例：

$$
F_{\text{抵抗}} = 1500 + 150v + 900v^2 \quad [\text{N}]
$$

* $A$：転がり抵抗
* $B$：機械損失
* $C$：空気抵抗

---

#### 【4】出力と速度・牽引力の関係式

$$
P = F \cdot v \quad \Rightarrow \quad F = \frac{P}{v}
$$

* $P$：出力（W）
* $F$：牽引力（N）
* $v$：速度（m/s）
  ※VVVF区間では $P = \text{一定} \Rightarrow F \propto 1/v$

---

#### 【5】橋梁上のたわみ（梁力学）

$$
\delta = \frac{F \cdot L^3}{48EI}
$$

* $\delta$：たわみ量（m）
* $F$：荷重（N）
* $L$：橋のスパン長（m）
* $E$：ヤング率（Pa）
* $I$：断面二次モーメント（m⁴）

---

#### 【6】横揺れの固有振動数

$$
f = \frac{1}{2\pi} \sqrt{\frac{k}{m}}
$$

* $f$：固有振動数（Hz）
* $k$：ばね定数（N/m）
* $m$：車体質量（kg）

---

#### 【7】三相交流電源（VVVF制御基礎）

（a）**時間領域での三相電圧**：

$$
\begin{aligned}
v_a(t) &= V_m \sin(\omega t) \\
v_b(t) &= V_m \sin\left(\omega t - \frac{2\pi}{3}\right) \\
v_c(t) &= V_m \sin\left(\omega t + \frac{2\pi}{3}\right)
\end{aligned}
$$

* $V_m$：最大電圧
* $\omega = 2\pi f$：角速度（rad/s）

---

（b）**三相交流の電力**：

$$
P = \sqrt{3} \cdot V_L \cdot I_L \cdot \cos\phi
$$

* $V_L$：線間電圧（V）
* $I_L$：線電流（A）
* $\phi$：力率角

---

（c）**VVVF制御とPWM周波数制御**

* PWM（Pulse Width Modulation）により、出力サイン波を近似
* 周波数 $f$ を変化させて回転数 $n$ を制御

$$
n = \frac{60f}{p}
$$

* $n$：モータ回転数（rpm）
* $f$：周波数（Hz）
* $p$：極数（磁極対数）

---

#### 【8】サイン波の工学的意義

* AC波形はサイン波が基本
* モータでは **回転磁界** を生成
* ベクトル合成・フーリエ級数の基底にもなる

---

```python
# Program Name: train_dynamics_185kei_colab.py
# Creation Date: 20250611
# Overview: Integrated dynamics simulation of 185 series train using physical formulas
# Usage: Paste and run in Google Colab. Calculates speed, traction, resistance, deflection, vibration, and motor RPM.

# --- Install necessary libraries ---
!pip install numpy matplotlib pandas

# --- Import libraries ---
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# --- Unified Parameters for 185 series EMU ---
r = 0.46                          # 車輪半径 [m]
N = np.linspace(400, 700, 100)   # 回転数 [rpm]
M = 300_000                      # 列車質量 [kg]
a = 0.8                          # 加速度 [m/s²]
P = 2_400_000                    # 出力 [W]

A = 1500                         # 転がり抵抗 [N]
B = 150                          # 機械損失 [N·s/m]
C = 900                          # 空気抵抗 [N·s²/m²]

L = 30                           # 橋梁スパン [m]
E = 2.0e11                       # ヤング率 [Pa]
I = 0.05                         # 断面二次モーメント [m^4]

k = 1.5e6                        # 横ばね定数 [N/m]
m_car = 35_000                   # 車体質量 [kg]

Vm = 1000                        # 三相最大電圧 [V]
f = np.linspace(10, 120, 100)    # インバータ周波数 [Hz]
p = 4                            # 極数

# --- Calculations ---
v = 2 * np.pi * r * (N / 60)                # 速度 [m/s]
F_resist = A + B * v + C * v**2             # 抵抗力 [N]
F_traction = P / v                          # 牽引力 [N]
F_load = (M / 10) * 9.81                    # 1両荷重 [N]
delta = (F_load * L**3) / (48 * E * I)      # 橋梁たわみ [m]
f_natural = (1 / (2 * np.pi)) * np.sqrt(k / m_car)  # 横揺れ固有振動数 [Hz]
n_motor = 60 * f / p                        # モータ回転数 [rpm]

# --- 表形式でデータ確認 ---
df = pd.DataFrame({
    "Rotational Speed [rpm]": N,
    "Speed [m/s]": v,
    "Resistance Force [N]": F_resist,
    "Traction Force [N]": F_traction
})
df.head()

# --- Plot: Resistance & Traction Force vs Speed ---
plt.figure(figsize=(10, 6))
plt.plot(v, F_resist, label='Resistance Force [N]')
plt.plot(v, F_traction, label='Traction Force [N]')
plt.xlabel('Speed [m/s]')
plt.ylabel('Force [N]')
plt.title('Force vs Speed (185 series EMU)')
plt.legend()
plt.grid(True)
plt.show()

# --- Plot: Motor RPM vs Frequency (VVVF Control) ---
plt.figure(figsize=(10, 6))
plt.plot(f, n_motor)
plt.xlabel('Inverter Frequency [Hz]')
plt.ylabel('Motor Speed [rpm]')
plt.title('Motor Speed vs Frequency (VVVF Controlled)')
plt.grid(True)
plt.show()

# --- Summary Outputs ---
print(f"Bridge Deflection (1 car load): {delta:.4f} m")
print(f"Natural Frequency of Lateral Vibration: {f_natural:.3f} Hz")
print("Motor RPM for f = [10, 20, 30, 40, 50] Hz:", 60 * np.array([10, 20, 30, 40, 50]) / p)


```