終端速度

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




這篇文章本來應該是前一篇文章〈自由落下〉最後一部分的內容，但由於這個程式的寫法比較不一樣，需要解釋的東西較多，所以獨立出來另外寫一篇。這個程式的目標：當小球從高空落下時，同時受到重力及空氣阻力的作用，試著找出小球的運動過程及終端速度，同時將得到的資料存成文字檔。

物理原理

假設空氣阻力

$$f = -bv$$

小球落下時所受合力（向下為正）

$$F = mg - bv = ma$$

小球剛開始運動時

$$ v = 0 ~~~~~ f = 0 ~~~~~ a = g $$

當小球速度增加時、增加、減小。當小球所受合力為零時，小球不會再加速，此時的速度稱為終端速度(terminal velocity)，通常代號為 v_t。

以下圖片是以小球質量 m = 1 kg、重力加速度 g = 9.8 m/s²、空氣阻力係數 b = 0.1 N s / m 模擬得到的結果，v_t = 97.99999000067568 m/s，理論值 = 98 m/s。

小球終端速度，b = 0.1， y-t 圖

小球終端速度，b = 0.1， v-t 圖

小球終端速度，b = 0.1， a-t 圖

程式 4-4：終端速度

取得程式碼

"""
 VPython教學: 4-4.自由落下, 有空氣阻力, 求終端速度
 Ver. 1: 2018/3/8
 Ver. 2: 2019/2/2
 Ver. 3: 2019/9/6
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size=1     # 小球半徑
m=1        # 小球質量
h=0        # 小球離地高度
g=9.8      # 重力加速度 9.8 m/s^2
b=0.1      # 空氣阻力 f=-bv
c1, c2=color.red, color.green
t=0        # 時間
dt=0.0001  # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Free Fall", width=400, height=400, x=0, y=0, center=vec(0, h/2, 0), background=vec(0, 0.6, 0.6))

# b1 with air drag, b2 without air drag 
b1 = sphere(pos=vec(2*size, h, 0), radius=size, color=c1, v=vec(0, 0, 0), a=vec(0, -g, 0))
b2 = sphere(pos=vec(-2*size, h, 0), radius=size, color=c2, v=vec(0, 0, 0), a=vec(0, -g, 0))

gd = graph(title="y-t plot", width=400, height=300, x=0, y=400, xtitle="t(s)", ytitle="y(m)")
gd2 = graph(title="v-t plot", width=400, height=300, x=0, y=700, xtitle="t(s)", ytitle="v(m/s)")
gd3 = graph(title="a-t plot", width=400, height=300, x=0, y=1000, xtitle="t(s)", ytitle="a(m/s<sup>2</sup>)")
yt1 = gcurve(graph=gd, color=c1)
yt2 = gcurve(graph=gd, color=c2)
vt1 = gcurve(graph=gd2, color=c1)
vt2 = gcurve(graph=gd2, color=c2)
at1 = gcurve(graph=gd3, color=c1)
at2 = gcurve(graph=gd3, color=c2)

# 開啟檔案 data.csv, 屬性為寫入, 先寫入欄位的標題
file = open("data.csv", "w", encoding="UTF-8")
file.write("t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n")
tp = 0

# 設定計算終端速度用的變數
eps = 0.0000000001
v1 = 0

"""
 3. 物體運動部分, 小球速度變化小於 eps 時停止運作
"""
while(True):
# 更新小球受力、加速度、速度、位置，畫 y-t 及 v-t 圖
    f = -b*b1.v
    b1.a = vec(0, -g, 0) + f/m
    b1.v += b1.a*dt
    b1.pos += b1.v*dt
    b2.v += b2.a*dt
    b2.pos += b2.v*dt
    yt1.plot(pos=(t, b1.pos.y))
    vt1.plot(pos=(t, b1.v.y))
    at1.plot(pos=(t, b1.a.y))
    yt2.plot(pos=(t, b2.pos.y))
    vt2.plot(pos=(t, b2.v.y))
    at2.plot(pos=(t, b2.a.y))
# 每隔 0.1 秒將資料轉成字串後寫入檔案
    tc = t
    if(tc == 0 or tc - tp >= 0.1):
        file.write(str(t) + "," + str(b1.pos.y) + "," + str(b2.pos.y) + "," + str(b1.v.y) + "," + str(b2.v.y) \
                   + "," + str(b1.a.y) + "," + str(b2.a.y) + "\n")
        tp = tc
# 計算終端速度, 當 abs(v2 - v1) > eps 時將 v2 的值指定給 v1, 當 abs(v2 - v1) > eps 時結束迴圈
    v2 = b1.v.y
    if(abs(v2 - v1) > eps): v1 = v2
    else: break
# 更新時間
    t += dt

print("t =", t, "vt =", v2)
file.close() # 關閉檔案

以下只說明這個程式與程式 4-3 的不同之處。

1. 為了使用 計算小球的加速度，需要定義小球質量 m。為了計算空氣阻力 f = -bv ，需要定義係數 b。在第55、56行計算空氣阻力 f 及加速度 b1.a。

f = -b*b1.v
b1.a = vec(0, -g, 0) + f/m

2. 第31~39行，由於 y、v、a 的數值差異很大，為了使 v 和 a 的曲線不會被壓扁，將 y-t、v-t、a-t 分開作圖。
3. 為了將資料存成文字檔，需要增加以下的程式碼，分別位於第42、43、70、71、81行。先用 open 開啟檔名為 data.csv 的文字檔，w代表寫入，若檔案不存在會新增檔案並寫入資料，若檔案已存在會覆寫檔案內容，最後的 encoding = "UTF-8" 是指定文字編碼格式為 UTF-8 (8-bit Unicode Transformation Format)。

file = open("data.csv", "w", encoding="UTF-8")

4. 再用 file.write 將字串 "t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n" 寫入檔案，最後的 \n 代表換行。

file.write("t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n")

5. 最後一定要用 file.close() 關閉檔案，否則程式會執行錯誤。

file.close() # 關閉檔案

6. 如果將所有的資料都寫入文字檔會使檔案太大，大約每隔 0.1 秒記錄一次資料應該就夠詳細了，因此在第41行定義了變數 tp、初始值為0，在第58行定義了變數 tc 並將此時的 t 數值指定給 tc。當 tc 等於 0 或 tc - tp >= 0.1（經過 0.1 秒）時，將 t、b1.pos.y、b1.v.y、b1.a.y、b2.pos.y、b2.v.y、b2.a.y 的值轉為字串並以逗號連接後寫入檔案。最後將 tc 的數值指定給 tp。

if(tc == 0 or tc - tp >= 0.1):
    file.write(str(t) + "," + str(b1.pos.y) + "," + str(b2.pos.y) + "," + str(b1.v.y) + "," + str(b2.v.y) \
             + "," + str(b1.a.y) + "," + str(b2.a.y) + "\n")
    tp = tc

7. 由於我只想趕快算出終端速度，並不想看小球落下過程的動畫，所以將 while 迴圈當中沒有設定 rate(1000) ，不限制每秒執行幾次，程式能跑多快就跑多快。
8. 為了計算終端速度，在第47、48行定義變數 eps = 0.0000000001、v1 = 0，在第74行定義變數 v2 = b1.v.y。若 abs(v2 - v1) > eps ，也就是速度變化仍大於我們所設定的精準度時，代表小球尚未達到終端速度，將 v2 的值指定給 v1，再執行一次 while 迴圈裡的程式碼。若條件不成立，也就是速度變化小於我們設定的精準度，代表小球已達到終端速度，使用指令 break 終止迴圈。

if(abs(v2 - v1) > eps): v1 = v2
else: break

9. 轉存成文字檔的資料格式如下（取得檔案）

---

t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)
0,-9.800000000000002e-08,-9.800000000000002e-08,-0.0009800000000000002,-0.0009800000000000002,-9.8,-9.8
0.10000000000000184,-0.04898368267428339,-0.049147097999999535,-0.9760913926464987,-0.9809799999999866,-9.702487885614207,-9.8
0.20009999999999428,-0.19518623057861884,-0.19649029399999726,-1.942460801550152,-1.9619599999999675,-9.60584997834477,-9.8
0.30019999999998326,-0.4376393325134774,-0.44202958800000336,-2.8992050587992453,-2.942940000000162,-9.510174595866035,-9.8
0.40029999999997223,-0.775384321771444,-0.78576498000003,-3.846420032023191,-3.9239200000003653,-9.415452151319196,-9.8
0.5003999999999612,-1.2074720800769563,-1.2276964700000768,-4.784200633998666,-4.904900000000568,-9.321673153331668,-9.8
0.6004999999999502,-1.732962942482813,-1.767824058000144,-5.71264083216006,-5.885880000000771,-9.228828205066046,-9.8
0.7005999999999392,-2.3509266032139267,-2.4061477440002315,-6.631833658015203,-6.866860000000974,-9.136908003278513,-9.8
0.8006999999999281,-3.0604420224488726,-3.1426675280003384,-7.541871216467292,-7.847840000001177,-9.045903337386646,-9.8
0.9007999999999171,-3.860597334029928,-3.9773834100004666,-8.442844695043977,-8.82882000000138,-8.955805088546489,-9.8

---

可以將資料匯入 LibreOffice Calc、MicroSoft Excel 之類的試算表軟體進行後續處理。

結語

在這個程式當中我們學到了將資料存成文字檔的方法，匯出的資料可以再利用其它的軟體處理並作圖。雖然在 Python 當中有一個很有名的繪圖套件 matplotlib，有興趣的同學可以參考〈Matplotlib 繪圖技巧：讀取數據及線性擬合〉、〈Matplotlib 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列〉。

VPython官方說明書

1. canvas: http://www.glowscript.org/docs/VPythonDocs/canvas.html
   
2. sphere: http://www.glowscript.org/docs/VPythonDocs/sphere.html
   
3. graph: http://www.glowscript.org/docs/VPythonDocs/graph.html
   

---

HackMD 版本連結：https://hackmd.io/@yizhewang/rJHr7-fG7