基于无网格方法的环境温度对乒乓球运动轨迹影响研究

丁勤卫，李 春，张俊伟，任 杰，施之皓，郝文星

乒乓球作为一项极具观赏性和竞技难度的运动，具有质量轻、速度快、旋转高以及战术多变性等特点[1-3]。乒乓球运动已成为我国普遍性运动且具有极为广泛的群众基础，并作为我国国球，继1996年奥特兰大奥运会以来已多次包揽奥运会冠军[2]。国内外学者对于乒乓球做了大量研究[4]，其中多数研究为战术理论研究[5]、球员打法研究[6-7]以及心理分析、乒乓球管理和教学管理研究[8]。

计算机技术及数值计算发展，为乒乓球运动学及动力学数值仿真提供了技术保证。例如:计算机仿真技术在乒乓球拍的强度研究、乒乓球碰撞研究以及乒乓球运动轨迹预测研究的应用等[9-13]。文献[9]基于ODE的交互式可视化仿真环境，着重于旋转球，考虑旋转过程中受到重力、Magnus升力以及空气阻力，借鉴网球升阻力系数，不考虑旋转速度的变化，求解运动方程得到乒乓球运动轨迹，验证了运动模型的准确性，为揭示乒乓球运动轨迹、乒乓球机器人的发展提供了科学依据。文献[10]基于模型参数学习和自适应模型的调整的乒乓球运动状态估计和轨迹预测算法，实验验证了算法的实时性和高效性。文献[11]采用改进的运动学模型的算法，考虑空气阻力对乒乓球速度的影响，将其速度作为动态反馈参数，结果表明效果准确性，并采用OPENGL的3D仿真平台直观分析预测轨迹。A NAKASHIMA等[12]实验研究了旋转乒乓球与球拍的碰撞模型，通过水平和垂直方向恢复力系数确定碰撞后乒乓球速度变化，并通过了其他实验数据验证碰撞模型和各方向的恢复力系数的可行性。

然而，乒乓球在空气中运动，实质上是一种流体的绕流过程，空气的粘性会影响着球体的升力和阻力，进而影响乒乓球的速度、旋转、轨迹及落点。在上述乒乓球轨迹预测仿真研究中为便于求解或因计算资源限制，多将乒乓球运动过程中的升、阻力系数作为常数以及乒乓球旋转速度作为定值。显然这与实际情况存在较大差距。为更加真实地研究乒乓球复杂的气动流场，OU KUI等[14]采用计算流体力学方法（CFD）对乒乓球轨迹进行数值模拟，得到乒乓球运动轨迹并分析其流场，仅求解了旋转速度大小及方向对乒乓球轨迹影响。

根据专业教练带队经验以及自身体验，乒乓球作为球员之间的快速回往复运动，不仅需要考虑乒乓球运动速度及旋转大小方向对乒乓球运动轨迹的影响，更需要研究环境因素对乒乓球运动轨迹影响。为此必须对乒乓球运动的三维非定常流场开展数值模拟，以便从流动机理角度研究球体不同运动状态的内在机理。1985年，R D MEHTA[15]提出环境温度对棒球、足球及网球的运动存在影响，但未涉及乒乓球。显然，轻而小的乒乓球比棒、足、网球轻得多的乒乓球受温度影响更大。因此，本文采用无网格格子Boltzmann方法（LBM）数值求解乒乓球三维非定常流场，研究不同温度对乒乓球运动轨迹的影响，分析乒乓球运动轨迹及运动参数随时间变化曲线，并且通过乒乓球运动过程中的瞬态流场细微结构得到影响球体运动诸参数的流动机理。

1 乒乓球分析基础

1.1 乒乓球理论基础

国际乒联规定比赛使用乒乓球质量m=2.7 g，半径为0.02 m，乒乓球桌尺寸分别为:长2.74 m、宽1.525 m、高0.76 m以及网高0.152 5 m[16-17]（见图1）。

根据乒乓球旋转所绕旋转轴不同可分为6大旋转[18]:上旋、下旋、左旋、右旋、顺旋以及逆旋（见图2）。图中，ω为乒乓球旋转角速度，单位是rad/s；v为乒乓球运动速度，单位是m/s。

图1 乒乓球球桌尺寸Figure1 Geometry Parameters of Table

图2 6大旋转Figure2 Six Spin

1.2 受力分析

Magnus效应是普遍存在的，在工程各方面都具有重要的作用[19-23]。在粘性流体中，旋转运动的球体会带动周围流体做旋转运动，当存在绕流时，球体表面将出现压力差，从而产生作用力即Magnus力。飞行中的乒乓球，自身旋转导致的Magnus效应可使其偏离原有轨迹。

在乒乓球比赛过程中，对于上旋球在Magnus效应下，乒乓球产生向下的Magnus作用力，对于下旋球则呈现相反的作用，同样侧向球旋转会导致球体产生侧向偏离[1，24]（见图3）。对于运动过程中的乒乓球，不同旋转方向因Magnus效应导致不同的轨迹偏离（见图4），图中，vx及vy分别为x及y轴方向速度，单位是m/s。

图3 下旋球Magnus效应Figure3 Magnus Effect of Back Spin

图4 不同旋转方向乒乓球运动轨迹Figure4 Trajectory of Table Tennis in Different Spin Direction

乒乓球在运动过程中受到自身重力、空气阻力以及Magnus力等综合作用。对不同旋转方向下乒乓球运动过程进行受力分析（见图5）。图中，Fg、FL以及FD分别为自身重力、Magnus力以及空气阻力。

图5 不同旋转方向乒乓球受力分析Figure5 Force of Table Tennis in Different Spin Direction

2 计算方法[25]

本文采用基于粒子和具有完整拉格朗日函数的方法对格子Boltzmann-BGK方程和不可压缩的Navier-Stokes方程进行求解。此方法无需对流体区域进行网格划分，可较好地体现乒乓球外形特性。对于运动中的乒乓球模拟免去了网格变形的问题，极大地提高了计算的可靠性。计算采用大涡模拟（LES）的方法对乒乓球流场中不同尺度的涡进行处理，对流速梯度较大的区域自动增加粒子数，提高了对乒乓球运动尾涡的捕捉精度。

2.1 格子Boltzmann-BGK方程

格子Boltzmann方法是以介观动理学模型为基础的无网格方法，通过求解描述粒子运动方程从而获得宏观信息，包括:速度、压力以及温度等，方程1为Boltzmann方程:

式中，r为自封闭系统内格子的一个格点；ei为流体第i个粒子的离散速度集合；δt为离散时间步长；t为当前时间步；fi是以速度ei运动的速度分布函数；Si(r,t)为碰撞算子项，表示分子碰撞对速度分布函数的影响。

碰撞算子反映了微观流体粒子的相互作用，对模型能否准确表示处流体系统的物理规律有着关键作用。由于碰撞算子涉及高度非线性积分以及分子间具有耦合作用，因此GROSS等[26]引入方程（2）所示的BGK简化线性化碰撞算子:

式中，τ称为松弛时间是平衡态分布函数。

格子Boltzmann-BGK方程是采用BGK线性碰撞算子简化的Boltzmann方程的一种离散形式，包括空间、速度以及时间的离散[25]。

2.2 基本格子Boltzmann模型

完整的格子Boltzmann模型包含:平衡态分布函数、离散速度模型以及分布函数的演化方程[27]。1992年，QIAN等[28]提出的DdQb系列模型作为格子Boltzmann方法的基本模型，如二维D2Q9[29]及三维D3Q27[30]。本文采用D3Q27模型（见图6）。

图6 D3Q27模型Figure6 D3Q27 Model

方程（3）为D3Q27系列模型的平衡态分布函数，方程（4）为其对应的三维离散速度矢量[30]:

式中，cs为离散系统的声速；u为宏观运动速度；wα为权系数，定义为:

eα为三维离散速度矢量，定义为:

式中:c为离散系统对应的格子速度。

3 参数与建模

图7为球桌及乒乓球的计算域，其大小为:3.5m×1.5m×2.5m。乒乓球在x轴方向（长度方向）的速度vx=6.0 m/s，在y轴方向（高度方向）的速度vy=1.6 m/s，在z轴方向（宽度方向）vz=0；乒乓球旋转速度:ωx=ωy=ωz=0 rad/s，乒乓球初始位置（0，0.2，0.7）m。

图7 计算域Figure7 Computational Domain

本文综合考虑环境温度变化导致空气粘性及密度的改变对乒乓球空气启动作用的影响，采用格子Boltzmann方法，研究不同环境温度对乒乓球轨迹影响。图8为环境温度对空气粘性及密度的影响曲线，表1为3种不同温度时相对应的空气密度粘度[31]。

图8 空气密度、粘性随温度变化Figure8 The Air Density and Viscosity Changes with Temperature

表1 不同环境温度对应的密度及粘度Table1 The Air Density and Viscosity in Different Temperature

4 结果与分析

图9为0、15和30℃3种不同情况下的乒乓球运动轨迹图。图9（a）和（b）分为xoy和xoz平面轨迹图，图9（c）和（d）分别为xoy和xoz平面轨迹局部放大图。

图9 乒乓球轨迹图Figure9 Trajectory of Table Tennis

图9 （a）及图9（c）可见在轨迹上升段，3种温度下乒乓球运动轨迹差别很小；在轨迹下降段，3种不同温度下轨迹出现微小差别；在乒乓球与桌面发生碰撞后，乒乓球轨迹产生较大偏差。乒乓球轨迹随温度升高，其最终点越低，分别为:（2.814 17，0.159 83）、（2.805 82，0.151 67）及（2.820 25，0.147 55），以乒乓球半径为基准，0℃与15℃轨迹x轴方向偏差为41.75%，在y轴方向偏差为40.8%；15℃与30℃轨迹x轴方向偏差为72.15%，在y轴方向偏差为20.6%。

图9（b）及图9（d）中可知在轨迹上升段，0、15以及30 ℃在z轴方向没有明显的偏转；但在轨迹下降段，三者相比产生微小的偏差；在与桌面碰撞后，乒乓球轨迹在3种不同温度情况下产生明显的偏差，并且随着乒乓球继续运动偏转距离越来越大。乒乓球在z轴方向的偏转随着温度升高而降低，轨迹最终点分别为:（0.159 83，0.695 08）、（0.151 67，0.698 93）及（0.147 55，0.699 9），相对于乒乓球半径，0℃与15℃轨迹z轴方向偏差为19.53%；15℃与30℃轨迹z轴方向偏差为4.85%。

图10及图11表明乒乓球运动速度及旋转速度在运动过程中随时间变化，各方向运动速度及旋转速度变化趋势相同，但数值有较大的差别。沿x轴方向，即球体前进方向速度因受到空气阻力作用在运动过程中逐渐减小，与桌面碰撞后速度发生突降；y轴运动速度因重力及Magnus作用力等多重作用力综合作用在运动过程中先减小后增大，在与桌面碰撞后乒乓球速度方向瞬时发生变化，并产生一定的速度损失；z轴运动速度有较小的变化，与桌面碰撞过程三者速度都发生突变。乒乓球旋转速度在x及y轴方向有一定变化，并与桌面碰撞后都发生突变。

虽然环境温度不同乒乓球的运动参数随时间和空间尺度变化的趋势基本相同，但随着时间或飞行距离的变化，尤其是与桌面碰撞反弹后其参数变化有着很大的差异。x轴方向运动速度在轨迹运动过程中3种不同温度条件下无明显差别，碰撞后，速度变化产生明显差别，15℃条件下乒乓球运动速度最小，其他两者无明显区别；在y轴方向运动速度碰撞前无明显差别，碰撞后乒乓球速度随温度升高而增大；在z轴方向运动速度因三者空气粘性阻力不同从而有着巨大的差别，运动速度大小随着温度增大而增大，但在碰撞后因碰撞复杂性，三者速度就很大的差别；旋转速度在时间t=0.15 s之前，在3个方向上不同温度情况下旋转速度无明显差别；旋转速度在时间t=0.15～0.41 s，3个方向有着巨大的差别，并随温度升高而减小；在时间t=0.41 s之后即碰撞后，乒乓球旋转速度发生突变，x、y轴方向随温度增大旋转速度突变量越小，在z轴方向0℃时其旋转速度相对于15℃和30℃条件下突变后的旋转速度较小以及15℃和30℃条件下突变后的旋转速度无明显差别。

乒乓球在空气中运动及与桌面碰撞，其运动轨迹和落点主要受速度和旋转速度的影响，如较之于15℃和30℃，0℃条件下乒乓球运动参数的变化较大，因此其轨迹相比于其他两种情况有很大差别；乒乓球碰撞后的运动参数的突变的差别使得乒乓球在不同温度条件下其轨迹在碰撞之后出现明显的差别。TANAKAK等[32]实验研究球体在不同速比SP（SP=ωr/v）下的升力系数CL（CL=2FL/πρr2v2）及阻力系数CD（CD=2FD/πρr2v2），研究结果表明具有不同的速比SP的乒乓球具有不同的升阻力系数；如图12所示为不同温度下运动过程中的乒乓球速比随时间变化曲线。

由图12可知在同一温度条件下，不同时刻下，乒乓球运动过程中其速比随时间发生变化。对于同一时刻下，乒乓球运动过程中速比随着温度减小而增大；时间t=0.2 s之前，速比差别小，时间t=0.2 s之后，速比差别大。

乒乓球运动中速度和旋转速度的变化，与其在不同时刻的瞬态流场密切相关。宏观上，空气分子的运动因温度升高速度加快，使得密度随温度升高而减小，粘性随温度升高而增大，空气中运动着的乒乓球的升力和阻力与密度和粘性系数成一定关系。因此，必须得到乒乓球运动中周围瞬态流场到微观结构，才能找到流体（空气）影响其各参数变化的流动机理。图13、图14以及图15分别为15℃条件下不同时刻乒乓球流场涡量、速度云图以及乒乓球表面压力图。

图10 乒乓球运动过程中速度随时间变化图Figure10 Velocity in Movement of Table Tennis

图11 乒乓球运动过程中旋转速度随时间变化图Figure11 Spin in Movement of Table Tennis

图12 不同温度下乒乓球速比SP随时间变化Figure12 SPChanges with Time in Different Temperature

图13 乒乓球xoy平面涡量场Figure13 Vorticity Field inxoyFace of Table Tennis

图13 为xoy平面内z轴方向涡量云图，时间t=0.1～0.4 s为碰撞前，t=0.5 s为碰撞后，图中给出的速度v表示该时刻乒乓球瞬时速度方向。由图可见，涡量分布的对称性取决于球体速度的方向，涡量大小决定运动速度。涡量为速度梯度，直接反应球体周围的压力变化。由各时间点的涡量图可看出，原不具有旋转作用的乒乓球由于其周围压力变化从而受到流场粘性阻力的作用产生了旋转作用，同时乒乓球运动速度因周围压力变化受到流场的作用其产生变化，不同温度下具有不同的密度和粘性作用使得乒乓球在运动过程中受到流场作用不同，与图10及图11所表现出球体在运动过程中不同温度下运动速度及旋转速度随时间变化相吻合。

图13表明原不具有旋转的乒乓球由于球体在运动过程中受到流场作粘性阻力作用产生一定的旋转，导致其上下部分周向速度与球体运动速度方向不同，造成乒乓球速度流场的不对称性，因此，如图14时间t=0.1～0.4 s为碰撞前，t=0.5 s为碰撞后，乒乓球因流场不对称性（关于瞬时速度方向）从而受到流场作用力，不同温度具有不同的粘性作用，造成这与图9所表现出不同温度下球体宏观轨迹发生变化相吻合。

图15为不同时刻下乒乓球表面静态压力分布，对于初始没有旋转的乒乓球在运动过程中其上、下侧存在压力差，且在乒乓球上、下侧表面分布不均；乒乓球前、后、上以及下侧的表面压力分布随着时间的变化产生巨大的变化，表明乒乓球在运动过程中其受到的空气阻力及升力发生变化的，并与图12及相关实验研究相吻合。

不同温度具有不同的密度以及粘性对乒乓球的作用能在轨迹特别是碰撞后的轨迹中明显表现出来，20℃和30℃两者在碰撞后z轴方向旋转速度突变量相当，但因密度便显出很大差别进而影响其碰撞后轨迹。在空气动力粘度方面粘性随温度升高而增大，并能在旋转速度中表现出来:在x轴方向温度小的即动力粘度小的总是先出现变化的，这也就是0℃运动参数特别其旋转速度变化相对于15℃和30℃较为提前，而后两者空气粘度相对差别不大也就是后两者变化具有相似性。

图14 乒乓球xoy平面速度场Figure14 Velocity Field inxoyFace of Table Tennis

图15 不同时刻下乒乓球表面压力图（xoy平面）Figure15 Pressure in Face of Table Tennis at Different Time（xoy Face）

5 结论与建议

5.1 结论

采用格子Boltzmann方法的数值算法，不同环境温度对乒乓球运动轨迹的影响，获取乒乓球运动轨迹、速度及旋转速度随时间的变化规律，并分析乒乓球运动过程中速度场、涡量场以及乒乓球表面压力分布，得到如下结论。

（1）环境温度影响空气粘性及密度；乒乓球在运动过程中因空气粘性作用其旋转速度及运动速度产生巨大的变化，进而影响运动过程中乒乓球气动力系数的变化，由乒乓球运动过程中受力分析可知气动力系数的变化会导致乒乓球在运动过程中受力是发生变化的。因此，运动学方法不能过于简化运动过程的受力作用，需考虑空气粘性作用。

（2）不同温度对乒乓球轨迹存在很大的影响。以乒乓球半径作为参考，在长度方向上，0℃与15℃偏差为41.75%，15℃与30℃偏差为72.15%；在高度方向，0℃与15℃偏差为40.8%，15℃与30℃偏差为20.6%；在宽度方向，0℃与15℃偏差为19.53%，15℃与30℃偏差为4.85%.

（3）不同温度对乒乓球运动参数具有很大的影响。0℃条件下其三个方向的运动速度和旋转速度与15℃和30℃有着很大的差别，而15℃与30℃在个别的方向上或与桌面发生碰撞才产生较大差别。

（4）乒乓球的流场结构及其表面压力分布随时间变化，对乒乓球的作用效果与乒乓球运动参数变化及轨迹偏转相吻合。研究乒乓球运动过程中流场细微结构的变化探究乒乓球流动机理，为乒乓球运动提供一定的理论指导。

5.2 建议

基于以上研究，提出以下建议:（1）乒乓球与桌面碰撞时间短，但碰撞后运动参数运动轨迹都有很大的变化，故球员要具有很敏捷的判断及击球反应；（2）就本文乒乓球运动轨迹而言，当环境温度较高，球员应当对于击球点做出稍向下移动且稍向右移动的微小动作；当环境温度较低时，球员应当对击球点做出稍向上移动且一定的右移动的较小动作。

[1]余万，李春，朱玲，等.不同比赛气候条件对乒乓球飞行轨迹影响研究[J].运动，2016，138（10）:23-24.

[2]窦远行.中国乒乓球队历届奥运战绩:里约第5次包揽金牌[EB/OL].（2016-08-18）[2016-12-17]http://news.sohu.com/20160818/n464841 590.shtml.

[3]曹犇，雷正方.中国乒乓球赛事产业化发展前景规划[J].西安体育学院学报，2017，34（3）:295-299.

[4]余万，李春，任杰，等.基于计算流体力学方法的乒乓球轨迹仿真[J].上海体育学院学报，2017，41（3）:89-94.

[5]唐建军.乒乓球技术学习:理论解释与实际运用[J].北京体育大学学报，2006，29（2）:259-261.

[6]虞丽娟，张辉，凌培亮.乒乓球比赛技战术分析的系统研究与应用[J].上海体育学院学报，2008，32（6）:39-43.

[7]魏利婕，史桂兰.乒乓球削球打法的现状与持续发展的可行性研究[J].北京体育大学学报，2006，29（12）:1709-1711.

[8]宋亚炳，王伟强，刘建进，等.竞赛教学模式在乒乓球选项教学中的实验研究[J].体育学刊，2014，21（2）:90-94.

[9]杨华，关志明.基于ODE的乒乓球运动轨迹仿真研究[J].计算机仿真，2011，28（9）:230-233.

[10]章逸丰，熊蓉.乒乓球机器人的视觉伺服系统[J].中国科学:信息科学，2012，42（9）:1115-1129.

[11]王奇志，杨晓晓.乒乓球轨迹预测的研究与仿真[J].计算机工程与科学，2013，35（2）:164-169.

[12]NAKASHIMA A，KOBAYASHI Y，OGAMA Y，et al.Modeling of Re⁃bound Phenomenon between Ball and Racket Rubber with Spinning Ef⁃fect in ICCAS-SICE 2009-ICROS-SICE International Joint Confer⁃ence[C].Japan:Fukuoka，2009:2295-2300

[13]MANIN L，POGGI M，HAVARD N.Vibrations of table tennis racket composite wood blades:modeling and experiments[J].Procedia Engi⁃neering，2012，34（2）:694-699.

[14]KUI OU，Patrice Castonguay，Antony Jameson.Computational Sports Aerodynamics of a Moving Sphere:Simulating a Ping Pong Ball in Free Flight in AIAA Applied Aerodynamics Conference[C].Hawaii:HEHTAulu，2011:1-16.

[15]MEHTA R D.Aerodynamics of sports balls[J].Annual Review of Fluid Mechanics，1985，17（1）:151-189.

[16]THE INTERNATIONAL TABLE TENNIS FEDERATION.ITTF Hand⁃book 2011/2012[EB/OL].（2011-12-25）[2016-12-17]http://www.ittf.com/ittf_handbook/ittf_hb.html.

[17]THE INTERNATIONAL TABLE TENNIS FEDERATION.ITTF Tech⁃nical Leaflet T1:The Table[EB/OL].（2013-05-12）[2016-12-17]http://ittf.com/stories/pictures/T1_The_Table_BoD2013.pdf.

[18]王正伦.教你打乒乓球[M].南京:江苏科学技术出版社，2012.

[19]MEININGER R C.ICAR 101[J].Rock Products，2004，（10）:16-16.

[20]BUTLER AMES.Whirling Spools Lift This Plane[N/OL].https://books.google.com/books?id=xSgDAAAAMBAJ&pg=PA26&dq=Popular+Sci⁃ence+1931+plane#v=onepage&q&f=false，2016-12-17/1930-09-17

[21]G0TTSCJALK M J.Self-steering axle suspension system having a rota⁃ry stabilizer:US，US7360773[P]:2008-05-13.http://www.google.com/patents/US7360773

[22]龚志斌，李杰，张辉.旋转圆柱对翼型气动性能特性影响的数值模拟研究[J].空气动力学学报，2015，33（2）:254-258.

[23]陈东阳，LAITH K.ABBS，芮筱亭，等.结构误差对旋转稳定弹丸气动特性影响的数值模拟[J].空气动力学学报，2014，32（5）:705-710.

[24]SARAFIAN H.Impact of the Drag Force and the Magnus Effect on the Trajectory of a Baseball[J].World Journal of Mechanics，2015，05（4）:49-58.

[25]郭照立.格子Boltzmann方法的原理及应用[M].北京:科学出版社，2009.

[26]BHATNAGAR P L，GROSS E P，KROOK M.A Model for Collision Processes in Gases.I.Small Amplitude Processes in Charged and Neu⁃tral One-Component Systems[J].Physical Review，1954，94（3）:511-525.

[27]何雅玲，李庆，王勇，等.格子Boltzmann方法的工程热物理应用[J].科学通报，2009，（18）:2638-2656.

[28]QIAN Y H，D'HUMIERES D，LALLEMAND P.Lattice BGK Models for Navier-Stokes Equation[J].Epl，1992，17（6BIS）:479.

[29]冉政，陈建.关于格子Boltzmann方法的几个问题（I）[J].空气动力学学报，2016，34（3）:334-339.

[30]张良奇.格子Boltzmann方法基础理论研究及其在不可压缩流动中的应用[D].重庆:重庆大学，2014:15-84.

[31]杨世铭，陶文铨.传热学[M].北京:高等教育出版社，2006.

[32]TANAKA K，FUKUIU T，Miyazaki T，et al.Aerodynamic Properties of a Table Tennis Ball[J].Journal of Japan Society of Fluid Mechanics，2014，33（1）:37-45.