使用 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}g{t}^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{g{x}^2}{2v_0^2}\cdot (1+{\tan }^2\theta )$$
若已知目標物所在的位置

$(x,y)$，則初速度仰角

$\theta$ 可以由

$\tan \theta$ 為變數的一元二次方程式解出。

$$g{x}^2\cdot {\tan }^2\theta -2v_0^{2}x\cdot \tan \theta +2v_0^{2}y+g{x}^2=0$$

$$\tan \theta =\frac{v_0^{2}x\pm \sqrt{v_0^4{x}^2-g{x}^2(2v_0^{2}y+g{x}^2)}}{g{x}^2}$$

用影片中給定的條件

$v_0=828\mathrm{m}/\mathrm{s}$、

$x=20\mathrm{km}$、

$y=500\mathrm{m}$、

$g=9.8\mathrm{m}/{\mathrm{s}}^2$ 代入上式可得

$$\tan \theta \approx 6.82\mathrm{或}0.17\theta \approx {81.66}^{\circ }\mathrm{或}{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 sv{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 sv{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}/\mathrm{s}$、質量

$m=43.5\mathrm{kg}$，目標物水平距離

$x=20\mathrm{km}$、高度

$y=500\mathrm{m}$，無因次的空氣阻力係數

$C=0.2$、空氣密度

$\rho =1.3\mathrm{kg}/{\mathrm{m}}^3$、砲彈截面積

$s=0.01815{\mathrm{m}}^2$，重力加速度

$g=9.8\mathrm{m}/{\mathrm{s}}^2$，時間間隔

$dt=0.01\mathrm{s}$，計算結果

$\theta \approx {24.6}^{\circ }$，與理論計算的結果

$\theta \approx {81.66}^{\circ }\mathrm{或}{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，這篇文章裡的程式就是由它改寫而成的。以後在講到斜向拋射時，可以拿這篇文章的主題作為例子，說明理想狀態與實際情境有多大的差異。

---

tags:Physics、VPython、Google 試算表