2020年第九届数学建模国际赛小美赛
A题 自由泳
原题再现：
在所有常见的游泳泳姿中，哪一种最快？哪个冲程推力最大？在自由泳项目中，游泳者可以选择他们的泳姿，他们通常选择前面的爬行。然而，游泳运动员是否有可能进一步提高他们的划水速度？请用一个数学模型来分析中风和指导运动员的训练。
整体求解过程概述(摘要)
游泳是生活中常见的一种运动形式，那么，在各种泳姿中，哪种泳姿能获得最大的速度或最大的推进力呢？是否有可能在爬行的基础上进行改进，进一步提高游泳速度？本文结合水动力公式对上述问题进行了讨论。
对于问题一，我们首先从几何上抽象出人体的形状，将手臂视为一个统一的圆柱体，将腿部视为两个圆柱体的组合。然后对四种泳姿的动作进行分析。根据流体力学中的阻力公式，将动阻力和静阻力分别表示，其中动阻力假定为全部推进力，腿的推进力作为补偿项。然后，根据牛顿第二定律的公式，建立了行驶速度随运动频率和时间变化的偏微分方程。在求解方程时，以运动频率为常数，将偏微分方程转化为常微分方程。然后利用四阶Runge-Kutta方法求解该方程，得到旅行速度与时间的定量关系。经过比较，我们发现自由泳是四种游泳方式中速度最快的。在30s内向前距离可达53.852m，比最慢蛙泳远14m左右。在推进力比较中，蝶泳在稳态下最大推进力可达150N左右，比最小仰泳大95N左右。
对于问题二，首先，在人类几何的抽象中，我们将手臂细化为上臂、前臂和手。以运动频率为常数，以上臂与肩线夹角、前臂与水平面夹角为变量，建立了运动速度与时间的偏微分方程。为了便于求解，我们假设两个角度在一次游泳中保持不变。然后将偏微分方程转化为常微分方程，采用固定步长搜索法遍历各角度求解，得到当两个角度均为0.5∏时的最大旅行速度，即当手臂伸直且手臂运动平面垂直于水平面时，可以得到最大旅行速度。最后，我们发现人体可以在27秒内以最大速度前进58.25米。为了保证模型的准确性，我们还进行了敏感性分析，证明了直臂划水是不同身体特征运动员的最佳划水。此外，我们还考虑了手指夹角对移动速度的影响，确定了当双手合拢时可以获得最大速度。
最后，我们将我们的结果汇编成一份提案，尽我们所能帮助运动员，我们期待着他们继续取得成功！
模型假设：
为了简化给定的问题，我们对我们的模型做出以下假设：
1、假设人体部分的某种几何近似对模型结果的影响可以忽略不计。
2、假设游泳时惯性阻力和空气阻力可以忽略不计。
3、假设动态阻力在游泳过程中提供完全推进力。
4、假设在游泳过程中身体各部分之间的角度保持不变，并且在一次游泳过程中人体运动的频率不变。
5、仰泳腿的动作可视为与自由泳相同，由于人体关节的限制，手不能达到与自由泳相同的运动幅度。因此，我们参考文献，并考虑由半臂冲程提供的身体前半部分的推进力。
问题重述：
问题背景
游泳是众所周知的一项竞技运动。在游泳项目中，运动员经常被要求使用特定的游泳方式进行比赛，在一定距离内能花最短时间的运动员将获胜。在游泳比赛中，目前常见的游泳方式有四种：爬泳、蛙泳、蝶泳和仰泳，它们都要求不同的动作和标准。在奥运会上，这四种风格都是合理的游泳方式，具有不同的竞技要求。
现在，有必要弄清楚四种方法中哪种最快，哪种游泳方法可以获得最大的推进力。此外，在自由泳项目中，由于没有规定的游泳泳姿，运动员往往选择爬行。那么，在自由泳中，是否有可能进一步提高泳姿，从而达到更快的旅行速度？这是本文要解决的问题。
此外，我们还有光荣的使命。通过数学模型对问题进行分析，提出改进游泳风格的方案，并提出建议，供运动员训练时参考，以达到更好的效果！
我们的工作
通过我们的游泳知识和流体动力学分析，我们开发了几个数学模型，可以说明该主题所需的问题。具体如下：
•通过对人体在四种泳姿中不同运动的几何近似，将手臂和手视为一个统一的圆柱体，腿和脚视为圆柱体，根据流体力学中的流体阻力公式计算四种泳姿的推进力，并根据推进力和阻力公式结合牛顿第二定律求解偏微分方程，得到各泳姿速度与行程时间的关系。这将导致最推进和最快的冲程。
•在上述结论的基础上，我们对爬行动作进行了较为详细的分解，将手臂分别视为上臂、前臂和手掌分别视为两个圆柱体和一个椭圆形圆柱体，设置游泳过程中手臂之间的夹角，并求出解出的方程的极值，得到能使速度最大化的身体角度。
•此外，我们推测手掌的手指张开度也会对有效速度产生影响，为此，我们模拟了与手指张开角度相关的推进力。根据该模型，我们可以得到最合适的角度。
•最后，我们将我们的研究成果汇编成一份详细的方案，供运动员在训练期间参考，希望该方案能帮助他们取得更好的效果。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
function dVb=he(t,Vb,gama,beita)
D1=0.083;
La1=0.34;
D2=0.054;
La2=0.3;
Ah=0.019;
Lh=0.19;
Ka=400;
Kb=31.487;
Kh=900;
f=1.11;
xita=mod(2*pi*f*t,pi);
dVb=((Ka*D1*(4/3*pi^2*f^2*La1^3*sin(gama)^2-2*pi*f*La1^2*Vb*sin(xita)*sin(gama)+
La1*Vb^2*sin(xita)^2)+Ka*D2*(4/3*pi^2*f^2*La2^3*sin(beita)^2+4*pi^2*f^2*La2^2*La1*
sin(beita)*sin(gama)+4*pi^2*f^2*La1^2*La2*sin(gama)^2-2*pi*f*La2^2*Vb*sin(beita)
*sin(xita)-4*pi*f*La1*La2*Vb*sin(gama)*sin(xita)+La2*Vb^2*sin(xita)^2)+Kh*Ah*
(4*pi^2*f^2*(La1*sin(gama)+La2*sin(beita)+Lh*sin(beita)/2)^2-4*pi*f*
(La1*sin(gama)+La2*sin(beita)+Lh*sin(beita)/2)*Vb*sin(xita)+Vb^2*sin(xita)^2 ))
*sin(xita)-sign(Vb)*Kb*Vb^2)/60;
end
S1=zeros(numel(0:pi/2/90:pi),numel(0:pi/2/90:pi));
max_S=0;
max_gama=0;
max_beita=0;
for gama=0:pi/2/90:pi
for beita=0:pi/2/90:pi
[t11,y11]=ode45(@(t,Vb)he(t,Vb,gama,beita),[0,27],0);
plot(t11,y11);
X1=[t11 y11];
S1(round(gama/(pi/2/90)+1),round(beita/(pi/2/90)+1))
=(numel(0:27)-1)/numel(X1(:,1))*sum(X1(:,2));
if(S1>=max_S)
max_S=S1;
max_gama=gama;
max_beita=beita;
end
end
end
max_gama=max_beita/pi*180;
max_beita=max_beita;
gama1=0:pi/2/90:pi;
beita1=0:pi/2/90:pi;
gama1=gama1/pi*180;
beita1=beita1/pi*180;
[gama2,beita2]=meshgrid(0:0.2:180,0:0.2:180);
S11=griddata(gama1,beita1,S1,gama2,beita2)
shading interp
mesh(gama2,beita2,S11)
clc
clear
T_t=10;
%f=1.667
%T=0.6s
T=0.9;
f=1.11;
A1=50.735;
A2=31.487;
A3=25.735;
L1=(533+90)/1000;%
%sign(2*pi*f*L1-y*sin(2*pi*f*t))*
L2=878/1000;
%%
% T_xs0=0;A_xs0=0;
% T_xs1=-0.1;A_xs1=0;
% T_xs2=-0.05;A_xs2=0;
% T_xs3=0.05;A_xs3=0;
% T_xs4=0.1;A_xs4=0;
%%
T_xs0=0;A_xs0=0;
T_xs1=0;A_xs1=0.2;
T_xs2=0;A_xs2=0.1;
T_xs3=0;A_xs3=-0.1;
T_xs4=0;A_xs4=-0.2;
[t0,y0]=ode45(@(t,y)free_Fenxi_fucntion(t,y,T_xs0,A_xs0),[0,T_t],0);
[t1,y1]=ode45(@(t,y)free_Fenxi_fucntion(t,y,T_xs1,A_xs1),[0,T_t],0);
[t2,y2]=ode45(@(t,y)free_Fenxi_fucntion(t,y,T_xs2,A_xs2),[0,T_t],0);
[t3,y3]=ode45(@(t,y)free_Fenxi_fucntion(t,y,T_xs3,A_xs3),[0,T_t],0);
[t4,y4]=ode45(@(t,y)free_Fenxi_fucntion(t,y,T_xs4,A_xs4),[0,T_t],0);
%%
X1=[t0 y0];
S0=(numel(0:T_t)-1)/numel(X1(:,1))*sum(X1(:,2));
S00=zeros(1,numel(t0));
S00(1)=y0(1);
for i0=2:numel(t0)
S00(i0)=S00(i0-1)+(numel(0:T_t)-1)/numel(X1(:,1))*y0(i0);
end
%%
X1=[t1 y1];
S1=(numel(0:T_t)-1)/numel(X1(:,1))*sum(X1(:,2));
S11=zeros(1,numel(t1));
S11(1)=y1(1);
for i11=2:numel(t1)
S11(i11)=S11(i11-1)+(numel(0:T_t)-1)/numel(X1(:,1))*y1(i11);
end
%%
X2=[t2 y2];
S2=(numel(0:T_t)-1)/numel(X2(:,1))*sum(X2(:,2));
S22=zeros(1,numel(t2));
S22(1)=y2(1);
for i22=2:numel(t2)
S22(i22)=S22(i22-1)+(numel(0:T_t)-1)/numel(X2(:,1))*y2(i22);
end
%%
X3=[t3 y3];
S3=(numel(0:T_t)-1)/numel(X3(:,1))*sum(X3(:,2));
S33=zeros(1,numel(t3));
S33(1)=y3(1);
for i33=2:numel(t3)
S33(i33)=S33(i33-1)+(numel(0:T_t)-1)/numel(X3(:,1))*y3(i33);
end
%%
X4=[t4 y4];
S4=(numel(0:T_t)-1)/numel(X4(:,1))*sum(X4(:,2));
S44=zeros(1,numel(t4));
S44(1)=y4(1);
for i44=2:numel(t4)
S44(i44)=S44(i44-1)+(numel(0:T_t)-1)/numel(X4(:,1))*y4(i44);
end
D_fenxi_picture(t0, S00, t1, S11, t2, S22, t3, S33, t4, S44, y0,y1,y2,y3,y4)