2022年第十一届数学建模国际赛小美赛A题翼龙如何飞行解题全过程文档及程序

2022年第十一届数学建模国际赛小美赛

A题翼龙如何飞行

原题再现：

翼龙是翼龙目中一个已灭绝的飞行爬行动物分支。它们存在于中生代的大部分时期：从三叠纪晚期到白垩纪末期。翼龙是已知最早进化出动力飞行的脊椎动物。它们的翅膀是由皮肤、肌肉和其他组织膜形成的，这些组织从脚踝延伸到显著延长的第四根手指[1]。
翼龙有两种主要类型。基底翼龙是体型较小的动物，通常有全齿颚和长长的尾巴。它们宽阔的翅膜可能包括并连接后腿。在地面上，它们会有一个尴尬的伸展姿势，但它们的关节解剖结构和强壮的爪子会使它们成为有效的攀爬者，而且它们可能生活在树上。基生翼龙是小型脊椎动物的食虫动物或捕食者。后来翼龙（翼龙目）进化出许多体型、形状和生活方式。翼龙的翅膀较窄，后肢自由，尾巴高度缩小，脖子较长，头部较大。在地面上，翼龙的四肢都走得很好，姿态直立，后脚直立，翼指向上折叠，在三指“手”上行走。化石轨道表明，至少有一些物种能够奔跑、涉水或游泳[2]。
翼龙身上有一层被称为比重瓶纤维的毛发状细丝，覆盖着它们的身体和翅膀的一部分[3]。在生活中，翼龙会有光滑或蓬松的皮毛，与鸟类的羽毛不同。早期的说法是，翼龙在很大程度上是冷血动物，像现代蜥蜴一样从环境中获得温暖，而不是燃烧卡路里。然而，后来的研究表明，它们可能是温血（吸热）、活跃的动物。呼吸系统使用气囊进行有效的单向“流通”呼吸，气囊将骨骼挖空到极致。翼龙的成年体型多种多样，从非常小的无颌类到已知最大的飞行生物，包括翼展至少9米的Quetzalcoatlus和Hatzegopteryx[4][5]。吸热、良好的氧气供应和强壮的肌肉相结合，使翼龙成为强大而能干的飞行者。
翼龙飞行的机制目前还没有完全理解或建模。佐藤胜富美利用现代鸟类进行了计算，得出结论认为翼龙不可能在高空停留[6]。在《翼龙的姿态、运动和古生态学》一书中，有理论认为翼龙之所以能够飞行，是因为白垩纪晚期富含氧气、稠密的大气层[7]。然而，佐藤和《翼龙的姿态、运动和古生态学》一书的作者的研究都基于翼龙类海鸟的过时理论，而体型限制不适用于陆生翼龙，如阿志达翼龙和带翅翼龙。此外，Darren Naish得出结论，巨型翼龙不需要现在和中生代之间的大气差异[8]。
另一个很难理解的问题是翼龙是如何起飞的。如果翼龙是冷血动物，那么体型巨大、冷血动物代谢效率低下的大型翼龙如何能够实现类似鸟类的起飞策略，只使用后肢产生推力，从而升空，目前尚不清楚。后来的研究表明，它们是温血动物，拥有强大的飞行肌肉，并使用飞行肌肉作为四足动物行走[9]。朴茨茅斯大学的Mark Witton和约翰·霍普金斯大学的Mike Habib认为翼龙使用跳跃机制来获得飞行[10]。它们翅膀前肢的巨大力量将使它们能够轻松起飞[9]。翼龙一旦升空，速度可达120公里/小时，飞行数千公里[10]。
您的团队被要求根据化石测量，为至少一种大型翼龙的飞行过程建立一个合理的数学模型，并回答以下问题。

1.对于您选择的翼龙物种，估计其正常飞行期间的平均速度。

2.对于你选择的翼龙物种，估计它在正常飞行中的翅膀拍打频率。

3.研究大型翼龙是如何起飞的；它们有可能像平地或水面上的鸟一样起飞吗？定量解释原因。

整体求解过程概述(摘要)

翼龙的飞行原理一直困扰着进化生物学家。巨大的面积，再加上化石保存不善，使得对翼龙飞行的解释成了问题，而且常常引起争议。
为了研究翼龙飞行的主题，我们创建了三个空气动力学模型，如下所示：

模型一：滑动模型；
模型二：机翼颤振模型；
模型三：四足发射起飞模型。

在模块I中，我们使用运动方程（EOM）来描述翼龙向下滑翔的动力学过程。该模型的约束条件主要是翼龙的重量和翼尖距地面或海面的高度。考虑到直接求解以函数为变量的优化问题具有挑战性，本文采用配点法。我们建立了步长为1s的离散点来描述翼龙在每个时刻的飞行状态。产生的滑动轨迹是一条平滑的曲线。
在模块二中，根据EOMs和叶片单元理论，建立了基于颤振角和扭转角的机翼颤振模型。考虑到周期边界飞行姿态的约束，我们假设翼龙将上升到其原始高度。我们计算了翼龙在拍动翅膀时的运动状态的变化。
基于相关文献和上述模型，本文提出了一个包含滑翔和颤振起飞的翼龙两级飞行模型。模型中的一个周期被假定为由滑动和连续拍打组成。通过计算机模拟，我们确定翼龙和格查尔亚特龙的平均速度分别为38.65m/s和42.49m/s，拍动频率分别为1.06Hz和0.83Hz。
在模块III中，我们假设以下起飞程序：翼龙使用四足发射，然后在达到一定高度后利用空气动力爬升。建立参数后，模型模拟翼龙拍打翅膀和起飞的过程，评估其飞行趋势，判断其是否具备起飞的必要条件。
模拟结果表明，翼龙和格查尔阿特拉斯四足动物在平地上发射后，可以借助上升气流起飞。而由于支撑力不足，水面阻力大，使得格查尔亚特鱼几乎不可能在水面上起飞，而翼龙则可以在特定角度上起飞。
为了测试对空气密度的敏感性，我们还评估了在高空气密度环境中的飞行性能。总体而言，我们的结果对空气密度不确定性具有鲁棒性。分析表明，该模型对升力系数等关键参数具有较强的鲁棒性。

模型假设：

通过对问题的全面分析，为了简化我们的模型，我们做出以下合理的假设。

•假设1：所有相关文献中的翼龙化石数据和根据化石估计的身体数据都是可靠的，并且在允许的误差范围内。
理由：所有化石数据都经过精确测量，所得数据可靠。现代三维重建方法[7]可以通过三维建模恢复恐龙的结构和外观，得到准确的恢复数据。

•假设2：风向沿y轴，风速仅取决于海拔高度
理由：在离地面不到400米的天空中，风向可以近似地视为与地面平行，风速取决于海拔高度[8]。然后我们可以人为地把风向作为y轴的正方向。

•假设3：在大型翼龙的正常飞行过程中，翼龙使用一种飞行策略，先拍打翅膀升起，然后滑动一段距离。以一次上升和一次下滑为周期，周期开始和结束时的高度基本不变。
理由：现存最大的鸟——秃鹫gryphus，翼展为5米，使用与上述类似的飞行策略，这是节省能源和增加飞行距离的最佳方式[9]。

•假设4：在滑翔和拍打阶段，关节点的姿态移动是瞬时的。
理由：扑翼阶段的过程比较复杂。翼龙会调整它的姿态来完成上升。这种姿态变化可视为飞行过程中的突然变化。

问题重述：

问题背景

翼龙是一种已经灭绝的爬行动物，有近100种。尽管翼龙和恐龙生活在同一个时代，但它们不是恐龙。
因为它们太大了，所以关于它们是否不能飞，或者如果能够飞，它们可能的功率能力、起飞方法和飞行特性可能是什么，也有争论。
本课题组选取了两种具有代表性的翼龙——无齿翼龙和无齿翼龙作为研究对象。Quetzalcatlus是一种存在于晚白垩世晚期的翼手目动物，是已知最大的飞行动物。无齿翼龙是翼展为7-9米的翼龙，生活在公元前70-65百万年前的晚白垩世[1]。
根据化石测量数据和建立的模型，可以更合理地量化两种翼龙的飞行能力。然后测定了这两种翼龙在正常飞行时的平均速度和翼瓣频率，判断它们是否有可能像鸟一样在平地或水面上起飞。
在这里插入图片描述
问题重述

考虑到问题陈述中确定的背景信息和限制条件，我们需要开发模型来评估两只翼龙的飞行能力，并解决以下问题：
•问题1：通过化石和相关研究恢复选定翼龙的身体数据，并建立模型预测翼龙正常飞行时的速度。
•问题2：根据前面问题中获得的模型和机翼颤振模型，计算所选翼龙的机翼襟翼频率。
•问题3：讨论选定的大型翼龙是如何起飞的，以及它们是否可以像鸟一样在平地或水面起飞。

模型的建立与求解整体论文缩略图

在这里插入图片描述

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：(代码和文档not free)

1 clear;
2 clc;
3 v=[]; v(1)=10; %initial speed
4 gma=[]; gma(1)=0;
5 psi=[]; psi(1)=0;
6 x=[]; x(1)=0;
7 y=[]; y(1)=0;
8 h=[]; h(1)=200;
9 z=[];
10 m=18.6; %quality
11 Sw=2.26; %wing area
12 b=5.34; %span
13 w10=25; %basic wind speed
14 CL=1.6;
15 rho=1.23;
16 phi=0;
17 L=[]; L(1)=abs(0.5*rho*CL*Sw*v(1)^2);
18 Cdpro=0.014;
19 g=9.8;
20 k=1.1;
21 pi=3.14;
22 hmin=0.03;
23 i=1;
24 t=0.01;
25 W=[];
26 D=[];
27 a=h(1);
28 while(L(i)<=(m*g*3) && h(i)>0.5 && a>0)
29 D(i) = 0.5*rho*Cdpro*Sw*v(i)^2+0.5*rho*0.01*Sw*v(i)^2+(rho*k*(CL*Sw*v(i
))^2)/(2*pi*b^2);
30 W(i) = w10*log(h(i)/hmin)/log(10/hmin);
31 x(i+1) = v(i)*cos(gma(i))*cos(psi(i))*t+x(i);
32 y(i+1) = (v(i)*cos(gma(i))*sin(psi(i))-W(i))*t+y(i);
33 z(i) = (-v(i)*sin(gma(i)))*t;
34 h(i+1) = h(i)+z(i);
35 gma(i+1) = (L(i)*cos(phi)-m*g*cos(gma(i))+m*(w10*(-v(i))*sin(gma(i))/(
log(10/hmin)*h(i)))*sin(gma(i))*sin(psi(i)))*t/(m*v(i))+gma(i);
36 psi(i+1) = (L(i)*sin(phi)+m*(w10*(-v(i))*sin(gma(i))/(log(10/hmin)*h(i)
))*cos(psi(i)))*t/(m*v(i)*cos(gma(i)))+psi(i);
37 v(i+1) = (-D(i)+m*g*sin(gma(i))+m*(w10*(-v(i))*sin(gma(i))/(log(10/hmin
)*h(i)))*cos(gma(i))*sin(psi(i)))*t/m+v(i);
38 L(i+1) = 0.5*rho*CL*Sw*v(i)^2;
39 i = i+1;
40 fprintf(’%d’, v(i));
41 a = h(i)-b/2*abs(sin(phi))*cos(gma(i));
42 end
43 T = i*t
44 scatter3(x, y, h, ’p’);
45 xlabel(’X’), ylabel(’Y’), zlabel(’Z’);
46 grid on;

1 clear;
2 clc
;
3 v=[]; 4 v(1)=15.8121;
%initial speed
5 rou=1.23; 6 Cdpro=0.014; 7 m=18.6; %quality 8 g=9.8; 9 Sw=2.26; %wing area
10 CL=1.6;
11 k=1.1;
12 pai=3.14;
13 b=5.34;
%span
14 hmin=0.03;
15 w10=25;
%basic wind speed
16 i=1;
17 t=0.01;
18 L=[];
19
20 W=[];
21 phi=0;
22 gma=[];
23 gma(1)=1.283;
24 psi=[];
25
%psi(1)=-0.0176;
26 psi=-0.0267;
27 x=[];
28 x(1)=22.5759;
29 y=[];
30 y(1)=-88.6534;
31 z=[];
32 h=[];
33 h(1)=189.752;
34 D=[];
35 a=h(1);
36
37 n=[-b/2:0.01:b/2];
38 ky=[];
39 sum=0;
40 for d=1:(b/(2
*0.01))+1
41 ky(d)=-0.2
*n(d)^2+0.2
*b^2;
42 sum=ky(d)+sum
;
43 end
44 lambda=2
*sum/(b/0.01+1);
45
46 t_pudong=0.5;
47 theta_Max=45;
48 theta0=0;
49 theta=[];
50 dao_theta=[];
51
52 beta_Max=15;
53 beta0=0;
54 beta=[];
55 dao_beta=[];
56 f=4.96;
57 omega=[];
58 v_xiangdui=[];
59 yingjiao=[];
60 z_yingjiao=[];
61 Fnc=[];
62 Fna=[];
63 Fdf=[];
64 bianhua=[];
65 bianhua(1)=0;
66
67 for i=1:f*t_pudong/0.01
68 for j=1:(b/0.01)+1
69 Fnc(i)=0;
70 Fna(i)=0;
71 Fdf(i)=0;
72 t=i*0.01;
73 theta(i,j)=45*cos(2*pai*f*t)+theta0;
74 dao_theta(i,j)=2*n(j)*pai*f*theta_Max*(2*pai*f*t);
75 beta(i,j)=n(j)*(beta_Max*cos(2*pai*f*t)+beta0)/b;
76 dao_beta(i,j)=-0.25*((-0.2)*n(j)^2+0.2*b^2)*(2*pai*f*abs(n(j))*
beta_Max*cos(2*pai*f*t))/b;
77 omega(i,j)=(2*(0+2+beta(i,j))*v(i))/(2+lambda);
78 v_xiangdui(i,j)=(((dao_theta(i,j)-dao_beta(i,j))*cos(2+beta(i,j))-
omega(i,j))^2+(v(i)+(dao_theta(i,j)-dao_beta(i,j))*sin(2+beta(i,j)))^2)
^(1/2);
79 yingjiao(i,j)=atan(((dao_theta(i,j)-dao_beta(i,j))*cos(2))/(v(i)+(
dao_theta(i,j)-dao_beta(i,j))*sin(2)));
80 z_yingjiao(i,j)=yingjiao(i,j)+theta(i,j)+2;
81 Nc(i,j)=0.5*rou*v_xiangdui(i,j)^2*2*pai*z_yingjiao(i,j)*(-0.2*n(j)
^2+0.2*b^2)*0.00001;
82 Fnc(i)=abs(Nc(i,j)*0.00001)+Fnc(i);
83 ax(i,j)=abs(n(j))*4*pai^2*f^2*45*cos(2*pai*f*t)+bianhua(i)*sin(2+
beta(i,j))+v(i)*cos(2+beta(i,j))+0.5*((-0.2)*n(j)^2+0.2*b^2)*(-n(j)*
beta_Max*4*pai^2*f^2*cos(2*pai*f*t)/b);
84 Na(i,j)=0.25*rou*pai*(-0.2*n(j)^2+0.2*b^2)^2*ax(i,j);
85 Fna(i)=Na(i,j)*0.001+Fna(i);
86 Df(i,j)=0.664*(z_yingjiao(i,j))*(v(i)*cos(2+beta(i,j))+omega(i,j)*
sin(2+beta(i,j)));
87 Fdf(i)=Df(i,j)*0.001+Fdf(i);
88 end
89 t=0.01;
90 W(i)=w10*log(h(i)/hmin)/log(10/hmin);
91 x(i+1)=v(i)*cos(gma(i))*cos(psi(i))*t+x(i);
92 y(i+1)=(v(i)*cos(gma(i))*sin(psi(i))-W(i))*t+y(i);
93 z(i)=(-v(i)*sin(gma(i)))*t;
94 h(i+1)=h(i)-z(i);
95
96 gma(i+1)=(Fnc(i)*cos(phi)-(m*g+Fna(i))*cos(gma(i))+(m+Fna(i)/9.8)*(w10
*(-v(i))*sin(gma(i))/(log(10/hmin)*h(i)))*sin(gma(i))*sin(psi(i)))*t/(m*
v(i))+gma(i);
97 psi(i+1)=(Fnc(i)*sin(phi)+(m+Fna(i)/9.8)*(w10*(-v(i))*sin(gma(i))/(log
(10/hmin)*h(i)))*cos(psi(i)))*t/(m*v(i)*cos(gma(i)))+psi(i);
98 v(i+1)=(Fdf(i)+(m*g+Fna(i))*sin(gma(i))+(m+Fna(i)/9.8)*(w10*(-v(i))*sin
(gma(i))/(log(10/hmin)*h(i)))*cos(gma(i))*sin(psi(i)))*t/m+v(i);
99 a=h(i)-b/2*abs(sin(phi))*cos(gma(i));
100 bianhua(i+1)=(v(i+1)-v(i))/0.01;
101 if(h(i)>200.5||a<0)
102 break;
103 end
104 end
105 scatter3(x,y,h,’p’);
106 xlabel(’X’),ylabel(’Y’),zlabel(’Z’);
107 grid on;

1 clc;
2 clear;
3 data=xlsread(’sensitivity analysis.xlsx’);
4 x=data(:,1);
5 y=data(:,2);
6 z=data(:,3);
7 scatter3(x,y,z)
8 figure
9 [X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))’,linspace(min(y),max(y)),’v4
’);
10 pcolor(X,Y,Z);shading interp
11 figure,contourf(X,Y,Z)
12 figure,surf(X,Y,Z);
13 figure,meshc(X,Y,Z)
14 view(0,0);
15 figure,meshc(X,Y,Z);