篮球运动数值模拟程序

1、气固两相流动与数值模拟篮球运动过程数值模拟程序辛传贤（东南大学化学化工学院，江苏 南京 211189）附录一functionSxyz,Vxyz,Wxyz,Impact_I1,Impact_I2,Impact_I12,Impact_I3,Impact_I4,warning=Basketball_Move(vf,vector_vf,vl,vector_vl,w,vector_w,s)%求解篮球运动的各种信息(主要为篮球的运动轨迹)%vf风速的大小,vector_vf风速的方向向量;vl篮球的初始速度大小,vector_vl篮球初始速度的方向向量;w篮球的初始旋转角速度大小,vector_w篮球的初

2、始角速度方向向量;s篮球的初始位置坐标%求解篮球运动的各种信息,Sxyz为轨迹信息,Vxyz为速度信息,Wxyz为角速度信息,Impact_I1、Impact_I2和Impact_I3分别为篮球与篮圈、篮板和地面碰撞的信息,warning为预警信息pg=1.185;%空气密度(298K)(国际单位制，下同)u=1.84*10-5;%空气黏度(298K)g=9.8;%重力加速度m=0.6;%篮球质量dp=0.246;%篮球直径rp=dp/2;%篮球半径Vb=4*pi/3*rp3;%篮球体积dl=0.45;%篮圈内径(直径)dt=0.02;%圈条直径A=(dl+dt)/2;%xz平面上圈条的环心线

3、的半径a1=dt/2;%圈条半径a2=a1+rp;%篮球与篮圈相撞的半径bx=7.5;%篮圈的圆心坐标by=3.05;bz=1.575;h1=2.75;%篮板下上界高度坐标(y轴),(注意x、y、z的方向笛卡尔坐标系)x向右，y向上，z向前),(下同)h2=3.8;t1=6.6;%篮板左右界宽度坐标(x轴)t2=8.4;z1=1.2;%篮板前位置坐标(z轴)rc=0.1;%篮球架立柱的半径vector_vl=vector_vl/norm(vector_vl);%投篮的单位方向矢量vx(1)=vl*vector_vl(1);%篮球初始三维速度vy(1)=vl*vector_vl(2);vz(1)

4、=vl*vector_vl(3);vector_vf=vector_vf/norm(vector_vf);%风场的单位方向矢量vfx=vf*vector_vf(1);%风场的三维速度vfy=vf*vector_vf(2);vfz=vf*vector_vf(3);vector_w=vector_w/norm(vector_w);%篮球的初始旋转方向矢量wx(1)=w*vector_w(1);%篮球的初始旋转三维分量wy(1)=w*vector_w(2);wz(1)=w*vector_w(3);sx(1)=s(1);%篮球初始球心位置坐标(9.3,2.0,5.8)，罚球线最外端(x轴正向)sy(1

5、)=s(2);sz(1)=s(3);I1=1;%分别记录篮球与篮圈、篮板、地面和篮球架立柱撞击信息I2=1;I12=1;%(篮圈同时与篮圈和篮板相撞)I3=1;I4=1;Impact_I1(1,1)=0;%分别记录篮球与篮圈、篮板、地面和篮球架立柱撞击次数信息,(注意第一个信息指标无意义)Impact_I2(1,1)=0;Impact_I12(1,1)=0;%(篮圈同时与篮圈和篮板相撞)Impact_I3(1,1)=0;Impact_I4(1,1)=0;Impact_I1(2,1)=0;%分别记录篮球与篮圈、篮板、地面和篮球架立柱撞击前经历的时间信息,(注意第一个信息指标无意义)Impact_

6、I2(2,1)=0;Impact_I12(2,1)=0;Impact_I3(2,1)=0;Impact_I4(2,1)=0;Impact_I1(3,1)=0;%分别记录篮球与篮圈、篮板、地面和篮球架立柱撞击过程经历的时间信息,(注意第一个信息指标无意义)Impact_I2(3,1)=0;Impact_I12(3,1)=0;Impact_I3(3,1)=0;Impact_I4(3,1)=0;earlywarning=rp/10;%预警距离为篮球半径的长度warning=0;%预警次数记录器normal_t=10-5;for i=2:100000%循环次数 earlywarning1=(sx(i-

7、1)-bx)2+(sy(i-1)-by)2+(sz(i-1)-bz)2-(A+earlywarning)2-(a2+earlywarning)2)/(2*(A+earlywarning)2+(sy(i-1)-by)2-(a2+earlywarning)2<=0);%判断篮球是否与篮圈达到预警距离 earlywarning2=(sx(i-1)>=(t1-rp-earlywarning)&sx(i-1)<=(t2+rp+earlywarning)&sy(i-1)>=(h1-rp-earlywarning)&sy(i-1)<=(h2+rp+ear

8、lywarning)&sz(i-1)>=(z1-rp-earlywarning)&sz(i-1)<=(z1+rp+earlywarning);%判断篮球是否与篮板达到预警距离earlywarning3=(sy(i-1)<=rp+earlywarning);%判断篮球是否与地面达到预警距离earlywarning4=(sx(i-1)-bx)2+sz(i-1)2<=(rc+rp+earlywarning)2);%判断篮球是否与篮球架立柱达到预警距离if earlywarning1|earlywarning2|earlywarning3|earlywarnin

9、g4%篮球在预警距离内时,迭代时间步长t为10-4秒;在预警距离外时,迭代时间步长t为10-3秒warning=warning+1;%预警次数记录器warning_t=10-5;t=10-5; elsenormal_t=0.3*10-4;t=0.3*10-4; endImpact_1=(sx(i-1)-bx)2+(sy(i-1)-by)2+(sz(i-1)-bz)2-A2-a22)/(2*A)2+(sy(i-1)-by)2-a22<=0);%判断篮球是否与篮圈相撞 Impact_2=(sx(i-1)>=(t1-rp)&sx(i-1)<=(t2+rp)&sy(i

10、-1)>=(h1-rp)&sy(i-1)<=(h2+rp)&sz(i-1)>=(z1-rp)&sz(i-1)<=(z1+rp);%判断篮球是否与篮板碰撞 h1,h2 篮板下上界高度坐标(y轴); t1,t2 篮板左右界宽度坐标(x轴); z1 篮板前位置坐标(z轴)(注意x、y、z的方向笛卡尔坐标系x向右，y向上，z向前)Impact_3=(sy(i-1)<=rp);%判断篮球是否与地面碰撞Impact_4=(sx(i-1)-bx)2+sz(i-1)2<=(rc+rp)2);%判断篮球是否与篮球架立柱碰撞if Impact_1=1&a

11、mp;Impact_2=0 %判断篮球是否仅与篮圈相撞Impact=1;%篮球与篮圈相撞的记录器I1=I1+1;%篮球与篮圈的撞击信息Impact_I1(1,I1)=i-1;%记录篮球与篮圈碰撞的次数信息Impact_I1(2,I1)=(i-1-warning)*normal_t+warning*warning_t;%记录篮球与篮圈碰撞前经历的时间信息if (sx(i-1)-bx)2+(sz(i-1)-bz)2)-A2>=0 %判断篮球是否在篮圈中心线的外侧jda=asin(sy(i-1)-by)/a2);%y方向的偏角(篮圈外侧)else jda=pi-asin(abs(sy(i-1)

12、-by)/a2);%y方向的偏角(篮圈内侧)(注意abs)endsinjdb=(sz(i-1)-bz)/(A+a2*cos(jda);%xz方向的偏角正弦值cosjdb=(sx(i-1)-bx)/(A+a2*cos(jda);%xz方向的偏角余弦值x0=(A+a1*cos(jda)*cosjdb+bx;%篮球与篮圈的接触点坐标(注意:此值为下述方eqa1,eqa2,eqa3的迭代初值y0=a1*(sy(i-1)-by)/a2+by;z0=(A+a1*cos(jda)*sinjdb+bz;ssxi=sx(i-1);%对篮球位置进行复制ssyi=sy(i-1);sszi=sz(i-1);ss=A,

13、a1,bx,by,bz,ssxi,ssyi,sszi;%Function求解碰撞公切点的系数矩阵xyz=fsolve(xx)Function(xx,ss),x0,y0,z0);%调用函数function F=Function(xx,ss)%xx为篮球与篮圈的碰撞公切点的迭代初始值(x0,y0,z0)，ss为Function求解碰撞公切点的系数矩阵x=xyz(1);y=xyz(2);z=xyz(3);%篮球与篮圈的公切点的坐标(x,y,z)afx=(x-bx)2+(y-by)2+(z-bz)2-A2-a12)/A2*(x-bx);%篮球与篮圈接触点处的法向量afy=(x-bx)2+(y-by)2

14、+(z-bz)2-A2-a12)/A2*(y-by)+2*(y-by);afz=(x-bx)2+(y-by)2+(z-bz)2-A2-a12)/A2*(z-bz);Function_A=afx,afy,afz;%Function_A为撞击点的法向量Function_vv=vx(i-1),vy(i-1),vz(i-1);%Function_vv碰撞时篮球的速度Function_vf=vfx,vfy,vfz;%Function_vf碰撞时风速Function_ww=wx(i-1),wy(i-1),wz(i-1);%Function_ww碰撞时篮球的角速度Function_s=sx(i-1),sy(

15、i-1),sz(i-1); %篮球撞击点的球心坐标V,W,S,rotation_t=Impact_Rotation_1(Function_A,Function_vv,Function_vf,Function_ww,Function_s,Impact);%求解篮球碰撞过程中的旋转角速度及反弹速度(完全非弹性碰撞过程)vx(i-1)=V(1);%碰撞后篮球的速度vy(i-1)=V(2);vz(i-1)=V(3);wx(i-1)=W(1);%碰撞后篮球的角速度wy(i-1)=W(2);wz(i-1)=W(3);sx(i-1)=S(1);%碰撞后篮球球心坐标sy(i-1)=S(2);sz(i-1)=S

16、(3);Impact_I1(3,I1)=rotation_t;%记录篮球与篮圈撞击过程经历的时间信息endif Impact_2=1&Impact_1=0%判断篮球是否仅与篮板碰撞Impact=2;%篮球与篮板相撞的记录器I2=I2+1;%篮球与篮板的撞击信息Impact_I2(1,I2)=i-1;%记录篮球与篮板碰撞的次数信息Impact_I2(2,I2)=(i-1-warning)*normal_t+warning*warning_t;%记录篮球与篮圈碰撞前经历的时间信息Function_A=0,0,1;%Function_A为撞击点的法向量Function_vv=vx(i-1),

17、vy(i-1),vz(i-1);%Function_vv碰撞时篮球的速度Function_vf=vfx,vfy,vfz;%Function_vf碰撞时风速Function_ww=wx(i-1),wy(i-1),wz(i-1);%Function_ww碰撞时篮球的角速度Function_s=sx(i-1),sy(i-1),sz(i-1); %篮球撞击点的球心坐标V,W,S,rotation_t=Impact_Rotation_1(Function_A,Function_vv,Function_vf,Function_ww,Function_s,Impact);%求解篮球碰撞过程中的旋转角速度及反

18、弹速度(完全非弹性碰撞过程)vx(i-1)=V(1);%碰撞后篮球的速度vy(i-1)=V(2);vz(i-1)=V(3);wx(i-1)=W(1);%碰撞后篮球的角速度wy(i-1)=W(2);wz(i-1)=W(3);sx(i-1)=S(1);%碰撞后篮球球心坐标sy(i-1)=S(2);sz(i-1)=S(3);Impact_I2(3,I2)=rotation_t;%记录篮球与篮圈撞击过程经历的时间信息endif Impact_2=1&Impact_1=1%判断篮球是否同时与篮圈与篮板相碰撞I12=I12+1;%篮球同时与篮圈和篮板撞击的信息Impact_I12(1,I12)=i

19、-1;%记录篮球同时与篮圈和篮板碰撞的次数信息Impact_I12(2,I12)=(i-1-warning)*normal_t+warning*warning_t;%记录篮球同时与篮圈和篮板碰撞前经历的时间信息if (sx(i-1)-bx)2+(sz(i-1)-bz)2)-A2>=0 %判断篮球是否在篮圈中心线的外侧jda=asin(sy(i-1)-by)/a2);%y方向的偏角(篮圈外侧)else jda=pi-asin(abs(sy(i-1)-by)/a2);%y方向的偏角(篮圈内侧)(注意abs)endsinjdb=(sz(i-1)-bz)/(A+a2*cos(jda);%xz方向

20、的偏角正弦值cosjdb=(sx(i-1)-bx)/(A+a2*cos(jda);%xz方向的偏角余弦值x0=(A+a1*cos(jda)*cosjdb+bx;%篮球与篮圈的接触点坐标(注意:此值为下述方eqa1,eqa2,eqa3的迭代初值y0=a1*(sy(i-1)-by)/a2+by;z0=(A+a1*cos(jda)*sinjdb+bz;ssxi=sx(i-1);%对篮球位置进行复制ssyi=sy(i-1);sszi=sz(i-1);ss=A,a1,bx,by,bz,ssxi,ssyi,sszi;%Function求解碰撞公切点的系数矩阵xyz=fsolve(xx)Function(x

21、x,ss),x0,y0,z0);%调用函数function F=Function(xx,ss)%xx为篮球与篮圈的碰撞公切点的迭代初始值(x0,y0,z0)，ss为Function求解碰撞公切点的系数矩阵x=xyz(1);y=xyz(2);z=xyz(3);%篮球与篮圈的公切点的坐标(x,y,z)afx=(x-bx)2+(y-by)2+(z-bz)2-A2-a12)/A2*(x-bx);%篮球与篮圈接触点处的法向量afy=(x-bx)2+(y-by)2+(z-bz)2-A2-a12)/A2*(y-by)+2*(y-by);afz=(x-bx)2+(y-by)2+(z-bz)2-A2-a12)/

22、A2*(z-bz);Function_A=afx,afy,afz;%Function_A为撞击点的法向量Function_AA=0,0,1;%Function_A为撞击点的法向量Function_vv=vx(i-1),vy(i-1),vz(i-1);%Function_vv碰撞时篮球的速度Function_vf=vfx,vfy,vfz;%Function_vf碰撞时风速Function_ww=wx(i-1),wy(i-1),wz(i-1);%Function_ww碰撞时篮球的角速度Function_s=sx(i-1),sy(i-1),sz(i-1); %篮球撞击点的球心坐标 V,W,S,rot

23、ation_t=Impact_Rotation_2(Function_A,Function_AA,Function_vv,Function_vf,Function_ww,Function_s);%求解篮球碰撞过程中的旋转角速度及反弹速度(完全非弹性碰撞过程)vx(i-1)=V(1);%碰撞后篮球的速度vy(i-1)=V(2);vz(i-1)=V(3);wx(i-1)=W(1);%碰撞后篮球的角速度wy(i-1)=W(2);wz(i-1)=W(3);sx(i-1)=S(1);%碰撞后篮球球心坐标sy(i-1)=S(2);sz(i-1)=S(3);Impact_I12(3,I12)=rotatio

24、n_t;%记录篮球与篮圈撞击过程经历的时间信息Endif Impact_3=1 %判断篮球是否与地面碰撞Impact=3;%篮球与篮球场地面相撞的记录器 I3=I3+1;%判断篮球与地面撞击的信息Impact_I3(1,I3)=i-1;%记录篮球与地面碰撞的次数信息Impact_I3(2,I3)=(i-1-warning)*normal_t+warning*warning_t;%记录篮球与地面碰撞前经历的时间信息Function_A=0,1,0;%Function_A为撞击点的法向量Function_vv=vx(i-1),vy(i-1),vz(i-1);%Function_vv碰撞时篮球的速度

25、Function_vf=vfx,vfy,vfz;%Function_vf碰撞时风速Function_ww=wx(i-1),wy(i-1),wz(i-1);%Function_ww碰撞时篮球的角速度Function_s=sx(i-1),sy(i-1),sz(i-1); %篮球撞击点的球心坐标V,W,S,rotation_t=Impact_Rotation_1(Function_A,Function_vv,Function_vf,Function_ww,Function_s,Impact);%求解篮球碰撞过程中的旋转角速度及反弹速度(完全非弹性碰撞过程)vx(i-1)=V(1);%碰撞后篮球的速度

26、vy(i-1)=V(2);vz(i-1)=V(3);wx(i-1)=W(1);%碰撞后篮球的角速度wy(i-1)=W(2);wz(i-1)=W(3);sx(i-1)=S(1);%碰撞后篮球球心坐标sy(i-1)=S(2);sz(i-1)=S(3);Impact_I3(3,I3)=rotation_t;%记录篮球与篮圈撞击过程经历的时间信息endif Impact_4=1 %判断篮球是否与篮球架立柱碰撞Impact=4;%篮球与篮球架立柱相撞的记录器 I4=I4+1;%判断篮球与篮球架立柱撞击的信息Impact_I4(1,I4)=i-1;%记录篮球与篮球架立柱碰撞的次数信息Impact_I4(2

27、,I4)=(i-1-warning)*normal_t+warning*warning_t;%记录篮球与篮球架立柱碰撞前经历的时间信息x=bx+rc*(1/(bx2-2*bx*sx(i-1)+sx(i-1)2+sz(i-1)2)(1/2)*(bx-sx(i-1);%求解篮球与篮球架立柱碰撞的公切点y=sy(i-1);z=-rc*sz(i-1)*(1/(bx2-2*bx*sx(i-1)+sx(i-1)2+sz(i-1)2)(1/2);afx=x-bx;%篮球于篮球架立柱碰撞公切点处的法向量afy=0;afz=z;Function_A=afx,afy,afz;%Function_A为撞击点的法向量

28、Function_vv=vx(i-1),vy(i-1),vz(i-1);%Function_vv碰撞时篮球的速度Function_vf=vfx,vfy,vfz;%Function_vf碰撞时风速Function_ww=wx(i-1),wy(i-1),wz(i-1);%Function_ww碰撞时篮球的角速度Function_s=sx(i-1),sy(i-1),sz(i-1); %篮球撞击点的球心坐标V,W,S,rotation_t=Impact_Rotation_1(Function_A,Function_vv,Function_vf,Function_ww,Function_s,Impact

29、);%求解篮球碰撞过程中的旋转角速度及反弹速度(完全非弹性碰撞过程)vx(i-1)=V(1);%碰撞后篮球的速度vy(i-1)=V(2);vz(i-1)=V(3);wx(i-1)=W(1);%碰撞后篮球的角速度wy(i-1)=W(2);wz(i-1)=W(3);sx(i-1)=S(1);%碰撞后篮球球心坐标sy(i-1)=S(2);sz(i-1)=S(3);Impact_I4(3,I4)=rotation_t;%记录篮球与篮圈撞击过程经历的时间信息end vvx(i-1)=vx(i-1)-vfx;%(V(篮球速度)-VF(风速度)矢量的分量vvy(i-1)=vy(i-1)-vfy;vvz(i-

30、1)=vz(i-1)-vfz;vv(i-1)=(vvx(i-1)2+vvy(i-1)2+vvz(i-1)2)0.5;%(V(篮球速度)-VF(风速度)矢量的模re=pg*dp*vv(i-1)/u;%雷诺数if re>1000|re=0%气动阻力系数CD的计算(当re=0时，假设CD=0.44)CD=0.44;%Newton公式elseif re<1000&re>1CD=30/re0.625;%Allen公式elseCD=24/re;%Stokes公式endG=-m*g;%篮球受到的重力(y轴负向)F1=pg*Vb*g;%篮球受到的浮力（y轴正向)fd(i-1)=-CD

31、*pg*vv(i-1)*pi*dp2/8;%气动阻力系数(注意正负号)fx(i-1)=fd(i-1)*vvx(i-1);%气动阻力fy(i-1)=fd(i-1)*vvy(i-1);fz(i-1)=fd(i-1)*vvz(i-1);fm=-pi*pg*dp3/8;%马格努斯力的系数(注意正负号)fmx(i-1)=fm*(vvy(i-1)*wz(i-1)-vvz(i-1)*wy(i-1);%马格努斯力fmy(i-1)=-fm*(vvx(i-1)*wz(i-1)-vvz(i-1)*wx(i-1);fmz(i-1)=fm*(vvx(i-1)*wy(i-1)-vvy(i-1)*wx(i-1);ax(i-

32、1)=(fx(i-1)+fmx(i-1)/m;%篮球的加速度ay(i-1)=(G+F1+fy(i-1)+fmy(i-1)/m;az(i-1)=(fz(i-1)+fmz(i-1)/m;vx(i)=vx(i-1)+ax(i-1)*t;%下一步长篮球的三维速度分量vy(i)=vy(i-1)+ay(i-1)*t;vz(i)=vz(i-1)+az(i-1)*t;awx(i-1)=0;%篮球的角加速度(注意，此时表示运动过程中没有改变角速度的力矩存在)awy(i-1)=0;awz(i-1)=0;wx(i)=wx(i-1)+awx(i-1)*t;%下一步长篮球的角速度wy(i)=wy(i-1)+awy(i-

33、1)*t;wz(i)=wz(i-1)+awz(i-1)*t;sx(i)=sx(i-1)+vx(i-1)*t+ax(i-1)*t2/2;%下一步长篮球球心的位置坐标sy(i)=sy(i-1)+vy(i-1)*t+ay(i-1)*t2/2;sz(i)=sz(i-1)+vz(i-1)*t+az(i-1)*t2/2;endSxyz=sx;sy;sz;Vxyz=vx;vy;vz;Wxyz=wx;wy;wz;end附录二function V,W,S,rotation_t=Impact_Rotation_1(A,vv,vf,ww,s,Impact)%求解篮球仅与一个碰撞客体相碰撞过程中的旋转角速度及反弹速度

34、(完全非弹性碰撞过程)%A为撞击点的法向量,vv为碰撞时篮球的速度,vf为风速,ww为碰撞时篮球的角速度,s为篮球撞击点的球心坐标%Impact为篮球与碰撞物碰撞的信息(Impact=1,篮球与篮圈相撞；Impact=2,篮球与篮板相撞；Impact=3,篮球与篮球场地面相撞；Impact=4,篮球与篮球立柱相撞t=10-6;%时间迭代步长pg=1.185;%空气密度(298K)(国际单位制，下同)u=1.84*10-5;%空气黏度(298K)g=9.8;%重力加速度m=0.6;%篮球质量dp=0.246;%篮球直径I=0.5*m*(dp/2)2;%篮球的转动惯量Vp=4/3*pi*(dp/2

35、)3;%篮球体积G=-m*g+pg*g*Vp;%篮球的重力(包括浮力)A=A/norm(A);%法向量A的单位向量An=A(1);%单位法向量A在x,y,z轴中的分量Bn=A(2);Cn=A(3);vx=vv(1);%篮球速度vv在x,y,z轴中的分量vy=vv(2);vz=vv(3);vfx=vf(1);%风速在x,y,z轴中的分量vfy=vf(2);vfz=vf(3);wx=ww(1);%角速度ww在x,y,z轴中的分量wy=ww(2);wz=ww(3);sx=s(1);%篮球来时的球心坐标sy=s(2);sz=s(3);if vx*An+vy*Bn+vz*Cn>0%使来时速度与法向

36、量之间的夹角为钝角An=-An;Bn=-Bn;Cn=-Cn;endif norm(cross(A,vv)=0%篮球碰撞时的来时速度与撞击的法向量共线时，切向量的求解if An=0Ct=0.4;%假设切向量的Ct=0.4At=-(Cn*Ct-(Bn*(An*(-An2*Ct2+An2-Bn2*Ct2+Bn2-Cn2*Ct2)(1/2)+Bn*Cn*Ct)/(An2+Bn2)/An;Bt=-(An*(-An2*Ct2+An2-Bn2*Ct2+Bn2-Cn2*Ct2)(1/2)+Bn*Cn*Ct)/(An2+Bn2);elseif An=0&Bn=0Ct=0.4;%假设切向量的Ct=0.4A

37、t=(-Bn2*Ct2+Bn2-Cn2*Ct2)(1/2)/Bn;Bt=-(Cn*Ct)/Bn;elseif An=0&Bn=0At=1;Bt=0;Ct=0;endelse%篮球碰撞时的来时速度与撞击的法向量不共线时，切向量的求解vxyz=(vx2+vy2+vz2)0.5;%篮球速度vv的模nv=An*vx+Bn*vy+Cn*vz;%(单位向量A)·v%求解切向量t,t=X*n+Y*v;t·n=0,|t|=1;equation_1=sym('X+Y*nv=0');equation_2=sym('X2+2*X*Y*nv+Y2*vxyz2=1&#

38、39;);X,Y=solve(equation_1,equation_2,'X','Y')X=-nv*(-1/(nv2-vxyz2)(1/2);Y=(-1/(nv2-vxyz2)(1/2);At=X*An+Y*vx;%t切向向量Bt=X*Bn+Y*vy;Ct=X*Cn+Y*vz;ABCt=(At2+Bt2+Ct2)0.5;%t切向向量的模At=At/ABCt;%单位切向向量Bt=Bt/ABCt;Ct=Ct/ABCt;endAnt=Bn*Ct-Cn*Bt;%法向量叉乘切向量Bnt=-(An*Ct-Cn*At);Cnt=An*Bt-Bn*At;NT=(Ant2+Bn

39、t2+Cnt2)0.5;%叉向量的模Ant=Ant/NT;%叉向量的单位向量Bnt=Bnt/NT;Cnt=Cnt/NT;A0=0;0;0;%原坐标系的坐标原点(0,0,0)A00=sx;sy;sz;%(n,t,nt,)坐标系的坐标原点在原坐标系的坐标A01=0;G;0;%重力在原坐标系的矢量坐标(注意x、y、z的方向笛卡尔坐标系x向右，y向上，z向前)A02=vfx;vfy;vfz;%风速在原坐标系的矢量坐标A03=wx;wy;wz;%篮球来时角速度在原坐标系中的矢量坐标A04=vx;vy;vz;%篮球来时速度在原坐标系中的矢量坐标B=An,At,Ant;Bn,Bt,Bnt;Cn,Ct,Cnt

40、;%坐标系变换方程的系数矩阵，A00+B·C=AQ,R=qr(B);%对B进行QR分解AA0=inv(R)*Q'*(A0-A00);%原坐标系的坐标原点在新坐标系中的坐标AA01=inv(R)*Q'*(A01-A00);%新坐标系下的重力矢量坐标AA01=AA01-AA0;%新坐标系中的重力矢量Gn=AA01(1);%重力(包括浮力)在n法向方向的分力Gt=AA01(2);%重力(包括浮力)在t切向方向的分力Gnt=AA01(3);%重力(包括浮力)在nt切向方向的分力AA02=inv(R)*Q'*(A02-A00);%新坐标系下的风速的矢量坐标AA02=AA

41、02-AA0;%新坐标系下的风速的矢量AA03=inv(R)*Q'*(A03-A00);%新坐标系下的角速度矢量坐标AA03=AA03-AA0;%新坐标系下的角速度矢量AA04=inv(R)*Q'*(A04-A00);%新坐标系下的速度矢量坐标AA04=AA04-AA0;%新坐标系下的速度矢量AA05=AA02-AA04;%新坐标系下风速与篮球速度的矢量差kn=28300;%kn,kt分别为篮球的n法向和t,nt切向方向钢度kt=630;sn=0;%篮球在n法向、t切向和nt切向方向的弹性变形st=0;snt=0;switch Impactcase 1uf=0.45;%篮球与篮

42、圈间的摩擦系数kkn=12.4973;%kkn,kkt分别为篮球与篮圈间的n法向和(t,nt)切向方向的阻尼系数kkt=1.8646;case 2uf=0.35;%篮球与篮板间的摩擦系数kkn=10.5958;%kkn,kkt分别为篮球与篮板间的n法向和(t,nt)切向方向的阻尼系数kkt=1.5809;case 3uf=0.57;%篮球与篮球场地面间的摩擦系数kkn=14.4415;%kkn,kkt分别为篮球与篮圈架立柱间的n法向和(t,nt)切向方向的阻尼系数kkt=2.1547;case 4uf=0.4;%篮球与篮球架立柱间的摩擦系数kkn=12.4973;%kkn,kkt分别为篮球与地

43、面间的n法向和(t,nt)切向方向的阻尼系数kkt=1.8646;endnumber=0;%记录篮球的碰撞时间参数while sn<=0number=number+1;re=pg*dp*norm(AA05)/u;%初始时篮球的雷诺数if re>=1000%篮球气动阻力系数CD的计算CD=0.44;%Newton公式elseif re<1000&re>1CD=30/re0.625;%Allen公式elseCD=24/re;%Stokes公式end FDn=CD*pg*norm(AA05)*pi*dp2/8*AA05(1);%气动阻力在n法向方向的阻力(n向上，t向

44、右，nt向里,下同)FDt=CD*pg*norm(AA05)*pi*dp2/8*AA05(2);%气动阻力在t切向方向的阻力FDnt=CD*pg*norm(AA05)*pi*dp2/8*AA05(3);%气动阻力在nt切向方向的阻力FMn=pg*pi*dp3/8*(AA05(2)*AA03(3)-AA05(3)*AA03(2);%马格努斯力在n法向方向的分力FMt=-pg*pi*dp3/8*(AA05(1)*AA03(3)-AA05(3)*AA03(1);%马格努斯力在t切向方向的分力FMnt=pg*pi*dp3/8*(AA05(1)*AA03(2)-AA05(2)*AA03(1);%马格努斯

45、力在nt切向方向的分力vn=AA04(1);%接触点的法向速度(注意新坐标体系中的n法向向量、t切向向量和nt切向向量的坐标分别为(1,0,0)、(0,1,0)和(0,0,1)vst=AA04(2)+dp/2*AA03(3);%接触点的滑移速度,切向方向vsnt=AA04(3)-dp/2*AA03(2);%接触点的滑移速度,nt切向方向fn=-kn*sn-kkn*vn;%篮球n法向接触力ft=-kt*st-kkt*vst;%篮球t切向接触力fnt=-kt*snt-kkt*vsnt;%篮球nt切向方向的接触力if (ft2+fnt2)0.5>uf*abs(fn)%判断篮球在接触表面是否发生

46、滑动(注意fn为norm(fn)vs=vst,vsnt;%切向速度矢量vs=vs/norm(vs);%单位切向速度的方向向量ft=-uf*abs(fn)*vs(1);%t切向方向的接触力fnt=-uf*abs(fn)*vs(2);%nt切向方向的接触力endFn=fn+Gn+FDn+FMn;%n法向方向的力Ft=ft+Gt+FDt+FMt;%t切向方向的力Fnt=fnt+Gnt+FDnt+FMnt;%nt切向方向的力 an=Fn/m;%篮球的n法向加速度at=Ft/m;%篮球的t切向加速度ant=Fnt/m;%篮球的nt切向方向加速度sn=sn+AA04(1)*t+0.5*an*t2;%下一步

47、长的n法向弹性变形位移st=st+AA04(2)*t+0.5*at*t2;%下一步长的t切向弹性变形位移snt=snt+AA04(3)*t+0.5*ant*t2;%下一步长的nt切向方向弹性变形位移AA04(1)=AA04(1)+an*t;%下一步长的n法向速度AA04(2)=AA04(2)+at*t;%下一步长的t切向速度AA04(3)=AA04(3)+ant*t;%下一步长的nt切向方向速度AA05=AA02-AA04;%下一步长的风速与篮球速度的矢量差Mn=0;%篮球n法向方向的力矩等于0(注意,力矩M=r×F)Mt=fnt*dp/2;%篮球t切向方向的力矩Mnt=-ft*dp

48、/2;%篮球nt切向方向的力矩awn=Mn/I;%篮球n法向方向的角加速度awt=Mt/I;%篮球t切向方向的角加速度awnt=Mnt/I;%篮球nt切向方向的角加速度 AA03(1)=AA03(1)+awn*t;%下一步长的n法向角速度AA03(2)=AA03(2)+awt*t;%下一步长的t切向角速度AA03(3)=AA03(3)+awnt*t;%下一步长的nt切向角速度endrotation_t=number*t;%记录篮球碰撞经历的时间S=sn;st;snt;%新坐标系下的变形矢量S=A00;%+B*S;%原坐标系下的位移矢量(注意加上篮球来时的位移矢量)V=B*AA04;%原坐标系下

49、的速度矢量W=B*AA03;%原坐标系下的角速度矢量end附录三function V,W,S,rotation_t=Impact_Rotation_2(A,AA,vv,vf,ww,s)%求解篮球同时与篮圈和篮板相碰撞过程中的旋转角速度及反弹速度(完全非弹性碰撞过程)%A为篮球与篮圈在撞击点处的法向量,AA为篮球与篮板在撞击点处的法向量，vv为碰撞时篮球的速度,vf为风速,ww为碰撞时篮球的角速度,s为篮球撞击点的球心坐标t=10-6;%时间迭代步长pg=1.185;%空气密度(298K)(国际单位制，下同)u=1.84*10-5;%空气黏度(298K)g=9.8;%重力加速度m=0.6;%篮球

50、质量dp=0.246;%篮球直径I=0.5*m*(dp/2)2;%篮球的转动惯量Vp=4/3*pi*(dp/2)3;%篮球体积G=-m*g+pg*g*Vp;%篮球的重力(包括浮力)A=A/norm(A);%法向量A的单位向量An=A(1);%单位法向量A在x,y,z轴中的分量Bn=A(2);Cn=A(3);AA=AA/norm(AA);%法向量AA的单位向量AAn=AA(1);%单位法向量AA在x,y,z轴中的分量BBn=AA(2);CCn=AA(3);vx=vv(1);%篮球速度vv在x,y,z轴中的分量vy=vv(2);vz=vv(3);vfx=vf(1);%风速在x,y,z轴中的分量vf

51、y=vf(2);vfz=vf(3);wx=ww(1);%角速度ww在x,y,z轴中的分量wy=ww(2);wz=ww(3);sx=s(1);%篮球来时的球心坐标sy=s(2);sz=s(3);if vx*An+vy*Bn+vz*Cn>0%使来时速度与法向量A之间的夹角为钝角(构建坐标体系1)An=-An;Bn=-Bn;Cn=-Cn;endif norm(cross(A,vv)=0%篮球碰撞时的来时速度与撞击的法向量A共线时，切向量的求解(构建坐标体系1)if An=0 Ct=0.4;%假设切向量At=-(Cn*Ct-(Bn*(An*(-An2*Ct2+An2-Bn2*Ct2+Bn2-Cn

52、2*Ct2)(1/2)+Bn*Cn*Ct)/(An2+Bn2)/An;Bt=-(An*(-An2*Ct2+An2-Bn2*Ct2+Bn2-Cn2*Ct2)(1/2)+Bn*Cn*Ct)/(An2+Bn2);elseif An=0&Bn=0Ct=0.4;%假设切向量At=(-Bn2*Ct2+Bn2-Cn2*Ct2)(1/2)/Bn;Bt=-(Cn*Ct)/Bn;elseif An=0&Bn=0At=1;Bt=0;Ct=0;endelse%篮球碰撞时的来时速度与撞击的法向量不共线时，切向量的求解vxyz=(vx2+vy2+vz2)0.5;%篮球速度vv的模nv=An*vx+Bn*v

53、y+Cn*vz;%(单位向量A)·v%求解切向量X=-nv*(-1/(nv2-vxyz2)(1/2);Y=(-1/(nv2-vxyz2)(1/2);At=X*An+Y*vx;%t切向向量Bt=X*Bn+Y*vy;Ct=X*Cn+Y*vz;ABCt=(At2+Bt2+Ct2)0.5;%t切向向量的模At=At/ABCt;%单位切向向量Bt=Bt/ABCt;Ct=Ct/ABCt;endB=At,Bt,Ct;Ant=Bn*Ct-Cn*Bt;%法向量叉乘切向量(构建坐标体系1)Bnt=-(An*Ct-Cn*At);Cnt=An*Bt-Bn*At;NT=(Ant2+Bnt2+Cnt2)0.5;

54、%叉向量的模Ant=Ant/NT;%叉向量的单位向量Bnt=Bnt/NT;Cnt=Cnt/NT;C=Ant,Bnt,Cnt;if vx*AAn+vy*BBn+vz*CCn>0%使来时速度与法向量AA之间的夹角为钝角(构建坐标体系2)AAn=-AAn;BBn=-BBn;CCn=-CCn;endif norm(cross(AA,vv)=0%篮球碰撞时的来时速度与撞击的法向量AA共线时，切向量的求解(构建坐标体系2)if AAn=0CCt=0.4;%假设切向量AAt=-(CCn*CCt-(BBn*(AAn*(-AAn2*CCt2+AAn2-BBn2*CCt2+BBn2-CCn2*CCt2)(1

55、/2)+BBn*CCn*CCt)/(AAn2+BBn2)/AAn;BBt=-(AAn*(-AAn2*CCt2+AAn2-BBn2*CCt2+BBn2-CCn2*CCt2)(1/2)+BBn*CCn*CCt)/(AAn2+BBn2);elseif AAn=0&BBn=0CCt=0.4;%假设切向量AAt=(-BBn2*CCt2+BBn2-CCn2*CCt2)(1/2)/BBn;BBt=-(CCn*CCt)/BBn;elseif AAn=0&BBn=0AAt=1;BBt=0;CCt=0;endelse%篮球碰撞时的来时速度与撞击的法向量AA不共线时，切向量的求解vxyz=(vx2+vy2+vz2)0.5;%篮球速度vv的模nv=AAn*vx+BBn*vy+CCn*vz;%(单位向量A)·v%求解切向量X=-nv*(-1/(nv2-vxyz2)(1/2);Y=(-1/(nv2-vxyz2)(1/2);AAt=X*AAn+Y*vx;%t切向向量BBt=X*BBn+Y*vy;CCt=X*CCn+Y*vz;ABCt=(AAt2+BBt2+CCt2)0.5;%t切向向量的模AAt=AAt/ABCt;%单位切向向量BB