使用 Google 試算表及 VPython 計算砲彈軌跡

作者：王一哲
日期：2022/3/2

前言

前幾天在 YouTube 上看到李永樂老師的影片「如何精準命中目標？戰爭到底帶給我們什麽？」，影片中提到理想狀態與實際情境下的斜向拋射運動差異，最後有提到計算道彈的數值方法。影片中李永樂老師是用 Excel 計算的，而我在高三多元選修中則是教學生用 VPython 計算，雖然工具不同，但是原理相同。

影片來源：https://youtu.be/nalzSo4f0iQ

理想狀態

若砲彈的初速度量值為 $v_0$、仰角為 $\theta$，重力加速度為 $g$，只考慮重力的作用，則砲彈的水平位移 $x$ 與時間 $t$ 的關係為 $$x = v_0 \cos \theta \cdot t ~\Rightarrow~ t = \frac{x}{v_0 \cos \theta}$$ 鉛直位移 $y$ 與時間 $t$ 的關係為 $$y = v_0 \sin \theta \cdot t - \frac{1}{2}gt^2$$
將 $t = \frac{x}{v_0 \cos \theta}$ 代入 $y$ 可得軌跡方程式 $$y = v_0 \sin \theta \cdot \frac{x}{v_0 \cos \theta} - \frac{1}{2}g \cdot \left( \frac{x}{v_0 \cos \theta} \right)^2 ＝ \tan \theta \cdot x - \frac{g}{2v_0^2 \cos^2 \theta} \cdot x^2$$ 由於 $$\frac{1}{\cos^2 \theta} = \sec^2 \theta = 1 + \tan^2 \theta$$ 可以將軌跡方程式改寫成 $$y = x \cdot \tan \theta - \frac{gx^2}{2v_0^2} \cdot (1 + \tan^2 \theta)$$ 若已知目標物所在的位置 $(x, y)$，則初速度仰角 $\theta$ 可以由 $\tan \theta$ 為變數的一元二次方程式解出。 $$gx^2 \cdot \tan^2 \theta - 2v_0^2 x \cdot \tan \theta + 2v_0^2 y + gx^2 = 0$$ $$\tan \theta = \frac{v_0^2 x \pm \sqrt{v_0^4 x^2 - gx^2 (2v_0^2 y + gx^2)}}{gx^2}$$
用影片中給定的條件 $v_0 = 828 ~\mathrm{m/s}$、$x = 20 ~\mathrm{km}$、$y = 500 ~\mathrm{m}$、$g = 9.8 ~\mathrm{m/s^2}$ 代入上式可得 $$\tan \theta \approx 6.82 ~或~ 0.17 ~~~~~ \theta \approx 81.66^{\circ} ~或~ 9.77^{\circ}$$

實際情境

若砲彈速度 $v$ 往右上方，速度與水平方向夾角為 $\theta$，則水平分量 $v_x = v \cos \theta$，鉛直分量 $v_y = v \sin \theta$。砲彈受到空氣阻力 $f$ 方向與速度相反，其量值為 $$f = \frac{1}{2} C \rho s v^2$$ 其中 $C$ 為無因次的阻力係數，$\rho$ 為空氣密度、$s$ 為物體截面積。因此砲彈受到向左下方的力量 $F$，其水平方向分量 $$F_x = f \cos \theta = \frac{1}{2} C \rho s v^2 \cos \theta = \frac{1}{2} C \rho s v v_x$$ 鉛直方向分量 $$F_y = mg + f \sin \theta = mg + \frac{1}{2} C \rho s v^2 \sin \theta = mg + \frac{1}{2} C \rho s v v_y$$ 由牛頓第二運動定律可以計算砲彈的加速度 $$a_x = \frac{F_x}{m} ~~~~~~~~~~ a_y = \frac{F_y}{m}$$
由於力量、加速度、速度會互相影響，很難將軌跡方程式寫出來，但是我們可以用數值方法計算軌跡。最簡單的作法是用時刻 $t$ 的速度計算力量，代入一小段時間 $dt$，計算砲彈在 $dt$ 內的速度變化以及時刻 $t + dt$ 時的速度；再用更新後的速度乘以 $dt$，計算砲彈在 $dt$ 內的位移以及時刻 $t + dt$ 時的位置；接著用更新後的速度計算時刻 $t + dt$ 時的力量；不斷重複以上的過程就可以算出軌跡。

使用數值方法計算軌跡的流程圖

我仿照影片中的作法用 Google 試算表做了以下的檔案，連結在此，有興趣的同學可以將檔案另存複本拿回去修改看看。如果採用影片中的數值，砲彈初速度 $v_0 = 828 ~\mathrm{m/s}$、質量 $m = 43.5 ~\mathrm{kg}$，目標物水平距離 $x = 20 ~\mathrm{km}$、高度 $y = 500 ~\mathrm{m}$，無因次的空氣阻力係數 $C = 0.2$、空氣密度 $\rho = 1.3 ~\mathrm{kg/m^3}$、砲彈截面積 $s = 0.01815 ~\mathrm{m^2}$，重力加速度 $g = 9.8 ~\mathrm{m/s^2}$，時間間隔 $dt = 0.01 ~\mathrm{s}$，計算結果 $\theta \approx 24.6^{\circ}$，與理論計算的結果 $\theta \approx 81.66^{\circ} ~或~ 9.77^{\circ}$ 差距相當大。

使用 Google 試算表及數值方法計算砲彈軌跡

VPython

接下來改用 VPython 計算砲彈，將初速度仰角 $\theta$ 由 $1^{\circ}$ 開始代入程式中計算砲彈與目標物最接近的距離，每次增加 $0.1^{\circ}$ ，直到 $29.9^{\circ}$ 為止。計算的結果為，當 $\theta = 24.6^{\circ}$ 時，砲彈與目標物的距離最近，量值約為 3.37 m。

模擬程式畫面截圖，為了使畫面較為清楚，每隔 2° 畫一條軌跡。

砲彈與目標物最接近距離與初速度仰角關係圖

"""
 VPython教學: 計算砲彈軌跡及砲彈與目標物最近的距離
 日期: 2022/3/2
 作者: 王一哲
 原始資料: https://youtu.be/nalzSo4f0iQ
"""
from vpython import *
import matplotlib.pyplot as plt

"""
 1. 參數設定, 設定變數及初始值
"""
v0 = 828             # 砲彈初速度量值
m = 43.5             # 砲彈質量
C = 0.2              # 無因次的空氣阻力係數
rho = 1.3            # 空氣密度
s = 0.01815          # 砲彈截面積
cof = 0.5*C*rho*s    # 空氣阻力係數
g = 9.8              # 重力加速度 9.8 m/s^2
L = 21000            # 地板長度
t = 0                # 時間
dt = 0.01            # 時間間隔
target = vec(20000, 500, 0)   # 目標物位置
ratio = 20           # 顯示箭頭的長度比例
deg_list, t_list, dis_list = [], [], []

"""
 2. 畫面設定
"""
scene = canvas(title="Cannon Trajectory", width=800, height=400, x=0, y=0,
               center=vec(0.5*L, 0.1*L, 0), background=color.black, range=0.3*L)
floor = box(pos=vec(0.5*L, -0.5, 0), size=vec(L, 1, 0.1*L), color=color.blue)
ball = sphere(pos=vec(0, 0, 0), radius=100, color=color.red)

# 開啟檔案 range.csv, 屬性為寫入, 先寫入欄位的標題
with open("CannonTrajectoryData.csv", "w", encoding="UTF-8") as file:
    file.write("theta(degree), t(s), distance(m)\n")

"""
 3. 物體運動部分
"""
def motion(degree):
    t = 0
    theta = radians(degree)
    ball = sphere(pos=vec(0, 0, 0), radius=1, color=color.red,
                  make_trail=True, v=vec(cos(theta), sin(theta), 0)*v0)
    arrow_mg = arrow(pos=ball.pos, shaftwidth=20, axis=vec(0, 0, 0), color=color.green)
    arrow_f = arrow(pos=ball.pos, shaftwidth=20, axis=vec(0, 0, 0), color=color.cyan)
    dis1 = mag(target - ball.pos)
    while True:
        rate(1000)
        f = -cof*ball.v.mag2*ball.v.norm()
        ball.a = f/m + vec(0, -g, 0)
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
        dis2 = mag(target - ball.pos)
        arrow_mg.pos = ball.pos
        arrow_mg.axis = vec(0, -m*g, 0) * ratio
        arrow_f.pos = ball.pos
        arrow_f.axis = f*ratio
        t += dt
        if dis2 < dis1: dis1 = dis2
        else:
            arrow_mg.visible = False
            arrow_f.visible = False
            return t - dt, dis1

for degree in arange(1, 30, 0.1):
    t, distance = motion(degree)
    deg_list.append(degree)
    t_list.append(t)
    dis_list.append(distance)
    print(degree, t, distance)
    with open("CannonTrajectoryData.csv", "a", encoding="UTF-8") as file:
        file.write(str(degree) + "," + str(t) + "," + str(distance) + "\n")

# 印出砲彈與目標物最近的距離及對應的初速度仰角
print(deg_list[dis_list.index(min(dis_list))], min(dis_list))

# 用 matplotlib.pyplot 繪圖
plt.figure(figsize=(6, 4.5), dpi=100)                        # 設定圖片尺寸
plt.xlabel(r'$\theta ~\mathrm{({}^{\circ})}$', fontsize=14)  # 設定坐標軸標籤
plt.ylabel(r'$distance ~\mathrm{(m)}$', fontsize=14)
plt.xticks(fontsize=12)                                      # 設定坐標軸數字格式
plt.yticks(fontsize=12)
plt.grid(color='grey', linestyle='--', linewidth=1)          # 設定格線顏色、種類、寬度
plt.plot(deg_list, dis_list, linestyle='-', linewidth=4, color='blue')   # 繪圖並設定資料點格式
plt.savefig('CanonTrajectoryPlot.svg')                       # 儲存圖片
plt.savefig('CanonTrajectoryPlot.png')
plt.show()                                                   # 顯示圖片

程式碼運作的流程大致上與 Google 試算表的計算流程相同。為了節省程式運作所需時間，我在自訂函式 motion 裡將前一個時刻 t 的砲彈與目標物距離儲存在變數 dis1，將現在的砲彈與目標物距離儲存在變數 dis2，若 dis2 < dis1 代表砲彈正在接近目標物，繼續運作 while 迴圈；若 條件不成立代表砲彈已經遠離目標物，執行 else 當中的程式碼，執行到第66行的 return 時回傳 前一個時刻 t - dt 及最接近的距離 dis1 並且結束 while 迴圈。如果只想要得到初速度仰角與最接近的距離，可以刪除第51行的 rate(1000)，但是模擬程式的畫面會變得很亂。

結語

我在之前的 VPython 講義中有寫過類似的程式：程式 6-3.斜向抛射, 使用for 迴圈改變仰角 theta, 空氣阻力係數 b，這篇文章裡的程式就是由它改寫而成的。以後在講到斜向拋射時，可以拿這篇文章的主題作為例子，說明理想狀態與實際情境有多大的差異。

---

HackMD 版本連結：https://hackmd.io/@yizhewang/SyJPMVie9