基于計算流體力學方法的乒乓球軌跡仿真

余 萬, 李 春, 任 杰, 朱 玲, 施之皓, 季云峰

(1.上海理工大學 能源與動力工程學院,上海 200093; 2.上海體育學院 中國乒乓球學院,上海 200438)

?

基于計算流體力學方法的乒乓球軌跡仿真

余 萬1, 李 春1, 任 杰2, 朱 玲2, 施之皓2, 季云峰2

(1.上海理工大學 能源與動力工程學院,上海 200093; 2.上海體育學院 中國乒乓球學院,上海 200438)

采用計算流體力學(CFD)方法獲取乒乓球在不同運動狀態下的流場及軌跡,與運動學求解方法進行比較驗證CFD方法的優越性。結果顯示:基于CFD方法的計算結果比運動學求解方法更符合實際;旋轉方向對乒乓球軌跡有很大影響,上旋球軌跡偏低,下旋球軌跡偏高,側旋球軌跡向左向偏轉;旋轉速度對乒乓球軌跡也有很大影響,上旋球乒乓球旋轉速度越大其落點距離越近;乒乓球的旋轉運動將導致更為復雜的流場狀態。

乒乓球; 運動仿真; 軌跡; 流場; 計算流體力學

Author’s address 1.School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China; 2.China Table Tennis College, Shanghai University of Sport, Shanghai 200438, China

在球類運動中,相對于棒球、網球、排球、足球,乒乓球特性鮮明:質量輕(2.7 g)、旋轉高(通常3 000 r/min)以及運動速度快(通常5 m·s-1)[1]。由于球體輕,其自身旋轉對運動軌跡影響極大。多數有關乒乓球的研究都是基于乒乓球戰術的理論研究,或者對于球員打法與心理進行分析,以及對于乒乓球管理、教學進行調查等。乒乓球的運動軌跡受運動速度和旋轉速度及方向的影響,要求球員具有高水平、高精度的預判能力。隨著近年乒乓球機器人技術的發展,乒乓球運動軌跡預測受到許多研究人員的關注。乒乓球運動軌跡預測的準確性在一定程度上決定乒乓球機器人對擊球點預測的精確性,從而影響接球動作。

Nakashima等[1]基于高速攝像機提出一種乒乓球旋轉速度和運動速度的評估方法,并采用動態仿真方法,對乒乓球轉速及速度評估方法進行驗證。Liu等[2]基于乒乓球動態模型和碰撞模型,預測乒乓球擊球點,研究控制乒乓球機器人的手臂球拍運動方法,仿真求解非線性方程組驗證方法的有效性。Nakashima等[3]提出碰撞過程中滾動及滑動摩擦的判斷條件,分析球桌以及球拍的碰撞模型,并對模型的有效性和準確性進行驗證。Kui Ou等[4]采用計算流體力學(Computational Fluid Dyamics,CFD)方法對乒乓球自由軌跡仿真,計算乒乓球運動軌跡并分析其流場問題,通過球體和空氣的空氣動力學耦合仿真求解球體自由運動軌跡以及不同轉速對球體運動軌跡的影響。

孫在等[5]在二維平面中建立并求解乒乓球弧圈球運動方程,表明乒乓球轉速對乒乓球軌跡有直接影響,不同轉速的弧圈球軌跡存在明顯的差異。楊華等[6]基于ODE(Open Dynamics Engine)交互式可視化仿真環境,并借鑒網球阻力系數和升力系數,在三維空間對乒乓球旋轉球軌跡進行仿真研究。王奇志等[7]在OpenGL-3D仿真平臺中分析乒乓球運動軌跡中的運動學模型以及乒乓球碰撞模型,改進運動模型并采用最小二乘法曲線擬合求解,得出乒乓球預測軌跡。楊春卉等[8]基于動量理論和動量矩理論建立乒乓球上、下旋轉球碰撞數學模型,探討摩擦力對乒乓球碰撞影響;采用Matlab編程仿真,并通過實例驗證數學模型的有效性及準確性。

無論是近年開展乒乓球機器人研究比較活躍的日本名古屋大學Nakashima教授,還是國內的一些學者,在分析乒乓球受力及運動情況時:一方面,將乒乓球運動中瑪格努斯升力系數及阻力系數通常視為常數或采用網球及其他類球體計算公式[5-6],在球體實際運動過程中,其瑪格努斯升力系數及阻力系數與旋轉速度和運動速度有關;另一方面,未考慮乒乓球在運動過程中受到空氣摩擦阻力及其他力作用,旋轉速度會發生一定的變化,瑪格努斯升力系數及阻力系數也隨之變化,從而導致軌跡預測不準確。

1985年,美國國家航天航空局(NASA)研究中心空氣動力學研究所科研人員研究了除乒乓球之外其他球體(如棒球、網球)的空氣動力學特性[8-9]。乒乓球作為我國“國球”,更應從空氣動力學方面對其運動和空氣動力學特性進行研究。因此,本文采用CFD方法無網格算法對乒乓球在不同轉速大小方向運動進行仿真計算,對比分析不同方法對于乒乓球運動軌跡的影響。

1 乒乓球理論基礎

1.1 乒乓球及球桌幾何參數 國際乒聯規定比賽使用乒乓球質量為2.7 g以及半徑為0.02 m,并規定恢復系數為0.89～0.92。乒乓球桌尺寸:長2.74 m、寬1.525 m、高0.76 m以及網高0.152 5 m[10-11]。

1.2 乒乓球旋轉及瑪格努斯效應

1.2.1 乒乓球旋轉球 乒乓球在運動過程中,其運動狀態各不相同,運動速度以及旋轉速度各不相同,應了解乒乓球旋轉速度大小與方向性。如圖1所示,乒乓球旋轉通常分為左旋、側旋、上旋、下旋、順旋以及逆旋,但在乒乓球運動過程中的旋轉通常是以上六大旋轉中多種旋轉綜合作用[12]。

左右旋:乒乓球繞著通過球心并垂直球桌面的軸旋轉的球稱作側旋球,即左右旋。當以發球員為參考基準,乒乓球向左旋轉稱作左旋球,向右旋轉稱作右旋球,如圖1(a)(d)所示。上下旋:乒乓球繞著通過球心垂直運動方向并平行球桌面的軸的旋轉,以發球員為參考基準,乒乓球繞該軸向前旋轉稱作上旋球,繞該軸向后旋轉稱作下旋球,如圖1(b)(e)所示。順逆旋:乒乓球繞著通過球心與球的運動方向平行的軸的旋轉,以發球員作為參考基準,乒乓球繞軸順時針旋轉稱作順旋球,繞軸逆時針旋轉時稱作逆旋球,如圖1(c)(f)所示。

注:v為乒乓球運動速度;ω為乒乓球旋轉角速度

1.2.2 瑪格努斯效應 在球體運動中,乒乓球的旋轉球、排球的側旋球以及足球的“香蕉球”有著相似的形成原理,即旋轉球在運動過程其運動軌道產生彎曲,這種現象稱作瑪格努斯效應,即旋轉球在運動過程因旋轉受到力的作用使其實際軌跡偏離其原有軌跡的現象,所受的力稱作瑪格努斯力。牛頓早在1671年就解釋了旋轉羽毛球運動軌跡偏離現象。瑪格努斯1852年實驗研究了旋轉圓柱體的瑪格努斯效應[13]。

圖2表明,流體(空氣)經過旋轉球下方因球體旋轉速度方向與流體速度方向相反,因球體表面黏性作用使得流體速度減少,在其上方旋轉速度方向與流體速度方向相同使得流體速度增大。根據伯努利方程[14]可知:球體下方壓力大于上方壓力,從而受到向上的作用力,該作用力大小以及作用力方向即瑪格努斯力大小和方向。

圖2 瑪格努斯效應示意

1.3 乒乓球受力分析 乒乓球在運動過程中受到重力、浮力、阻力以及瑪格努斯力等一系列力的綜合作用。重力大小與乒乓球質量有關,方向為豎直向下。浮力大小與乒乓球體積有關,作用力方向與重力向反[1]。阻力大小與乒乓球速度有關,阻力方向與乒乓球運動速度方向相反?，敻衽沽Υ笮∨c旋轉速度以及運動速度大小有關,其作用力方向與旋轉速度和運動速度方向有關。圖3為下旋和側旋乒乓球在運動中所受的作用力示意。

圖3 不同旋轉乒乓球受力情況

2 CFD無網格理論基礎

無網格方法分為:分子水平仿真,如Direct Simulation Montecarlo;宏觀層面仿真,如Vortex Particle Method和Smoothed Particle Hydrodynamics;基于介觀層面上的仿真,如Lattice Gas Automata和Lattice Boltzmann Method[15]。本文采用基于介觀層面的格子Boltzmann方法,其方法基本思路是以Boltzmann方法為基礎,通過碰撞遷移獲得宏觀流體的基本信息,該方法比微觀尺度的直接數值模擬節省時間,比宏觀尺度的粒子法具有更高的精度。

2.1 格子Boltzmann方法 格子Boltzmann方法(Lattice Boltzmann Method,LBM)是相對新的復雜流體仿真技術,并吸引了許多研究人員的關注。傳統CFD方法是宏觀性質上方程守恒(質量、動量以及能量)數值求解,而LBM是假想流體通過粒子組成,并且這些粒子在一個離散的格子執行連續的繁殖和碰撞過程。因為這些粒子的性質從而使得LBM比傳統CFD方法有著一定的優勢,尤其是在復雜的邊界處理、結合微觀相互作用以及計算的并行化方面。LBM來源于格子氣自動機(Lattice Gas Automata,LGA)方法,該方法被認為是一個簡化的虛擬分子動力學模型,其空間、時間以及粒子速度都是離散的[16]。LGA在用于流體力學仿真模擬時存在一些缺陷:統計噪聲和低雷諾數縮放時格子大小[17-19]。從LGA到LBM的轉變主要是通過使用總體平均函數,即密度分布函數替換格子方向布爾粒子數以除去統計噪聲。在LBM發展中一個重要的BGK(Bhatnagar-Gross-Krook)近似簡化,LBGK(Lattice Bhatnagar-Gross-Krook)模型的使用使得LBM方法在仿真時更加高效。

2.2 Boltzmann方程及BGK近似

2.2.1 Boltzmann方程 Boltzmann方程為:

(1)

該方程求解是很困難的,僅在特殊情況下才能求解。Boltzmann求解難點也就是方程等式右邊的Ω碰撞項,該項同時是非線性并且還與分子間的相互作用力存在耦合作用[21]。

2.2.2 BGK近似 為簡化求解Boltzmann方程,Mohamad等[20-23]進而引入如下方程替代碰撞項,從而得到Boltzmann方程的BGK近似方程:

(2)

平衡分布函數為[19]:

O(u3)+ωiρ

(3)

式中:cs為格子聲速;u為宏觀速度;ωi為權系數;ρ為流體宏觀密度;O(u3)為無窮小量。

2.3 格子Boltzmann方程 格子Boltzmann方程是Boltzmann方程基于BGK近似在空間、速度以及時間上的離散得來的[22-23]?；贐GK近似的Boltzmann方程在空間、速度以及時間上的離散,通過數學轉換可以得到如下含有外力項的格子Boltzmann方程[20]:

fi(r+eiδt,t+δt)-fi(r,t)=

(4)

式中:ei為空間離散量;δt為時間離散步長;Fi(r,t)為外力項。

本文采用簡化外力項格子Boltzmann方程,即在仿真計算時不考慮外力項。簡化的格子Boltzmann方程如下所示。

(5)

3 乒乓球運動仿真結果與分析

3.1 方法學對比 根據文獻[5]中乒乓球初始參數,采用計算流體力學方法求解乒乓球運動軌跡并與文獻[5]采用運動學方法求解得到的運動軌跡對比。圖4(a)(b)(c)(d)分別為運動軌跡對比及乒乓球各坐標軸方向旋轉速度變化曲線對比。球桌長度方向為X軸方向,寬度方向為Z軸方向,高度方向為Y軸方向。

圖4(a)為運動學求解軌跡和計算流體力學求解軌跡對比分析。計算流體力學求解得到的軌跡在Y軸和X軸均大于運動學求解得到的軌跡,高度方向約高0.05m,相對誤差約為25%,長度方向長約0.25m,相對誤差約為8.5%。文獻[5]中簡單地考慮乒乓球升力系數為1,阻力系數為0.44,相比文獻[4]中對于乒乓球靜止和旋轉時升阻力系數,文獻[5]中對于升阻力系數設置不合理。在圖4(b)(c)(d)中對比了運動學方法和計算流體力學求解過程中乒乓球旋轉速度隨時間變化曲線,運動學求解過程中3個方向的旋轉速度均為定值,而在計算流體力學求解中考慮空氣阻力作用,乒乓球在運動過程中受到阻力作用使得旋轉速度減小,在與球桌面發生碰撞各方向的旋轉速度均發生一定的變化。綜上,乒乓球運動軌跡的計算流體力學求解比運動學求解更為準確。

3.2 計算流體力學結果分析

3.2.1 條件設置 設置計算區域長3m、寬2m、高1m;對于乒乓球真實環境建模,乒乓球、球桌以及球網均按照國標所給的幾何參數建模以及求解尺度如圖5所示;外圍流場求解尺度為0.1m,尾跡加密尺度為12.5mm,乒乓球邊界求解尺度為6.25mm,球桌以及球網表面求解尺度為12.5mm。仿真計算條件:密度ρ=1.225 kg·m-3,動力黏度μ=178.94 μPa·s-1,D為乒乓球直徑(40 mm);乒乓球質量為2.7 g,重力加速度為9.81 m·s-2。乒乓球定義初速度設定:X軸方向分速度為6.2 m·s-1,Y軸方向分速度為1.6 m·s-1。

圖4 運動學及計算流體力學方法計算結果

Figure 4. The results of methods of kinematics and computational fluid dynamics

分別對無旋、下旋、側旋及上旋乒乓球運動進行仿真,并對不同旋轉速度的乒乓球運動進行仿真(表1)。

圖5 乒乓球、球桌、球網建模及求解尺度

表1 乒乓球旋轉方向及對應的旋轉速度

3.2.2 軌跡分析 由圖6可知,3種不同旋轉速度使得乒乓球軌跡落點距離隨著旋轉速度越快變得越小。由于旋轉速度越快旋轉產生的瑪格努斯力越大,軌跡向下偏轉越大,落點越近。下旋球和側旋球,旋轉速度越快,所受瑪格努斯力越大,下旋球軌跡越高,側旋球軌跡及落點越偏左。

圖6 不同旋轉速度上旋球XY平面軌跡

Figure 6. The topspin trajectory in different rotation speeds in face ofXY

圖7為分析不同旋轉方向在相同運動速度和旋轉速度條件下乒乓球運動軌跡。相比于無旋轉乒乓球運動軌跡,下旋、側旋以及上旋球旋轉效應受到向上、向左以及向下的瑪格努斯力作用。這種作用力使得下旋球軌跡向上偏轉在桌面上無落點,側旋球軌跡偏左落點偏左,但對乒乓球在X軸方向的落點影響不大;上旋球軌跡向下偏轉使得軌跡落點距離變小。

3.2.3 流場分析 圖8表明,在乒乓球運動過程中,在乒乓球后方造成速度最大區域,對圖9該區域壓力最低;上旋及下旋球相對于無旋球,乒乓球后高速區域分別向上及向下移動;在t=0.16 s時可看出,無旋轉乒乓球趨近軌跡最高點,下旋則處在軌跡上升段,上旋處在軌跡下降段;側旋乒乓球速度流場中高速區域受旋轉效應影響發生偏移;在與桌面發生碰撞即乒乓球落點處,側旋球因受側向旋轉作用與其他2種情況下落點速度流場范圍小。根據圖9壓力流場可知旋轉使得乒乓球周圍流暢變得更為復雜,無旋轉球的壓力流場趨近對稱,而具有旋轉的乒乓球壓力流場不具有對稱性;高壓區處于乒乓球運動前端,低壓區處于球運動尾部。

圖7 相同旋轉速度不同旋轉方向的乒乓球運動軌跡

Figure 7. The trajectory of table tennis in the same speed but different rotation directions

圖8 不同旋轉方向的乒乓球速度流場

Figure 8. The speed flow fields of table tennis in different rotation directions

圖9 不同旋轉方向的乒乓球壓力流場

Figure 9. The pressure flow fields of table tennis in different rotation directions

4 結論

(1) 對比計算乒乓球運動軌跡的2種方法發現,相較于運動學求解方法,計算流體力學計算結果更符合實際。

(2) 分析計算流體力學求解乒乓球運動軌跡可知,不同旋轉方向對乒乓球球軌跡有很大影響:乒乓球在向下旋轉時,乒乓球的旋轉產生的瑪格努斯力使得軌跡偏高,軌跡落點超出球桌桌面區域;乒乓球在向側旋轉時,乒乓球軌跡向左側偏轉,軌跡落點亦偏左;乒乓球在向上旋轉時,乒乓球受到向下的瑪格努斯力,該作用力使得乒乓球軌跡偏低,落點距離小。旋轉速度對乒乓球軌跡亦有很大影響:對于向上旋轉的乒乓球,旋轉速度越大,乒乓球軌跡落點距離越小。

(3) 對乒乓球運動流場分析可知,相比于無旋轉乒乓球,旋轉乒乓球流場較為復雜,無旋轉乒乓球速度及壓力流場具有對稱性,旋轉乒乓球因受旋轉作用其流場不具有對稱性。

[1] Nakashima A,Okamoto T,Hayakawa Y.An online estimation of rotational velocity of flying ball via aerodynamics [C]∥IFAC World Congress.South Africa:Cape Town,2014:7176-7181

[2] Liu C,Hayakawa Y,Nakashima A.Racket control and its experiments for robot playing table tennis [C]∥IEEE International Conference on Robotics and Biomimetics.China:Guangzhou,2012:241-246

[3] Nakashima A,Kobayashi Y,Ogawa Y,et al.Modeling of rebound phenomenon between ball and tacket tubber with spinning effect [C]∥ICCAS-SICE 2009 - ICROS-SICE International Joint Conference.Japan:Fukuoka,2009:2292-2300

[4] Kui Ou,Patrice C,Antony J.Computational sports aerodynamics of a moving sphere:Simulating a ping pong ball in free flight [C]∥Aiaa Applied Aerodynamics Conference.Hawaii:Honolulu,2011:1-16

[5] 孫在,余廣鑫,郭美,等.乒乓球弧圈球的空氣動力學原理及其運動軌跡的仿真分析 [J].體育科學,2008,28(4):69-71

[6] 楊華,關志明.基于ODE的乒乓球運動軌跡仿真研究 [J].計算機仿真,2011,28(9):230-233

[7] 王奇志,楊曉曉.乒乓球軌跡預測的研究與仿真 [J].計算機工程與科學,2013,35(2):164-168

[8] 楊春卉,袁志華,梁振剛.乒乓球反彈動態特性的仿真研究 [J].計算機仿真,2014,31(10):281-285

[9] Mehta R D.Aerodynamics of sports balls [J].Annual Review of Fluid Mechanics,1985,17(1):151-189

[10] The international table tennis federation.ITTF Handbook 2011/2012 [EB/OL].[2016-08-07].http://www.ittf.com/ittf_handbook/ittf_hb.html

[11] The international table tennis federation.ITTF Technical Leaflet T1:The Table [EB/OL].[2016-08-07].http://ittf.com/stories/pictures/T1_The_Table_BoD2013.pdf

[12] 王正倫.教你打乒乓球[M].南京:江蘇科學技術出版社,2012:1-163

[13] Seifert J.A review of the magnus effect in aeronautics [J].Progress in Aerospace Sciences,2012,55(5):17-45

[14] 潘慧炬.瑪格努斯效應的力學模型 [J].浙江體育科學,1995,17(3):16-19

[15] 劉桂榮.光滑粒子流體動力學 [M].長沙:湖南大學出版社,2005:1-148

[16] Mcnamara G R,Zanetti G.Use of the boltzmann equation to simulate lattice-gas automata [J].Physical Review Letters,1988,61(20):2332-2335

[17] Higuera F J,Jiménez J.Boltzmann approach to lattice gas simulations [J].Europhysics Letters,1989,9(7):663

[18] Chen S,Doolen G D.Lattice boltzmann method for fluid flows [J].Annual Reviews of Fluid Mechanics,2003,30(1):329-364

[19] Succi S.The lattice boltzmann equation:For fluid dynamics and beyond [M].Oxford:Oxford University Press,2001:1-213

[20] Mohamad A A.Lattice boltzmann method:Fundamentals and engineering applications with computer codes [M].Berlin:Springer,2011:1-324

[21] He X,Luo L.A Priori derivation of the lattice boltzmann equation [J].Physical Review E Statistical Physics Plasmas Fluids & Related Interdisciplinary Topics,1997,55(6):6333-6336

[22] He X,Luo L S.Theory of the lattice boltzmann method:From the boltzmann equation to the lattice boltzmann equation [J].Physical Review E:Statistical Physics,Plasmas,Fluids,and Related Interdisciplinary Topics,1997,56(6):6811-6817

[23] Qian Y H,Humieres D,Lallemand P.Lattice BGK models for navier-stokes equation [J].Europhysics Letters,1992,17(6):479

Simulation of Table Tennis Trajectory Based on Computational Fluid Dynamic Method

YU Wan1, LI Chun1, REN Jie2, ZHU Ling2, SHI Zhihao2, JI Yunfeng2

The trajectories of table tennis under different motions are obtained based on Computational Fluid Dynamics (CFD) method. Compared with the kinematics method, the superiority of CFD method is verified. The study analyzes the trajectories and flow fields of the table tennis on different rotational speeds in different directions and draws the following conclusions: The CFD method is more credible than kinematic method. The rotation directions have a significant impact on the trajectory, with a low trajectory for topspins, high trajectory for backspins, and a leftward trajectory for side spin. Rotating speed also has a great impact on the trajectory of table tennis: the greater the rotation speed of the topspin table tennis is, the closer the landing point is. The rotational movement of table tennis will lead to more complicated flow field state.

table tennis; motion simulation; trajectory; flow field; computational fluid dynamics

2016-09-20；

2016-11-25

國家自然科學基金資助項目(51676131,51176129);上海市教育發展基金“晨光計劃”項目(13CG55)

余萬(1994-),男,江西九江人,上海理工大學碩士研究生;Tel.:15821303829,E-mail:15821303829@163.com

李春(1963-),男,北京人,上海理工大學教授,博士生導師;Tel.:(021)55271729,E-mail:lichunusst@163.com

G847

A

1000-5498(2017)03-0089-06

DOI 10.16099/j.sus.2017.03.014