2019年认证杯SPSSPRO杯数学建模
纸飞机在飞行状态下的运动模型
A题 好风凭借力，送我上青云
原题再现：
纸飞机有许多种折法。世界上有若干具有一定影响力的纸飞机比赛，通常的参赛规定是使用一张特定规格的纸，例如 A4 大小的纸张，折成一架纸飞机。大多数比赛都不允许使用剪刀和胶水，有的比赛中可以少量使用胶带以调整纸飞机的重心。折好后在无风的室内环境投掷，并测量其留空时间和飞行距离等成绩。在 2012 年，一种叫做 Suzanne 的折法曾创下飞行距离的吉尼斯世界纪录，飞行了超过 69 米。Suzanne 的折法参见：https://v.qq.com/x/page/o0636km20dz.html
第一阶段问题： 如果纸飞机 (以 Suzanne 为例) 的翼面是精确的平面，不允许有弯曲，请建立合理的数学模型，根据纸飞机的三视图和重心位置，计算投掷出手的合理角度，以使飞行距离尽可能远。
整体求解过程概述(摘要)
世界上有若干个具有一定影响力的纸飞机比赛，其中有一种叫 Suzanne 的折法的纸飞机创造了纸飞机飞行距离最远的世界纪录。本文通过对纸飞机在飞行状态下的受力分析，通过纸飞机的三视图与重心位置，计算出了纸飞机飞行最远距离时的最佳投掷角度，同时分析重力位置对纸飞机飞行状态的影响，简要概述了纸飞机稳定性较好时的重心位置，同时对模型进行了推广应用。
首先，对一些基本的斜抛运动做分析，通过运用牛顿力学定律，建立了理想状态下的斜抛运动模型，计算出了在不考虑纸飞机初始掷出高度与阻力的情况下的最佳的投射角度为 45°。
接着，考虑了纸飞机掷出时存在初始高度的情况，改进了理想状态下的斜抛运动模型，通过运用实际的数据进行计算，得到了在不同初始高度与初始速度的情况下飞机的掷出的最佳角度与纸飞机的最远飞行距离，详细数据见附录。然后，又加入了空气对纸飞机存在浮力与阻力的因素，改进了前面两个模型，建立了非理想状态下的斜抛运动模型，通过对不同阻力系数的研究，得到了阻力系数对纸飞机飞行最远距离的影响，同时也求出了不同阻力系数下的最佳掷出角度，详细数据见正文。
最后，又分析了飞机自身重心位置与空气对纸飞机作用力不在同一位置的情况，通过建立方程组，进一步优化了飞机在非理想状态下的斜抛运动模型。
本文的优点在于通过简单的运动模型，逐渐深化改进，完善模型，最终建立一个比较完整的纸飞机运动模型。同时考虑了实际的情况下纸飞机的不同受力情况，使得建立的模型真实有效。最后对模型进行了推广与应用，在生产生活的其他方面具有较大的参考意义。
问题分析：
飞机一开始经过人对它的做功，使其产生一个初速度，接着在飞机飞行过程中，人对其不产生作用力，此时的作用力可以由中学物理知识得到，飞机在飞行状态中只受自身的重力与空气对它的阻力还有一部分升力。由于飞机是在无风的测试环境下进行飞行测试，所以在解题过程中不考虑复杂的阻力环境。
与真正的飞机相比，纸飞机是没有提供动力系统的装置，所以纸飞机在整个飞行状态下是呈滑翔式的飞行。纸飞机机翼与尾翼的形状在飞机的飞行状态下有着一定的重要作用，如果飞机的侧翼的面积很小，在飞行过程中纸飞机与空气的接触面积就会很小，所受的阻力同时也会减小，在下降的过程中下落的速度就会变快，落地时间减短，即使初速度比较大，但是飞机飞行的距离还是不会太远。所以飞机机翼的接触面积直接对飞机的飞行距离产生影响，纸飞机机翼面积越大，飞行时间可能会越长，飞行距离就会越远。
其次影响纸飞机飞行的因素还有飞机的重心点的位置。如果纸飞机的重心位置合适， 在其飞机飞行的过程中，机翼就会受到比较均匀的空气阻力，保证了飞机飞行的稳定性，才能飞得更久。机翼在在飞机的质量上具有很大的贡献值，所以机翼在机身上的位置就几乎可以决定纸飞机重心点的位置，如果重心太过靠前，空气阻力对与飞机的作用点为飞机的几何重心上，两个作用点就会有很大的差别，作用效果不会抵消，反而叠加，从而有可能加剧了飞机的下落速度，飞机在滑翔时间就会减短。如果能控制机翼的位置使重心点位置准确，就会增加纸飞机的稳定性。由于飞机在叠的过程中前后的纸张叠加，导致飞机的重心会在其几何中心之前，这必然会导致受力不均，由于飞机的对称性，飞机的重心点在中轴线上。
飞机在初始状态下的初始角度对飞机的飞行状态具有一定的影响。如果飞机的角度过大，可能倒是在竖直方向上的速度会变大，水平方向上的速度就小，虽然飞机飞行时间可能会增加，但是水平方向上的位移也不会很大。如果飞机初始角度很小，在竖直方向上的速度小，飞机在空中的飞行时间变短，即使速度很大，飞行的距离也不会太远。所以飞机的飞行角度一定在这个角度范围内会由最优的角度。
纸飞机的质量也与飞机飞行的时间有关系，在考虑阻力的情况下，如果飞机质量越大，在飞行过程中，重力就可能是阻力的多倍，此时飞机向下加速度要比质量小时候的加速度大。飞机飞行时间就会变短，飞机也不会飞太远。
模型假设：
（1）.假设纸飞机在飞行过程中，不受任何其他风力的影响，无不确定因素的影响。
（2）.假设纸飞机飞行路线的投影为一条直线，不存在拐弯飞行的状况。
（3）.假设纸飞机所受的空气对它的受力大小与速度成正比，比例系数不变。
（4）.假设重力加速度 g 恒为 9.8m/s^2 。
（5）.受力分析过程符合基本的牛顿定律。
论文缩略图：
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
%最佳角度
h=1.7:0.05:2.1;
h'
v=10:1:30;
a=asind(v(1)./sqrt((2*(v(1)*v(1)+9.8*h))))
a'
a=asind(v(2)./sqrt((2*(v(2)*v(2)+9.8*h))))
a'
a=asind(v(3)./sqrt((2*(v(3)*v(3)+9.8*h))))
a'
a=asind(v(4)./sqrt((2*(v(4)*v(4)+9.8*h))))
a'
a=asind(v(5)./sqrt((2*(v(5)*v(5)+9.8*h))))
a'
a=asind(v(6)./sqrt((2*(v(6)*v(6)+9.8*h))))
a'
a=asind(v(7)./sqrt((2*(v(7)*v(7)+9.8*h))))
a'
a=asind(v(8)./sqrt((2*(v(8)*v(8)+9.8*h))))
a'
a=asind(v(9)./sqrt((2*(v(9)*v(9)+9.8*h))))
a'
a=asind(v(10)./sqrt((2*(v(10)*v(10)+9.8*h))))
a'
a=asind(v(11)./sqrt((2*(v(11)*v(11)+9.8*h))))
a'
a=asind(v(12)./sqrt((2*(v(12)*v(12)+9.8*h))))
a’
a=asind(v(13)./sqrt((2*(v(13)*v(13)+9.8*h))))
a’
a=asind(v(14)./sqrt((2*(v(14)*v(14)+9.8*h))))
a’
a=asind(v(15)./sqrt((2*(v(15)*v(15)+9.8*h))))
a’
a=asind(v(16)./sqrt((2*(v(16)*v(16)+9.8*h))))
a’
a=asind(v(17)./sqrt((2*(v(17)*v(17)+9.8*h))))
a’
a=asind(v(18)./sqrt((2*(v(18)*v(18)+9.8*h))))
a’
a=asind(v(19)./sqrt((2*(v(19)*v(19)+9.8*h))))
a’
a=asind(v(20)./sqrt((2*(v(20)*v(20)+9.8*h))))
a’
a=asind(v(21)./sqrt((2*(v(21)*v(21)+9.8*h))))
a’
%最远距离
s=v/9.8;
s1=s.*sqrt(v.*v+2*9.8*h(1))
s1=s.*sqrt(v.*v+2*9.8*h(2))
s1=s.*sqrt(v.*v+2*9.8*h(3))
s1=s.*sqrt(v.*v+2*9.8*h(4))
s1=s.*sqrt(v.*v+2*9.8*h(5))
s1=s.*sqrt(v.*v+2*9.8*h(6))
s1=s.*sqrt(v.*v+2*9.8*h(7))
s1=s.*sqrt(v.*v+2*9.8*h(8))
s1=s.*sqrt(v.*v+2*9.8*h(9))
%不同阻力系数运动轨迹
h=2.1;
m=0.1;
g=9.8;
v=20;
th=40;
k=[0.01,0.1,0.5,1,1.5];
x=0:0.01:2;
y=h+((m*m*g)/(k(1)*k(1))+(m*v*sin(th))/k(1))*(k(1)*x)/(m*v*cos(th))+(m*m*g*log(1-(k(1)*x)/(m*v*c
os(th))))/(k(1)*k(1));
plot(x,y,'r');
hold on;
y2=h+((m*m*g)/(k(2)*k(2))+(m*v*sin(th))/k(2))*(k(2)*x)/(m*v*cos(th))+(m*m*g*log(1-(k(2)*x)/(m*v*
cos(th))))/(k(2)*k(2));
plot(x,y2,'b');
hold on;
y3=h+((m*m*g)/(k(3)*k(3))+(m*v*sin(th))/k(3))*(k(3)*x)/(m*v*cos(th))+(m*m*g*log(1-(k(3)*x)/(m*v*
cos(th))))/(k(3)*k(3));
plot(x,y3,'g');
hold on;
y4=h+((m*m*g)/(k(4)*k(4))+(m*v*sin(th))/k(4))*(k(4)*x)/(m*v*cos(th))+(m*m*g*log(1-(k(4)*x)/(m*v*
cos(th))))/(k(4)*k(4));
plot(x,y4,'y');
hold on;
y5=h+((m*m*g)/(k(5)*k(5))+(m*v*sin(th))/k(5))*(k(5)*x)/(m*v*cos(th))+(m*m*g*log(1-(k(5)*x)/(m*v*
cos(th))))/(k(5)*k(5));
plot(x,y5,'k')
title('不同阻力系数运动轨迹');
xlabel('水平距离 x/m');
ylabel('高度 y/m');
legend({'k=1.5','k=1','k=0.5','k=0.1','k=0.01'})
axis([0,1.9,0,2.2])
%阻力系数与最佳角度关系
x=[0 0.01 0.1 0.5 1 1.5 2];
th=[45 44.98 44.43 40.36 36.94 34.13 32.04];
plot(x,th);
xlabel('阻力系数 k/kg/s');
ylabel('最佳角度/°');
title('阻力系数与最佳角度关系');
grid on
