一個地球人在臺北: 馬克士威速率分布圖形

2020年7月19日星期日

馬克士威速率分布圖形

作者：王一哲
日期：2020/7/19

原理

全名為馬克士威 - 波茲曼速率分布 (Maxwell–Boltzmann distribution)，氣體分子的移動速率不會完全相同，氣體分子處於某個速率的機率密度函數為

$$
P(v) = 4\pi \left( \frac{m}{2\pi kT} \right)^{\frac{3}{2}} v^2 e^{-\frac{mv^2}{2kT}}
$$

上式中 $m$ 是原子質量，$k$ 是波茲曼常數 (Boltzmann constant) $1.38 \times 10^{-23} ~\mathrm{J/K}$，$T$ 是溫度，$v$ 是速率。下圖是維基百科上的圖片。不過我們只要知道函數的樣子，應該就能用函數繪圖軟體畫出一樣的圖片。

Wikipedia 原圖，圖片來源 https://en.wikipedia.org/wiki/File:MaxwellBoltzmann-en.svg

如果想要用 VPython 模擬出類似的結果，請參考 Glowscript 網站範例。

Glowscript 網站範例，圖片來源 https://www.glowscript.org/#/user/GlowScriptDemos/folder/Examples/program/HardSphereGas-VPython

使用 Matplotlib 繪圖

import numpy as np
import matplotlib.pyplot as plt

vmin, vmax, num = 0, 2500, 100   # 速度最小值, 最大值, 數量
v = np.linspace(0, 2500, 1000)   # 速度
amu = 1.66053904E-27             # 原子質量單位
mass = {'He': 4*amu, 'Ne': 20*amu, 'Ar': 40*amu, 'Xe': 132*amu}   # 惰性氣體質量
T = 298.15     # 溫度
k = 1.38E-23   # 波茲曼常數

# 產生機率分布資料
He = 4*np.pi*(mass['He']/(2*np.pi*k*T))**(3/2)*v**2*np.exp(-(mass['He']*v**2)/(2*k*T))
Ne = 4*np.pi*(mass['Ne']/(2*np.pi*k*T))**(3/2)*v**2*np.exp(-(mass['Ne']*v**2)/(2*k*T))
Ar = 4*np.pi*(mass['Ar']/(2*np.pi*k*T))**(3/2)*v**2*np.exp(-(mass['Ar']*v**2)/(2*k*T))
Xe = 4*np.pi*(mass['Xe']/(2*np.pi*k*T))**(3/2)*v**2*np.exp(-(mass['Xe']*v**2)/(2*k*T))

plt.figure(figsize=(6, 4.5), dpi=100)                # 設定圖片尺寸
plt.subplots_adjust(left=0.2)                        # 設定圖片左側位置
plt.xlabel('Speed (m/s)', fontsize=14)               # 設定坐標軸標籤
plt.ylabel('Propability Density (s/m)', fontsize=14)
plt.xticks(fontsize = 12)                            # 設定坐標軸數字格式
plt.yticks(fontsize = 12)
plt.grid(color='cyan', linestyle='--', linewidth=1)  # 設定格線顏色、種類、寬度
plt.xlim(vmin, vmax)                                 # 設定x軸繪圖範圍
plt.ylim(0, 0.0045)                                  # 設定y軸繪圖範圍
# 繪圖並設定線條顏色、寬度、圖例
line1, = plt.plot(v, He, color='blue', linewidth=3, label=r'$\mathrm{{}^4 He}$')
line2, = plt.plot(v, Ne, color='green', linewidth=3, label=r'$\mathrm{{}^{20} Ne}$')
line3, = plt.plot(v, Ar, color='red', linewidth=3, label=r'$\mathrm{{}^{40} Ar}$')
line4, = plt.plot(v, Xe, color='orange', linewidth=3, label=r'$\mathrm{{}^{132} Xe}$')
plt.legend(handles = [line1, line2, line3, line4], loc='upper right', fontsize=12)
plt.savefig('MaxwellDistribution.svg')   # 儲存圖片
plt.savefig('MaxwellDistribution.png')
plt.show()                               # 顯示圖片