1、1数学建模：二轮滑板的最佳运动方式摘要“游龙板”是近年来我国较为流行的一种新型滑板。该项运动可以提高人身体的柔软性和平衡感，是一项深受青少年喜爱的体育运动。与普通滑板相比，操作时人不用下板，利用双脚前后摆动就可以驱动游龙板前进。许多人感觉游龙板的前进似乎不需借助于外力，运动原理难以想象。游龙板结构虽然简单，但设计得却相当巧妙。它有两个结构特殊的轮子，可以将板的横向摆动转换成向前的推动力。滑板在运动过程中主要受到地面给轮子的静摩擦力和人对板的蹬力，其中人对板的蹬力使活力板转动，使轮架偏离一定角度，同时地面对轮子的静摩擦力使得滑板向前运动，而人体的扭动角度对板的转动起着至关重要的作用。 平地滑板前
2、进速度与水平推进力相关，在一个摆动周期内，水平推进力越大，活力板获得的加速度最大，即在一个周期内速度的最大值尽可能的大。而水平推进力和轮架旋转的角度 有关，建立 关于 的模型，解得当 为水 平F时，地面对滑板的摩擦力完全提供水平推进力，此时 最大。人对滑板90 水 平F的作用力做功转化为板子的转动动能和轮子的转动动能，由刚体转动的动能定理建立 ， ， ， 的关系，假设 ，化简得到 和 的关系，板轮 wtsin0板t进而建立微分方程，通过龙格库塔法解得四分之一周期内满足微分方程的 ，板再对 积分得到板子最佳旋转角度为 37 ，即最佳扭动角度。板 关键词： 游龙板 扭动角度 龙格库塔法 2一、问题
3、重述“二轮滑板”又称作为“两轮滑板” 、 “陆地冲浪板”是一种新型的滑板，和普通的四轮的滑板不一样，它只有两个轮。因其运动起来像龙和蛇一样扭动，所以国内多大称为“游龙板”或“活力板” 。它以人体运动理论和巧妙的力学原理，主要利用身体（腰部及臀部） 、双脚扭动及手的摆动来驱动活力板前进。使滑板运动得到淋漓尽致的发挥。现请建立人身体扭动角度与平地滑板前进速度的关系模型，并找到最佳扭动角度。二、问题分析2.1 活力板的构造活力板有前后两个面板，中间用轴连接，两个面板可以绕中间的连接轴相对转动面板下方安装有轮架和车轮，轮架可以绕倾斜的转轴自由转动这个构造的巧妙之处在于，轮架转轴的延长线并不通过车轮与地
4、面的接触点，而是有一小段垂直距离如图 1 所示。2.2 活力板的受力当人对板施加力 时，车板会有一个下角度的倾斜，地面对车轮的静摩脚F擦力使轮架绕转轴旋转，这时静摩擦力会产生水平当向上的分力推动活力板向前运动。人通过用身体左右晃动来使活力板运动，板面是在水平方向上左右摆动的，所以可以将人体的扭动角度近似的看成是板身摆动的角度，板身摆动的角度与轮轴摆动的角度存在一定的关系，轮轴摆动角度与水平方向上的分力有关，根据水平分力算出滑板速度，最终可以找到速度与扭动角度的关系。三、模型假设（1）假设人体扭动角度就是板的转动角度；（2）假设轮架的转动随时间成正弦函数关系；（3）假设人施加给滑板的力是均匀大小
5、的力；（4）假设运动过程不受空气阻力影响；（5）假设支持力不随人的扭动角度变化而变化。3四、符号说明符号 说明v 滑板前进速度F 水平 滑板受到水平方向的分力F 脚 使滑板作用摆动的力转轴与竖直方向的夹角轮架的旋转角度滑板转过的角度板板的旋转角速度轮滑轮绕轮轴旋转的角速度L脚到滑板中心的距离L 轮轴到滑轮的距离R 滑轮半径N 地面对滑轮与人整体支持力板m板的质量轮轮的质量五、模型的建立与求解5.1 模型的建立（一）轮子的结构由于游龙板前后两轮结构相同，只分析其一即可。轮子是通过轮架和一倾斜的转轴安装在板 E，轮子和轮架可以绕转轴的轴线 自由转动。为分析方p便，将轮子结构简化如图 1 所示。4运
6、动方向地面板面轮子转轴轮架pqAONDCqpB图 1 轮子结构简化侧视图图中用线段 AB 和 BO 代表转轴和轮架，O 点为轮子的圆心， D 点为轮子的触地点，C 点为 AB 的延长线与地面的交点。轮子的特殊结构体现在两个地方：第一，转轴 AB 是倾斜的，设与竖直方向的夹角为 ；第二， C 点与 D 点是不重合的，而是落在 D 点的前面。当人站在板上时，地面对后轮 D 点的支持力设为 N。人的双脚在板面上沿与前进方向相垂直的方向摆动，轮子和轮架就会绕转轴 摆动，进而产生向前的推动力。下面对推动力产生的原因进行具体分p析。（二）推动力产生原因分析如图 2 所示，设轮架绕转轴 摆过任一角度 ，轮子
7、摆过一定角度后，pC、D 蹲点的位置也发生了变化，设为 和 。轮子受地面支持力 N 的方向竖CD直向上， 点与地面问静摩擦力设为 F，与水平前进方向夹角为 ，做直 90线 ， 与滑板前进方向一致。TpS(1)支持力 N 对转轴 的转矩产生转矩的力的计算轮子支持力 N 产生使轮子绕转轴 回摆的转矩。N 在垂直于 方向的p TD分力 (也垂直于 )是使轮子绕 转动的有效力，其大小为： 。p sinN5运动方向 Aqpqp脚FBFT NNCS D90O图 2 轮子绕转轴转过任意角度 B 的立体图pTNNPSHDC90B图 3 对 转动力矩的计算Np产生转矩的力臂的计算将图 2 对应局部放大如图 3
8、所示。做 垂直于滑板前进方向的水平线HD于 H 点，则 面 B 。CDC由于两平面 和 B 平行，所以 平面 B ，那么，分力NTNPC所在直线到转轴 的距离就等于 点到平面 B 的距离，即 的长度。N p HD由图 3 可知：= (1)HDsinC6即为 对转轴 的力臂。Np则支持力 N 对转轴 的转矩 L 为：(2)sinsisinDCNH(2)静摩擦力 F 对转轴 的转矩p产生转矩的力的计算如图 4 所示，静摩擦力 F 产生阻止轮子回摆的力矩， F 在垂直于 方向TD的分力 (也垂直于 )为产生转矩的有效力。pNTpGPSFDCF 90R图 4 对 转动力矩的计算Fp做辅助线 于 S 点
9、，TR 于 R 点，连接 SR，可证明：SRTDD，从而可以求出 F 与 的夹角 的大小，设其大小为 ，可求出：FD T (3) sin)90cos(in)90cos(cos TR所以：(4)22si1cs1sin(5)niF产生转矩的力臂的计算如图 4 所示，转轴 平行于平面 ，所以， 对 转动力臂为 平pDTFpp行于平面 的距离。过 点做平面 的垂线，交平面 于点 G，连FDTCDT接 G， 两点， 即为所求的力臂。7可以证明： = SRT，设其大小为 ，由图可得出：GDC (6)cosincosinsin DTRDSTRSTtg所以：(7)22cosincossi(8)22ii DCC
10、G于是可求出力 F 阻止轮架回摆的力矩 ：L(9)2cos1intgFL(3)推动力的大小令(1)(2)中的 L= 。可求出：(10)2cos1intgNF在轮子绕转轴摆过任意角度 时，滑板受到的推动力实际上是 F 在水平前进方向的分力。(11)2costg1in-isinNF水 平综上所述，使游龙板前进的推动力来源于轮子与地面间的静摩擦力。轮子的结构决定了只要轮子绕转轴转过一个角度 ，轮子与地面间就会有静摩擦力F 存在， F 沿运动方向的水平分量驱动游龙板前进。所以在光滑的路面上，静摩擦力不存在，游龙板无法前进。在做该项运动时，操作者的双脚在板面上要不断地用力前后摆动，目的就8是为了保持角
11、一直存在，即角 必须依靠外力才能保持。如图 2 所示。此时人脚施力方向为垂直纸面向外，大小设为 ，则 应脚F脚该与力 F 在垂直于运动方向的分力相平衡，即有：(12)cosF脚可见，力 是力 F 赖以存在的保证，只要 存在， F 就存在，所以，系脚 脚统所受的推动力最终还是来源于人。尽管脚的施力方向与滑板前进的方向是垂直的，但是由于轮子结构的特殊性，可以将其转换成向前的推动力，使滑板前进。系统的动能还是来源于人的做功。分析速度与水平方向的力的关系可知，在一个周期内当水平方向的力达到最大时，其加速度最大，当力为零时，加速度为零，速度最大，假设 与时间t 的关系呈正弦规律变化幅值为 ，角频率为 ，
12、如下式：0(13)twsin0另外，轮子的结构使得 变化时，也会牵扯到系统重心的高低变化。时系统重心最低，处于稳定平衡状态； 时系统重心最高，处于0 18不稳定平衡状态。所以，游龙板只能朝一个方向运动，不能倒滑，因此 取0.9由（13）式水平方向的力与 的关系，将 代入(11)式。在 matlab 中6画得图形如图 5 所示，由图求得最优解为当 时，水平力最大。本文将人90身体扭动角度近似的看做滑板平面转动的角度 。0 0.5 1 1.5 2 2.5 3 3.500.10.20.30.40.50.60.7beltaforce9图 5 轮子旋转角度与水平推进力的关系人对板施加力使得轮轴控制轮架与
13、轮子进行滑动，首先，为了找到 与的关系，建立了刚体转动的动能定理（忽略转动过程中阻力做功） ，即人对滑板做的功等于滑板动能的增量与轮子动能增量的和，如（14）式所示，其中是关于时间 t 的函数：脚轮板 ， F(14)dLFJ0221脚轮轮板板因为 与 的关系如（12）式，将 代入上式解得此时的 ，脚 9脚F与 可以通过转动惯量计算公式解出：板J轮(15)21LmJ板板 (16)(rl轮轮将（12） （15） （16）式代入（14）式，两边同时对 t 求导可得： (17)板板轮脚板板 213101cosin)(2dt LmtRLF最终建立扭动角度与平地滑板速度的模型如：(18)t0223020d
14、sincosg1i-si1tcostin)(m2dt costgin-1isi dtg-板脚 板板轮脚板板水 平 总 总水 平wNFLRLNFtavt将数据代入微分方程，使用 matlab 中的函数 ode45，通过龙格库塔法，求得 关于时间 的变化关系，如图 6 所示。因为 时扭动角度最佳，将其板t 9010代入（13）式得到此时的 t=0.25，代入（18）式即得到当 时扭动角度最37佳。使用 matlab 中的函数 ode45，通过龙格库塔法，求得 关于时间 的变化关系板t如图 6 所示：0 0.01 0.02 0.03 0.04 0.05 0.06 0.0705101520253035
15、404550t/somega/rad*s-1图 6 关于时间 的变化关系板六、模型的评价与推广本文将模糊的最佳扭动角度的求解问题转化为滑板旋转角度的问题，这样就可以很清楚地求出滑板前进速度与滑板旋转角度之间的函数关系，将复杂的问题通过转化变为简单的问题；对于刚体的转动过程只考虑初始状态和末状态，不考虑过程中作用力的具体情况，是问题的求解变得简单；在求解微分方程是运用了龙格库塔法，计算过程中可以改变步长，不需要计算高阶导数具有精度高、收敛、稳定等优点。在列刚体转动的动能关系时没有考虑摩擦力做功；在模型求解时很多参数都是估计得来的，存在一定的误差。七、参考文献1 刘卫国. MATLAB 程序设计与
16、应用 M. 北京：高等教育出版社，2006.112 卓金武. MATLAB 在数学建模中的应用M. 北京：北京航空航天大学出版社，2011.3 活力板、游龙板视频教学.http:/ 2014.8.194 游龙板教学 . http:/ 2014.8.195 张晓娜，梁亚波. 活力板的力学原理 J.用科学的眼光看世界,2012(7):46-47.1八、附录图 1 程序 1：%main.cclc;clear;close all;global j_boad j_wheel lp alpha N belta0 jiaopinlvm_wheel=0.4;m_boad=3;lp=56*0.01;l=2*0.
17、01;r=4*0.01;alpha=pi/6;N=700;belta_row=best_belta(alpha);belta0=max(belta_row);j_boad=1/12*m_boad*lp2;j_wheel=m_wheel*(2*l*r+r2);jiaopinlv=8*pi;T=pi/2/jiaopinlv;%t,omega=ode45(myfun,0,T,0.1);%angle=para_integar(0,T,omega,length(t)*180/pi;belta=belta0*sin(jiaopinlv*t);figure;plot(t,omega);xlabel(t/s)
18、;2ylabel(omega/rad*s-1);grid onfigure;plot(t,belta);grid on程序 2：function b,f = best_belta( alpha )%UNTITLED2 请在此处输入函数概要% 请在此处输入详细说明belta=linspace(0,pi,1000);tmp1=sin(alpha)*sin(belta).*sin(belta);tmp2=1-sin(alpha)*sin(alpha)*sin(belta).*sin(belta);tmp3=1+tan(alpha)*tan(alpha)*cos(belta).*cos(belta);
19、force=tmp1./(sqrt(tmp2./tmp3);f=max(force);n=find(force=f);b=belta(n);figure;plot(belta,force);xlabel(belta);ylabel(force);grid onend程序 3：function dy = myfun( t,omega )%UNTITLED2 请在此处输入函数概要% 请在此处输入详细说明% global temp1 temp2 temp3global j_boad j_wheel lp alpha N belta0 jiaopinlvbelta=belta0*sin(16*pi*t
20、);%tmp1=N*sin(alpha)*sin(belta)*cos(belta);tmp2=1-sin(alpha)*sin(alpha)*sin(belta)*sin(belta);tmp3=1+tan(alpha)*tan(alpha)*cos(belta)*cos(belta);force=tmp1/(sqrt(tmp2/tmp3);temp1=lp*force/j_boad;temp2=2*j_wheel*(pi/2*jiaopinlv*cos(jiaopinlv*t)*(-pi/2*jiaopinlv2*sin(jiaopinlv*t);%dy=temp1-temp2/(j_boad*omega);3end程序 4：function s = para_integar(a,b,y,n )%对离散点的抛物线积分数值解法%a 积分下限值%b 积分上限值%n 分割区间数h=(b-a)/n;temp1=0;temp2=0;for i=1:n-1if rem(i,2)=0temp1=temp1+y(i+1);endif rem(i,2)=0temp2=temp2+y(i+1);endendtemp=y(1)+y(n)+4*temp1+2*temp2;s=h/6*temp;end
