終端速度 - HackMD

<style> .markdown-body table{ display: unset; } </style> # 終端速度 > 作者：王一哲 > 第1版：2018年3月20日 > 第2版：2024年3月21日 <br /> 這篇文章本來應該是前一篇文章〈[自由落下](https://hackmd.io/s/S1e8LxzGQ)〉最後一部分的內容，但由於這個程式的寫法比較不一樣，需要解釋的東西較多，所以獨立出來另外寫一篇。這個程式的目標：當小球從高空落下時，同時受到重力及空氣阻力的作用，試著找出小球的運動過程及終端速度，同時將得到的資料存成文字檔。 <br /> ## 物理原理假設空氣阻力 $$f = -bv$$ 小球落下時所受合力（向下為正） $$F = mg - bv = ma$$ 小球剛開始運動時 $$ v = 0 ~~~~~ f = 0 ~~~~~ a = g $$ 當小球速度 $v$ 增加時、$f$ 增加、$a$ 減小。當小球所受合力為零時，小球不會再加速，此時的速度稱為終端速度 (terminal velocity)，通常代號為 $v_t$，由上式可知 $$ v_t = \frac{mg}{b} $$ 以下圖片是以小球質量 $m = 0.1 ~\mathrm{kg}$、重力加速度 $g = 9.8 ~\mathrm{m/s^2}$、空氣阻力係數 $b = 0.1 ~\mathrm{kg/s}$ 模擬得到的結果，$v_t = 9.799001457836292 ~\mathrm{m/s}$，$理論值 = 9.8 ~\mathrm{m/s}$。 <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/jZW99jv.png"> <div style="text-align:center">小球終端速度，b = 0.1， y-t 圖</div> <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/XE0pzKw.png"> <div style="text-align:center">小球終端速度，b = 0.1， v-t 圖</div> <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/6GkhYff.png"> <div style="text-align:center">小球終端速度，b = 0.1， a-t 圖</div> <br /> ## 程式 4-4：終端速度（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/04.%E8%87%AA%E7%94%B1%E8%90%BD%E4%B8%8B/4-4_airdrag_vt.py)） ```python= """ VPython教學: 4-4.自由落下, 有空氣阻力, 求終端速度 Ver. 1: 2018/3/8 Ver. 2: 2019/2/2 Ver. 3: 2019/9/6 Ver. 4: 2024/3/21 修改 while 迴圈的停止條件，修改參數減少程式運作時間作者: 王一哲 """ from vpython import * """ 1. 參數設定, 設定變數及初始值 """ r = 1 # 小球半徑 m = 0.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.001 # 時間間隔 """ 2. 畫面設定 """ scene = canvas(title="Terminal Velocity", 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*r, h, 0), radius=r, color=c1, v=vec(0, 0, 0), a=vec(0, -g, 0)) b2 = sphere(pos=vec(-2*r, h, 0), radius=r, color=c2, v=vec(0, 0, 0), a=vec(0, -g, 0)) gd = graph(title="<i>y</i> - <i>t</i> plot", width=600, height=450, x=0, y=400, xtitle="<i>t</i> (s)", ytitle="<i>y</i> (m)", fast=False) gd2 = graph(title="<i>v</i> - <i>t</i> plot", width=600, height=450, x=0, y=700, xtitle="<i>t</i> (s)", ytitle="<i>v</i> (m/s)", fast=False) gd3 = graph(title="<i>a</i> - <i>t</i> plot", width=600, height=450, x=0, y=1000, xtitle="<i>t</i> (s)", ytitle="<i>a</i> (m/s<sup>2</sup>)", fast=False) 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 = 1E-6 v1 = 0 v2 = -1E6 """ 3. 物體運動部分, 小球速度變化大於 eps 時繼續運作 """ while abs(v2 - v1) > eps: rate(1000) # 更新小球受力、加速度、速度、位置，畫 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 # 更新 v1 為 v2 目前的值，更新 v2 為 b1.v.y v1, v2 = v2, b1.v.y # 更新時間 t += dt print("t =", t, "vt =", v2) file.close() # 關閉檔案 ``` <br /> 以下只說明這個程式與程式 4-3 的不同之處。 1. 為了使用計算小球的加速度，需要定義小球質量 m。為了計算空氣阻力 f = -bv ，需要定義係數 b。在第60、61行 ```python f = -b*b1.v b1.a = vec(0, -g, 0) + f/m ``` 計算空氣阻力 f 及加速度 b1.a。 2. 第38~43行，由於 y、v、a 的數值差異很大，為了使 v 和 a 的曲線不會被壓扁，將 y-t、v-t、a-t 分開作圖。 3. 為了將資料存成文字檔，需要增加以下的程式碼，分別位於第46 ～ 48、73 ～ 78行。第46行，先用 open 開啟檔名為 data.csv 的文字檔，**w** 代表寫入，若檔案不存在會新增檔案並寫入資料，若檔案已存在會覆寫檔案內容，最後的 **encoding = "UTF-8"** 是指定文字編碼格式為 UTF-8 (8-bit Unicode Transformation Format)。 ```python file = open("data.csv", "w", encoding="UTF-8") ``` 第47行，用 **file.write** 將字串 "t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n" 寫入檔案，最後的 **\n 代表換行**。 ```python file.write("t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n") ``` 第85行，最後一定要用 **file.close() 關閉檔案**，否則程式會執行錯誤。我們在之後的課程〈[使用 For 迴圈計算水平抛射資料](https://hackmd.io/@yizhewang/SJoo2fGfQ)〉會介紹另一個寫法，使用 with 開啟檔案，使用這個寫法就不需要關閉檔案。 ```python file.close() # 關閉檔案 ``` 4. 如果將所有的資料都寫入文字檔會使檔案太大，大約每隔 0.1 秒記錄一次資料應該就夠詳細了，因此在第48行定義了變數 tp、初始值為0，在第73行定義了變數 tc，將此時的 t 數值指定給 tc。接著第74 ~ 77行，當 tc 等於 0 或 tc - tp >= 0.1（經過 0.1 秒）時，用 **str** 將格式為浮點數的變數 t、b1.pos.y、……等數值轉成字串，再用 **+** 接成一個較長的字串，最後用 file.write 寫入檔案。第78行，將 tc 的數值指定給 tp。 ```python 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 ``` 5. 如果只想要趕快算出終端速度，並不想看小球落下過程的動畫，可以刪除 while 迴圈當中的 rate(1000) ，不限制每秒執行幾次，程式能跑多快就跑多快。如果覺得刪除 rate(1000) 之後程式運算速度還是太慢，可以再刪除與繪製函數圖形有關的程式碼。 6. 為了計算終端速度，在第51 ~ 53行定義變數 eps = 1E-6、v1 = 0、v2 = -1E6。第57行，while 迴圈運作的條件為 ```python abs(v2 - v1) > eps ``` 也就是速度變化仍大於我們所設定的精準度時，代表小球尚未達到終端速度，需要繼續運算。接下來在第80行，更新 v1 為 v2 目前的值，更新 v2 為 b1.v.y。 ```python v1, v2 = v2, b1.v.y ``` 7. 轉存成文字檔的資料格式如下（[取得檔案](https://raw.githubusercontent.com/YiZheWangTw/VPythonTutorial/master/04.%E8%87%AA%E7%94%B1%E8%90%BD%E4%B8%8B/data.csv)）： ```csv= t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2) 0,-9.800000000000001e-06,-9.800000000000001e-06,-0.009800000000000001,-0.009800000000000001,-9.8,-9.8 0.10000000000000007,-0.048837982593961625,-0.0504798,-0.9419039213273656,-0.9898000000000009,-8.86696304171435,-9.8 0.20000000000000015,-0.18632103459137758,-0.19894980000000026,-1.7852642296382608,-1.969800000000004,-8.02275852889063,-9.8 0.3000000000000002,-0.40401831622843054,-0.44541979999999975,-2.548330013785355,-2.9497999999999855,-7.258928915129776,-9.8 0.4000000000000003,-0.6942928029302021,-0.7898897999999974,-3.23874594301281,-3.9297999999999664,-6.567821878866057,-9.8 0.5000000000000003,-1.050234574833993,-1.2323597999999951,-3.8634288540200266,-4.909799999999988,-5.942513659639614,-9.8 0.6000000000000004,-1.4655915907108468,-1.7728297999999951,-4.428637046335492,-5.8898000000000135,-5.376739693357867,-9.8 0.7000000000000005,-1.9347070527533745,-2.411299799999998,-4.940032980226851,-6.869800000000039,-4.864831851624774,-9.8 0.8000000000000006,-2.452462734727637,-3.147769800000003,-5.402740005277641,-7.849800000000064,-4.40166165637874,-9.8 0.9000000000000007,-3.0142277057300206,-3.982239800000011,-5.821393687957935,-8.82980000000009,-3.9825889009430098,-9.8 ... ``` 可以將資料匯入 LibreOffice Calc、MicroSoft Excel 之類的試算表軟體進行後續處理。 ## 結語在這個程式當中我們學到了將資料存成文字檔的方法，匯出的資料可以再利用其它的軟體處理並作圖。雖然在 Python 當中有一個很有名的繪圖套件 matplotlib，有興趣的同學可以參考〈[Matplotlib 繪圖技巧：讀取數據及線性擬合](https://hackmd.io/s/ByN63TEC7)〉、〈[Matplotlib 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列](https://hackmd.io/s/Bkd0EMyCm)〉。 <br /> ## 2022年9月29日補充：使用 NumPy 以及 Matplotlib 繪圖如果想要使用 Python 讀取文字檔的資料並繪圖，可以使用 NumPy 以及 Matplotlib，程式碼如下。 <br /> ```python= import numpy as np import matplotlib.pyplot as plt # 從 data.csv 讀取資料 t, y1, y2, v1, v2, a1, a2 = np.loadtxt('data.csv', skiprows=1, delimiter=',', unpack=True) # 繪圖 plt.figure(figsize=(8, 5), dpi=100) # 設定圖片尺寸 plt.xlabel(r'$t ~\mathrm{(s)}$', fontsize=16) # 設定坐標軸標籤 plt.ylabel(r'$v ~\mathrm{(m/s)}$', fontsize=16) plt.xticks(fontsize=12) # 設定坐標軸數字格式 plt.yticks(fontsize=12) plt.grid(color='0.5', linestyle='--', linewidth=1) # 設定格線顏色、種類、寬度 # 繪圖並設定線條顏色、寬度、圖例 line1, = plt.plot(t, v1, color='red', linewidth=3, label=r'$v_1$') line2, = plt.plot(t, v2, color='blue', linewidth=3, label=r'$v_2$') plt.legend(handles=[line1, line2], loc='lower left', fontsize=16) plt.savefig('v-t_plot.svg') # 儲存圖片 plt.savefig('v-t_plot.png') plt.show() # 顯示圖片 ``` <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/bQeNGnE.png"> <div style="text-align:center">使用 Matplotlib 繪製的 v-t 圖</div> <br /><br /> 假設有以下的實驗數據，若推測 y 與 x 成反比，可以使用以下的程式碼作圖並計算最接近直線的斜率及截距。 <div style="text-align:center"> | x | 0.2 | 0.4 | 0.6 | 0.8 | 1.0 | | - | --- | --- | --- | --- | --- | | y | 10.1| 4.9 | 3.3 | 2.6 | 1.9 | </div><br /><br /> ```python= import numpy as np import matplotlib.pyplot as plt x = np.asarray([0.2, 0.4, 0.6, 0.8, 1.0]) y = np.asarray([10.1, 4.9, 3.3, 2.6, 1.9]) z = 1/x # 產生計算最接近直線需要用到的陣列 A ，計算直線的斜率 m 和截距 c A = np.vstack([z, np.ones(len(z))]).T m, c = np.linalg.lstsq(A, y, rcond = -1)[0] print(m, c) plt.figure(figsize=(8, 6), dpi=100) # 設定圖片尺寸 plt.xlabel(r'$1/x$', fontsize=16) # 設定坐標軸標籤 plt.ylabel(r'$y$', fontsize=16) plt.xticks(fontsize=12) # 設定坐標軸數字格式 plt.yticks(fontsize=12) plt.grid(color='0.5', linestyle='--', linewidth=1) # 設定格線顏色、種類、寬度 # 繪圖並設定線條顏色、寬度、圖例 data, = plt.plot(z, y, color='blue', marker='o', markersize=10, linestyle='', label="Data") fitline, = plt.plot(z, m*z + c, color='orange', linewidth=2, label="Fitted Line") plt.legend(handles=[data, fitline], loc='upper left', fontsize=16) plt.savefig('fit_plot.svg') # 儲存圖片 plt.savefig('fit_plot.png') plt.show() # 顯示圖片 ``` <br /> 最接近直線的斜率、縱軸截距分別為 $slope = 2.028088701161563$、$intercept = -0.07080253431890154$，以下是繪圖成果。 <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/aQKfx36.png"> <div style="text-align:center">使用 Matplotlib 繪製資料點及最接近直線</div> <br /><br /> ## 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 --- ###### tags: `VPython`