落體的運動

 

一、電腦模擬實驗:自由落體問題

問題:想知道炮彈打到那

物理量:軌跡隨時間的變化 x(t), y(t)

物理定律(數學模型):f = m a ,在此 f 純粹來自於重力 (粗體字代表是向量)

粒子在地表附近所受的重力是 g (g = GmM/R),g = 9.8 牛頓/米

f = ma = m ( ax ex + ay ey ) ,其中 ax = dvx/dt 、 ay = dvy/dt

只有重力時,f = fy ey = g m ey ,也就是說 fy = g m

演算法:

利用 Euler 演算法

位置 xn+1 = xn + vx,n Δt ; yn+1 = yn + vy,nΔt

速度 vx,n+1 = vx,n + ax,n Δt ; vy,n+1 = vy,n + ay,nΔt

 

進階補充:對計算穩定性及精確度,有更好的選擇

Euler-Richarson 演算法, 先求出

ax,n = fx(rn,vn,t)/m

vx,mid = vx,n + (1/2) ax,nΔt

xmid = xn + (1/2) vx,nΔt

ax,mid = fx( xmid, vx,mid, t+(1/2)Δt ) / m

然後才做

vx,n+1 = vx,n + ax,midΔt

xn+1 = xn + vx,midΔt

註:為什麼 Euler-Richarson 演算法 比 Euler 的好,可參考 Gould 課本上的推導: p1p2

 

程式流程:

(0) 開始,宣告變數

(1) 問質量、初速度(含初速率與角度)

(2) 建立位置 x 與 y、速度 vx 與 vy 的初始值

(3) 用演算法求下一個小時段 Δt 後的位置與速度

(4) 畫出粒子位置並檢查高度是否己達地面,未達地面則回到 (3)

(5) 問要不要再射一次,若要則回到 (1)

(6) 結束

 

你等不及了嗎?先偷玩一下老師己經做好的可執行檔 cannon.x

你做不出來嗎?先偷看一下老師己經寫好的範圍程式 cannon.f cannon.f.txt

非不得已請勿先看

 

動動腦、動動手,自己再改造(或回家做):

1. 畫個大砲把炮管的角度表現出來

2. 以火藥包裝填數表現炮彈動能,轉換成初速率

3. 物體被炸破片的彈出(配合亂數產生器 ran3.f使用方法

4. 試寫一個煙火程式,會多段開花變色的

 

分析與討論:

想預測炮彈打到那堙A主要是基於那一條物理定律?

什麼樣的先決(初始)條件決定了炮彈的軌跡(和落點)?

同樣的火藥包充填數,怎麼打炮彈才飛得最遠?

 

 

二、電腦模擬實驗:受空氣阻力的落體

問題:在有空氣阻力的情況下,炮彈軌跡變得如何?

物理量:一樣是軌跡隨時間的變化 x(t), y(t)

物理定律(數學模型):f = m a,受力 f 同時包含了重力及空氣阻力

空氣阻力非常複雜(因此飛機或高速火車的阻力設計需要靠風洞實驗或電腦模擬來進行),基本上它是與速率有關,速率越高阻力越大,靜止的物體則沒有空氣阻力。可見它是一個與速率有關的量,常見的簡化公式,有隨速率一次方呈正比的,也有隨兩次方呈正比的,方向是速度反向(負號)的

Fd = - k1 v

Fd' = - k2 v2

由於有 x 與 y 兩個分量,假設速度與水平線的夾角為 θ,即 θ = tan-1(vy/vx) [ Fortran 程式寫成 theta = atan(vy/vx) ] ,

Fd,x = Fd cosθ

Fd,y = Fd sinθ

(或者是,直接用 vx/√(vx2+vy2) 及 vy/√(vx2+vy2) 會比先求 θ = atan(vx/vy) 取角度才去求 cos(θ) 及 sin(θ) 更好)

因此,在有關力之部分的公式,

fx = Fd cosθ

fy = - m g + Fd sinθ

其他的部分,與上一個議題相同

新議題:如何求得室氣阻力係數 k?

利用終端速度,即垂直落下過程中重力與阻力達平衡,使合力為零不再加速,而呈等速運動時的速度。

合力 = 0 = - m g + Fd = - m g + k1 vy2

即 k1 = vy2 / (m g),此終端速度可藉由實驗量得

一個粉筆頭大小的石子其終端速度大約是每秒三十公尺。另外,一個重 0.254 克 、半徑 2.54 公分的保利龍球有以下的實測數據(引用 Physics Teacher 24, 153 (1986) )

時刻(秒) 位置(公尺)
-0.132 0.0
0.0 0.075
0.1 0.260
0.2 0.525
0.3 0.870
0.4 1.27
0.5 1.73
0.6 2.23
0.7 2.77
0.8 3.35

想想看,我們要怎麼從上列的數據求得空氣阻力係數 k?

 

老師寫的範例程式:styrofoam_k.f.txt styrofoam_k.f styrofoam_k.x

上面表格數據的資料檔:styrofoam.txt

儘量自己想,不會再參考

演算法:

同樣採用 Euler 或 Euler-Richardson 演算法

 

程式設計:

與上一個議題相同

參考範例程式

cannon_drift.f cannon_drift.f.txt cannon_drift.x

 

質量與空氣阻力係數對拋射體軌跡的影響

mass = 1.0, k = 2.0

mass = 0.5, k = 1.0

mass = 2.0, k = 2.0

 

 

 

分析與討論:

以同樣的力道向上拋一個小石子,在有空氣阻力與無室空氣阻力的兩種情況比較下,那一種會在空中停留較久?

那一種軌跡像羽毛球?那一種像鉛球?

在有空氣阻力的情況下,拋射物的角度要如何設定其距離才會最遠?

 

 

挑戰:

試著做出你自己的憤怒鳥遊戲 (技術上的細節歡迎與教授討論 ) 。

 

 

作業:彈珠台程式

彈珠與光滑釘子的交互作用模式

N = Nx ex + Ny ey

| N | = 1

V dot N = Vx * Nx + Vy * Ny

碰撞後 V 平行於 N 上的分量因反彈而方向相反,與 N 垂直的分量則繼續維持,因此碰撞後

V = (-1)* (V dot N) N

V// = V - (V dot N) N = V + V

新的速度則為 V// + V

 

老師提供的的範例程式 pinball.f pinball.f.txt pinball.x ,輸入 範例

 

動手改造與討論:

將你的程式改為一次釋放兩顆、三顆彈珠。

如果是對彈珠有摩擦力的釘子,則彈珠與之交互作用模式如何?

 

 

 

PGPLOT: PGPT (n,x,y,SYMBOL)

 

 

亂數產生器的用法

ran3.f 內含的 ran3(iseed) 函式為例:

iseed = 7
x = ran3(iseed)
y = ran3(iseed)

準備任意一個整數值放在某個整數變數 iseed 的值內,依上述的方法每次使用 ran3( ) ,則 x 與 y 都會獲得一個 0.0 到 1.0 之間的數值,每次的數值雖然不一樣,但這些數值在 0 到 1 之間出現的機率是均等的。另外,直接把整數填入 ran3( ) 卻並不支援,如 xx = ran3(7),會導致 "區段錯誤 (Sgementation Fault)"

試試看,寫一個丟骰子程式。例範:dice.f ,記得與 ran3.f 一起編譯,如 gfortran diec.f ran3.f <Enter>

另外 gfortran 支援內建副程式 random_number(實數變數),就不必另外用外部的亂數函式了,可參考另一骰子程式範例 dice_gfortran.f。重設亂數種子的副程式是 random_seed 細節請見: http://www.nsc.liu.se/~boein/f77to90/a5.html#section21c

 

進階閱讀:

Gould and Tobochnik, An Introduction to Computer Simulation Methods -- Applications to Physical Systems, Addison Wesley (1996) Chapter 2, Chapter 3