宋家慶，鮑四元，吉 陳

(蘇州科技大學(xué)，江蘇 蘇州 215011)

空竹拋體運(yùn)動(dòng)的數(shù)值解法和動(dòng)畫模擬

宋家慶，鮑四元*，吉 陳

(蘇州科技大學(xué)，江蘇 蘇州 215011)

研究了馬格努斯效應(yīng)下的旋轉(zhuǎn)拋體運(yùn)動(dòng)，考慮了空氣阻力作用于質(zhì)心上的情況。對(duì)于不同初速度和拋出角度，基于 Mathematica 的計(jì)算功能計(jì)算得到運(yùn)動(dòng)方程的數(shù)值解，然后利用軟件的可視化功能畫出運(yùn)動(dòng)軌跡曲線。在Mathematica高版本中還通過(guò)控件改變參數(shù)，得到動(dòng)畫模擬。使學(xué)生更加深入、直觀地模擬空竹拋體的運(yùn)動(dòng)，同時(shí)提高了學(xué)習(xí)的主動(dòng)性和創(chuàng)造性。

馬格努斯效應(yīng)；旋轉(zhuǎn)拋體；計(jì)算機(jī)模擬；動(dòng)態(tài)可視化；Mathematica

抖空竹是漢族傳統(tǒng)文化中的一種游戲，空竹一般為木質(zhì)或竹質(zhì)，中空，因而得名，用線繩抖動(dòng)空竹可使其高速旋轉(zhuǎn)而發(fā)出的響聲。雖為玩具卻蘊(yùn)含著如下物理原理：物體旋轉(zhuǎn)時(shí)會(huì)產(chǎn)生馬格努斯效應(yīng)。足球運(yùn)動(dòng)中的“香蕉球”[1,2]也應(yīng)用了該原理，當(dāng)球員用右腳內(nèi)側(cè)“搓”球時(shí)引起摩擦，致使足球在球門附近時(shí)作逆時(shí)針旋轉(zhuǎn)，形成橫向作用力(即馬格努斯力)，使原本向右飛行的球逐漸向左偏轉(zhuǎn)。文獻(xiàn)[3-4]討論了在空氣阻力影響下空竹以不同初速度和拋出角度時(shí)的下落或上拋運(yùn)動(dòng)。在大學(xué)力學(xué)的教學(xué)和研究中，引入軟件進(jìn)行數(shù)值模擬可以達(dá)到事半功倍的效果[5，6]。故本文使用Mathematica軟件數(shù)值求解了空竹拋體運(yùn)動(dòng)的微分方程，繪出其軌跡。然后通過(guò)控件改變初速度及拋出角度的數(shù)值，動(dòng)態(tài)可視化模擬空竹的拋體運(yùn)動(dòng)。

1 旋轉(zhuǎn)拋體運(yùn)動(dòng)的理論計(jì)算

建立如圖1所示Oxyz直角坐標(biāo)系，Oxy平面垂直于地面，設(shè)旋轉(zhuǎn)拋體的自轉(zhuǎn)軸于Oz軸平行。

根據(jù)質(zhì)心運(yùn)動(dòng)定理[2]得旋轉(zhuǎn)拋體運(yùn)動(dòng)的微分方程如下：

(1)

其中m是拋體質(zhì)量， 是旋轉(zhuǎn)拋體的角速度，其矢量方向指向Oz軸正向，v是拋體質(zhì)心的速度，vx是Ox方向分速度， 是Oy方向分速度，k是空氣阻力對(duì)應(yīng)的系數(shù)， 是與流體性質(zhì)與物體大小幾何形狀有關(guān)的常數(shù)[3-4]。

圖1 空竹受力的示意圖

式(1)中的參數(shù)選取如下值：

μ=0.01，ω=12.56 rad/s，m=0.02 kg，g=9.8 m/s2，k=0.05 kg/s 。

微分方程式(1)的初始條件：t=0時(shí)，x=y=0，且初速度為v0，矢量方向與x軸成α角。

2 基于Mathematica的空竹運(yùn)動(dòng)軌跡

對(duì)于空竹給定的初速度和拋出角度可以運(yùn)用Do語(yǔ)句循環(huán)計(jì)算，運(yùn)用NDSolve命令求微分方程組數(shù)值解。建立微分方程組的語(yǔ)句如下：

mj=0.01;w=12.56;m=0.02;k=0.05;g=9.8；

v0=8(* 初速度 *)； (* 初速度的傾角 *);

dx[t_]=D[x[t],t];dx2[t_]=D[x[t],{t,2}](* 分別定義x[t]的一階、二階導(dǎo)數(shù) *);

dy[t_]=D[y[t],t];dy2[t_]=D[y[t],{t,2}](* 分別定義y[t]的一階、二階導(dǎo)數(shù) *);

eqs={-mj*w*dy[t]-k*dx[t]= =m*dx2[t],mj*w*dx[t]-k*dt[t-m*g= =m*dy2[t],

(dx[t]/.t->0)= =v0*Cos[ ],dy[t]/.t->0)= =v0*Sin[ ]}

對(duì)于上述常微分方程組，直接求解較為復(fù)雜。這里采用數(shù)值方法求解，在Mathematica中使用NDSolve命令。即輸入并執(zhí)行

subs=NDSolve[eqs,{x,y},{t,0,10}]

輸出為插值函數(shù)x[t]和y[t]。根據(jù)軌跡的參數(shù)方程,使用ParametricPlot命令可繪出xy平面上的軌跡，其中運(yùn)動(dòng)時(shí)間取為 秒，則輸入

ParametricPlot[Evaluate[{x[t]/.subs[[1]],y[t]/.subs[[1]]},{t,0,3},

AxesLable->{“x”,”y”}]；

運(yùn)行后可得對(duì)應(yīng)于輸入的初速度和拋出角度的運(yùn)動(dòng)軌跡，見(jiàn)圖2。

圖2 υ0=8 m/s和α=0.6rad時(shí)的軌跡圖

類似地，初速度和拋射角分別取若干組其他數(shù)值，可得出對(duì)應(yīng)的運(yùn)動(dòng)軌跡圖像，參見(jiàn)圖3-6，其中圖6中自轉(zhuǎn)角速度較大，取為ω=28.26 rad/s。

圖3 υ0=3.8 m/s和α=0.6rad時(shí)的軌跡圖

圖4 υ0=8 m/s和α=0時(shí)的軌跡圖

圖5 υ0=8 m/s和α=1.397 5 rad時(shí)的軌跡圖

圖6 υ0=6 m/s和α=0.785 rad時(shí)的軌跡圖其中ω=28.26 rad/s

上述圖3-圖6與文獻(xiàn)[3]中使用解析法所得的圖像一致，且有著較高的計(jì)算效率。利用文中方法可得任何初速度和拋出角度條件下的運(yùn)動(dòng)軌跡，避免了繁復(fù)的解析過(guò)程，且數(shù)值模擬有著足夠的精度，便于模擬運(yùn)動(dòng)軌跡。

進(jìn)一步，依照如下命令

ParametricPlot[Evaluate[{dx[t]/.subs[[1]],dy[t]/.subs[[1]]},{t,0,3},

AxesLable->{“vx”,”vy”}]；

可得運(yùn)動(dòng)過(guò)程的x方向速度和y方向速度的變化情況。由于以上多種情況的結(jié)果類似，這里只顯示一組結(jié)果，見(jiàn)圖7。另外，基于關(guān)系

(2)

類似地，可由軟件繪出空竹的速度-時(shí)間曲線，見(jiàn)圖8。由圖8可見(jiàn)，t=1.25s前后，空竹出現(xiàn)最小速度，大小約為1m/s。

圖7 υ0=8 m/s和α=0.6 rad時(shí)的x、y方向速度圖

圖8 υ0=8 m/s和α=0.6 rad時(shí)的速度-時(shí)間變化圖

3 空竹運(yùn)動(dòng)的動(dòng)畫模擬

在Mathematica的后期版本[5-8]中(如Mathematica9.0)，引入了控件命令，便于模擬動(dòng)畫，即使用Manipulate命令和Animator選項(xiàng)。

如果給出初速度和拋出角度的取值范圍，可使用Animator控件以滾動(dòng)條的形式改變初速度和拋出角度，從而實(shí)時(shí)模擬出空竹拋體運(yùn)動(dòng)的軌跡方程，代碼只需稍作修改如下：

Clear[v,];

Manipulate[ss=NDSolve[eqs,{x,y},{t,0,4}];

ParametricPlot[Evaluate[{x[t]/.First[ss],y[t]/.First[ss]}],{t,0,4}],

{v,0,10,Animator},{,0,Pi/2,Animator}]

運(yùn)行以上兩段代碼后則可生成一張隨初速度和拋出角度滾動(dòng)條控制的動(dòng)態(tài)圖像。截取其中兩幅如圖9和圖10所示。

圖9 實(shí)時(shí)動(dòng)態(tài)軌跡之一

圖10 實(shí)時(shí)動(dòng)態(tài)軌跡之二

上述方法可以實(shí)時(shí)顯示任意初速度和拋出角度下的運(yùn)動(dòng)軌跡，從而更加直觀的認(rèn)識(shí)空竹拋體運(yùn)動(dòng)的變化規(guī)律。

4 結(jié) 論

空竹斜拋運(yùn)動(dòng)方程的解析解的形式較復(fù)雜，需要大量的公式推導(dǎo)及運(yùn)算，而不少科學(xué)計(jì)算軟件具有數(shù)值求解微分方程的功能，且實(shí)現(xiàn)的過(guò)程簡(jiǎn)單，易于操作。利用了Mathematica軟件的強(qiáng)大數(shù)值計(jì)算能力及可編程擴(kuò)展功能實(shí)時(shí)生成軌跡和動(dòng)畫，從而通過(guò)軟件模擬再現(xiàn)了空竹的各類運(yùn)動(dòng)形式。

[1] 周道祥.從“野渡無(wú)人舟自橫”到“香蕉球”技術(shù)[J].力學(xué)與實(shí)踐,2005,27(3):94-95.

[2] 趙致真.從香蕉球說(shuō)開(kāi)去[J].力學(xué)與實(shí)踐,2008,30(3)：112-114.

[3] 郝成紅,黃耀清,王歡,等.考慮空氣阻力時(shí)空竹的斜拋運(yùn)動(dòng)[J].大學(xué)物理,2016(3): 15-17,26.

[4] 于鳳軍.馬格努斯效應(yīng)與空竹的下落運(yùn)動(dòng)[J].大學(xué)物理,2012,31(9):19-21.

[5] 張韻華.Mathematica7實(shí)用教程[M].合肥：中國(guó)科學(xué)技術(shù)大學(xué)出版社,2011.

[6] 鮑四元,孫洪泉,陳旭元.Mathematica在振動(dòng)波問(wèn)題中的應(yīng)用[J].物理與工程,2010,20(4):49-51.

[7]http://reference.wolfram.com/language/tutorial/IntroductionToManipulate.html.

[8] 張舜.關(guān)于物理實(shí)驗(yàn)在計(jì)算機(jī)上的應(yīng)用[J].大學(xué)物理實(shí)驗(yàn)，1999(4)：54-57.

A Numerical Solution for Motion of a Rotating Projectile and Its Animation Simulation

SONG Jia-qing，BAO Si-yuan，JI Chen

(Suzhou Science and Technology University，Jiangsu Suzhou 215011)

We studied the rotating projectile motion under the Magnus effect and considered the effect of the air resistance against the direction of movement.For different initial speeds and throwing angles,the numerical solutions of the motion equation are obtained based on the calculation function of Mathematica.Then we draw the trajectory curves based on the visual function of Mathematica.By changing the control parameter,the animation simulation is obtained.This can make students more deeply and directly simulate the diabolo projectile motion,at the same time,improve the learning initiative and creativity.

Magnus effect；rotating projectile；computer simulation；dynamic visualization；Mathematica

2016-06-24

1007-2934(2016)06-0056-04

O 4-39

A

10.14139/j.cnki.cn22-1228.2016.006.014

*通訊聯(lián)系人