2018年3月29日 星期四

斜向抛射

作者：王一哲
日期：2018/3/29




將一個小球由地面以初速 v0、仰角 𝜃 抛出，小球受到重力作用向下加速，計算小球的飛行時間及水平射程。分為以下3種不同的狀況：

  1. 只考慮重力的作用，小球撞到地板時停止運動。　（GlowScript 網站動畫連結）
  2. 同時考慮重力及空氣阻力 $f = -bv$，同時畫出考慮空氣阻力與不考慮空氣阻力的小球。　（GlowScript 網站動畫連結）
  3. 使用 for 迴圈，分別計算不同仰角 𝜃、不同的空氣阻力係數 b 對應的飛行時間和水平射程。　（GlowScript 網站動畫連結）

成果如下：


只考慮重力的作用，小球撞到地板時停止運動






同時考慮重力及空氣阻力，b = 0.1




同時考慮重力及空氣阻力，b = 0.1




使用 for 迴圈改變仰角 𝜃 = 15º ~ 90º




使用 for 迴圈改變空氣阻力係數 b = 0 ~ 0.9



程式 6-1：斜向抛射, 球落地停止運作

取得程式碼

"""
 VPython教學: 6-1.斜向抛射, 球落地停止運作
 Ver. 1: 2018/2/21
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1              # 小球半徑
v0 = 30               # 小球初速
theta = radians(30)   # 小球抛射仰角, 用 radians 將單位轉為弧度
L = 100               # 地板長度
g = 9.8               # 重力加速度 9.8 m/s^2
t = 0                 # 時間
dt = 0.001            # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Projection", width=800, height=400, x=0, y=0,
               
center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)
ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.red, make_trail=True,
              v=vec(v0*cos(theta), v0*sin(theta), 0), a=vec(0, -g, 0))

"""
 3. 物體運動部分
"""
while(ball.pos.y - floor.pos.y >= size):
    rate(1000)
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    t += dt

print(t, ball.pos.x + L/2)




參數設定


在此定義的變數有 size、v0、theta、L、g、t、dt，用途都已經寫在該行的註解中。這裡出現了一個新的函式

radians(x)

這個函式可以將 x 的單位由角度轉換為弧度(rad)，radians 包含於 math 函式庫中，使用前需要引入函式庫，但是 vpython 已經包含了 math，不需要額外再引入 math。需要這樣做是因為 python 的三角函數輸入值皆以 rad 為單位。另一個類似的函式為

degrees(x)

這個函式可以將 x 的單位由 rad 轉換為角度。




畫面設定


畫面設定部分與程式 5-1 非常類以，不同之處在於小球的初位置及初速

ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.red, make_trail=True,
              v=vec(v0*cos(theta), v0*sin(theta), 0), a=vec(0, -g, 0))

小球出發時高度為0，初速度 $v_x = v_0 \cos \theta$、$v_y = v_0 \sin \theta$。



物體運動


物體運動部分與程式 5-1 非常類以，不同之處在於小球碰到地板時停止動畫，因此 while 迴圈中設定的條件為

ball.pos.y - floor.pos.y >= size

在 while 迴圈結束之後，印出飛行時間 t 及水平射程 R。可以手動修改 v0 或 theta 的數值，觀察 t 和 R 的變化。



程式 6-2.斜向抛射, 球落地停止運作, 有空氣阻力


取得程式碼

"""
 VPython教學: 6-2.斜向抛射, 球落地停止運作, 有空氣阻力
 日期: 2018/2/21
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1              # 小球半徑
m = 1                 # 小球質量
v0 = 30               # 小球初速
theta = radians(30)   # 小球抛射仰角, 用 radians 將單位轉為弧度
L = 100               # 地板長度
b = 0.1               # 空氣阻力 f=-bv
g = 9.8               # 重力加速度 9.8 m/s^2
t = 0                 # 時間
dt = 0.001            # 時間間隔
t1, t2 = [], []       # 小球落地時間

"""
 2. 畫面設定
"""
scene = canvas(title="Projection with Air Resistance", width=800, height=400,
               x=0, y=0, center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)
# ball: 有空氣阻力; ball2: 沒有空氣阻力
ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.blue, make_trail=True,
              v=vec(v0*cos(theta), v0*sin(theta), 0))
ball2 = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.red, make_trail=True,
               v=vec(v0*cos(theta), v0*sin(theta), 0), a=vec(0, -g, 0))

"""
 3. 物體運動部分, 小球觸地停止運動, 記錄落地時間, 畫出兩者軌跡
    印出落地時間及水平射程, 由於第一次記錄的才是落地時間, 故只印出第一次的值
"""
while(ball.pos.y - floor.pos.y >= size or ball2.pos.y - floor.pos.y >= size):
    rate(1000)
    f = -b*ball.v
    ball.a = vec(0, -g, 0) + f/m
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    ball2.v += ball2.a*dt
    ball2.pos += ball2.v*dt
    if(ball.pos.y - floor.pos.y <= size):
        ball.v = vec(0, 0, 0)
        t1.append(t)
    if(ball2.pos.y - floor.pos.y <= size):
        ball2.v = vec(0, 0, 0)
        t2.append(t)
    t += dt

print(t1[0], ball.pos.x + L/2, t2[0], ball2.pos.x + L/2)



參數設定


程式 6-2 與 6-1 很像，但是增加了小球質量 m 及空氣阻力係數 b。如果對於加上空氣阻力的寫法已經很熟悉，其實可以跳過 6-1 直接寫 6-2。



畫面設定


為了與不考慮空氣阻力的理想狀況對照，需要畫出兩個球，ball 是考慮空氣阻力的狀況、ball2 是理想狀況。由於 ball 的加速度與速度有關，放在 while 迴圈裡再設定即可。另外為了記錄落地時間，需要新增兩個空白的串列 t1、t2。



物體運動


  1. 為了讓兩個小球都碰到地板後才停止動畫，因此 while 迴圈中設定的條件為
    2. ball.pos.y - floor.pos.y >= size or ball2.pos.y - floor.pos.y >= size
  2. 由於 ball 需要考慮空氣阻力，需要加上以下的程式碼，用原來的速度計算空氣阻力，再代入 $F = ma$ 更新加速度。
    f = -b*ball.v
ball.a = vector(0, -g, 0) + f/m

  3. 為了記錄飛行時間，在小球碰到地板之後用 t1.append(t) 及 t2.append(t) 將時間 t 新增到串列的最後面。但其實只有第一個儲存的值才是真正的飛行時間，所以最後只印出串列當中索引值為 0 的元素。這是個偷懶的作法，比較好的寫法應該是只記錄小球剛碰到地板時的時間 t 就好。


程式 6-3.斜向抛射, 使用for 迴圈改變仰角 theta, 空氣阻力係數 b

取得程式碼 
"""
 VPython教學: 6-3.斜向抛射, 使用 for 迴圈改變仰角 theta, 空氣阻力係數 b
 Ver. 1: 2018/3/29
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1                 # 小球半徑
m = 1                    # 小球質量
v0 = 30                  # 小球初速
degree = 30              # 小球抛射仰角, 單位為角度
theta = radians(degree)  # 小球抛射仰角, 用 radians 將單位轉為弧度
L = 100                  # 地板長度
#b = 0.0                  # 空氣阻力 f=-bv
g = 9.8                  # 重力加速度 9.8 m/s^2
t = 0                    # 時間
dt = 0.001               # 時間間隔
n = 10                   # 於 for 迴圈當中代入數值的次數, 改變 theta 時設為16, 改變 b 時設為10

"""
 2. 畫面設定
"""
scene = canvas(title="Projection with for loop", width=800, height=400, x=0, y=0,
               center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)

# 開啟檔案 data.csv, 屬性為寫入, 先寫入欄位的標題
file = open("data.csv", "w", encoding="UTF-8")
file.write("v0(m/s), theta(degree), b, t(s), R(m)\n")

"""
 3. 物體運動部分
"""
for i in range(0, n):
# 改變 theta 使用以下 2 行
#    degree = 15 + 5*i
#    theta = radians(degree)
# 改變 b 使用以下 1 行
    b = 0.1*i
    ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=vec((n-i)/n, 0, i/n),
                  make_trail=True, v=vec(v0*cos(theta), v0*sin(theta), 0))
    while(ball.pos.y - floor.pos.y >= size):
        rate(1000)
        f = -b*ball.v
        ball.a = vec(0, -g, 0) + f/m
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
        t += dt
    print(v0, degree, b, t, ball.pos.x + L/2)
    file.write(str(v0) + "," + str(degree) + "," + str(b) + "," + str(t) + "," + str(ball.pos.x + L/2) + "\n")

file.close() # 關閉檔案

程式設計部分

程式 6-3 是將 5-5、6-1、6-2 合為一體的產物，有兩個使用方法：
  1. 改變仰角 theta，範圍為 15º ~ 90º，每次增加 5º，總共執行 16 次。為了記錄代入的角度數值，需要增加變數 degree，再用 radians 把單位換成 rad、量值存入 theta。
    2. degree = 15 + 5 * i
theta = radians(degree)
  2. 改變空氣阻力係數 b，範圍為 0 ~ 0.9，每次增加 0.1，總共執行 10 次。
    3. b = 0.1 * i
  3. 為了方便使用者從圖中看出執行的順序，我將小球的顏色加上一點變化，一開始小球是紅色的，隨著執行次數增加會慢慢變為藍色。
    4. color = vector((n-i)/n, 0, i/n)

數據處理部分

若以 v0 = 30 m/s、b = 0、改變 theta 產生的數據如下 
v0(m/s), theta(degree), b, t(s), R(m)
30,15,0.0,1.5839999999999363,45.900795265255965
30,20,0.0,3.677999999999706,59.03149043776981
30,25,0.0,6.265000000000427,70.3385473519133
30,30,0.0,9.32600000000027,79.5271128295263
30,35,0.0,12.836999999998325,86.28128482496055
30,40,0.0,16.771999999997515,90.43154651019775
30,45,0.0,21.101000000002806,91.83195767269348
30,50,0.0,25.791000000008538,90.44021668289648
30,55,0.0,30.806000000014667,86.29457484901386
30,60,0.0,36.10800000000655,79.53000000000226
30,65,0.0,41.65599999999362,70.34058348411925
30,70,0.0,47.40899999998021,59.02925653658275
30,75,0.0,53.32199999996643,45.9119104107409
30,80,0.0,59.35099999995238,31.407745894621137
30,85,0.0,65.44999999994846,15.946886250519292
30,90,0.0,71.57199999997769,0.0
將數據匯入 SciDAVis 作圖得到的成果如下
t - 𝛳 關係圖
R - 𝛳 關係圖
若以 v0 = 30 m/s、theta = 30、改變 b 產生的數據如下 
v0(m/s), theta(degree), b, t(s), R(m)
30,30,0.0,3.060999999999774,79.5271128295263
30,30,0.1,5.980000000000332,65.7683163073318
30,30,0.2,8.780000000000573,55.694591159308345
30,30,0.30000000000000004,11.478999999999077,48.05540251489901
30,30,0.4,14.091999999997629,42.10134734984786
30,30,0.5,16.62899999999734,37.33313950127536
30,30,0.6000000000000001,19.09900000000036,33.44825688170242
30,30,0.7000000000000001,21.510000000003306,30.2339703048744
30,30,0.8,23.868000000006187,27.533728914252674
30,30,0.9,26.17800000000901,25.238123689064427
將數據匯入 SciDAVis 作圖得到的成果如下
t - b 關係圖
R - b 關係圖

結語

使用 VPython 可以避開很多不容易用手算的數學，例如在空氣阻力不可忽略的前提下，計算斜向抛射的飛行時間和水平射程。雖然是一次將一小段時間代入程式中運動，難免會有一點誤差，但只要將代入的時間縮短一點，計算結果就不致於偏差太多。

VPython官方說明書

  1. canvas: http://www.glowscript.org/docs/VPythonDocs/canvas.html
  2. box: http://www.glowscript.org/docs/VPythonDocs/box.html
  3. sphere: http://www.glowscript.org/docs/VPythonDocs/sphere.html
  4. math.radians: https://docs.python.org/3/library/math.html
HackMD 版本連結：https://hackmd.io/@yizhewang/HJ5pArGMm
於

沒有留言:

張貼留言