昆虫飞行振翅模拟.doc

北京大学政学者论文集（2001年） 昆虫振翅飞行原理的数值模拟研究昆虫振翅飞行原理的数值模拟研究The Numerical Simulation of Insects Flight北京大学 力学与工程科学系理论与应用力学专业 98级 程暮林摘要昆虫经过漫长的进化过程，获得了令人惊叹的飞行技巧，其高度的飞行机动性令现代飞行器黯然失色。随着非定常流理论的不断完善，人们逐渐认识到非定常效应在昆虫飞行中起着重要作用。前人采用实验和无粘假定下的理论分析的方法对此进行了一些研究。然而由于该问题本身的复杂性，使得理想化的理论分析和试验的研究手段难于得到全面和深入的结果。本文采用流体力学数值模拟方法研究本问题，通过采取转换参考系的方法对昆虫的振翅飞行所扰动而产生的流场进行了数值模拟，克服了直接对此问题进行数值模拟所必然遇到的困难“动边界问题”，缩短了计算时间。文中通过改变控制飞行的各个参数，从而比较详细的研究了这些参数对昆虫飞行的影响，为研究昆虫的飞行控制机制提供了依据。AbstractInsects have gained very surprising ability of flight after thousand years of evolution. People have come to realize the important effect of unsteady flow in getting enough lift for insects with the improvement of Unsteady Flow Theory. We have got some useful results with experimental and theoretical methods, but now it is difficult to go futhure with such methods. In this article , I use the numerical simulation method (the standard k-turbulence model) in a non-inertia reference linked to wings, thus avoid the difficulty from moving boundary and save the time of computation. Then I summarize the effect of different parameters controlling the flight of insects with this computing method.一 前言经过几千万年漫长的进化过程，昆虫获得了令人惊叹的飞行本领，有的昆虫能够悬停于空中，有的能够以很小的曲率半径急速转弯等等。这种高度灵活的飞行控制能力是人类现有的飞行机械所无法比拟的。因此，关于昆虫振翅飞行原理的深入研究将有助于提高现有的飞行器的性能，推动微型飞行器的开发和研制。昆虫学家们很早就对昆虫的飞行发生了浓厚的兴趣，随着技术的发展，他们采用高速摄影，摄像技术对各种昆虫的飞行进行了大量的拍摄，并且同时测得了一些与昆虫飞行相关的参数（如振翅频率，振翅的振幅等等）【14】，通过对这些数据进行数据拟合等处理后，他们发现这些参数间存在一定的关系，例如，翅的表面积正比于昆虫质量的2/3方，而振翅频率正比于昆虫质量的1/4方。虽然这些规律均是比较粗糙的，但能够证明这样一个观点：虽然经过上千万年的漫长进化过程，由于环境的差异导致各种昆虫彼此的形态，大小，习性均有较大的差异，然而由于受到“振翅飞行”这一物理条件的限制，各种昆虫的身体结构，翅的形态，翅的运动模式等飞行参数必然符合一定的力学规律。然而生物学家的研究受到一定的限制，未能对这个问题进行深入的数理分析，需要力学工作者的参与和合作。传统的空气动力学只注意翼型在定常气流中的升力，它不能解释昆虫为什么能够通过扑翅获得足够的升力而浮于空中：按传统空气动力学计算出来的升力远小于昆虫的重力【58】，因此产生了“大黄蜂飞不起来”的谬误。随着非定常流理论的完善，人们逐渐认识到非定常效应在昆虫飞行升力的产生中起重要作用。T.Weis-Fogh的研究解开“昆虫升力之谜”，在研究了大量的昆虫飞行的资料以后，提出了振翅飞行的ClapFling机制【9】。所谓的ClapFling指的是：昆虫的一对翅上抬至最高点，两翅面保持对应翅脉平行的合在一起（即Clap），然后从翅前缘（the leading-eadge）绕翅后缘(the trailing-eadge)快速的剥开。当两翅面张开到一定程度时，两翅面彻底分开，分别继续作平动。出于理论分析和实验上的考虑，这种运动引起的流动被简化二维流动。右图是这种运动的示意图（fig.1）。fig.1 开始阶段两翅面靠在一起（fig.1.a），然后两翅面以某一角速度绕后缘匀速转动（fig.1.b），翅间的张角增大，空气挤入两翅面所张开的区域中。当两翅面张开的角度增大到一定值时（通常认为是/3）两翅面脱离接触(fig.1.c)，各自以常速度平动(fig.1.d)。（取自T.Weis-Fogh的论文）在此之后，M.J.Lighthill【10】, R.H.Edwards和H.K.Cheng【11】采用无粘位势流的假设从理论上应用复变函数论的方法，求得了此问题的理论解。实验方面，T.Maxworthy【12】和G.R.Spedding【13】分别完成了相似的二维流动显示实验，其结果大致上验证了理论上的解析解，同时也显示了由于昆虫飞行时的平均雷诺数（Re）处于1002000这一范围内，所以空气的粘性作用是不可忽略的。理论分析中的无粘假定导致了实验与理论之间的差别。近几十年来，由于流体实验技术的发展，使得对由于昆虫振翅飞行所扰动而产生的流场的观察和测量更为细致。C.P.Ellington【14】采用活的吊飞昆虫进行流场显示实验，其结果表明：昆虫向下拍翅过程中，在翅的前缘将产生一个涡，被称为前缘涡（the leading eadge vortex）。由于这个涡的存在，导致翅面的上方产生一个低压区，因而昆虫获得升力。由于在实验中采取活的昆虫作为实验对象使测量受到很大的限制，例如无法测定翅根部的受力情况等等。人们开始制造微型振翅机械来代替活的昆虫作为实验对象。这样做法的另一个好处是可以模拟空中自由飞行的昆虫的振翅模式，从而克服测量吊飞昆虫中所无法避免的误差。M.H.Dickinson等人按照这一思路完成了试验【15】，取得了很大的成功（fig.2）。fig.2 升力和阻力随时间周期性变化的曲线。红线是实验测量值，蓝线是翅平动产生的升力和阻力（采用机翼理论计算）。（取自M.H.Dickison的论文）由于近几年来计算流体力学的飞速发展，数值模拟逐渐成为流体力学研究的另一重要手段。Z.JangWang通过利用涡流法模拟昆虫的振翅飞行的二维流场【16】发现：由二维效应产生的升力(fig.3)足够使昆虫浮于空中。fig.3 采用涡流法模拟昆虫振翅飞行所得到的升力和阻力随时间变化的曲线图。可以看出与前面提到的模型实验结果十分相似。（取自Z.JangWang的论文）研究昆虫飞行的根本目的是建立一套完整的振翅飞行的理论，掌握昆虫飞行和控制的机理，并便指导微型飞行器的研究和开发。以前的研究多集中于讨论振翅飞行的机理性问题，如昆虫是如何通过振动翅膀而获得足够的升力，但对于控制飞行的各个参数（如振翅频率，振翅的振幅等等）的变化对升力和阻力的影响进而控制飞行的方向，飞行的速度，飞行的姿态这一问题讨论得较少。实验上，设计实验的难度很大，并且试验要求较高的测量精度和实验技巧，这就极大的限制了实验手段在这一问题上的应用。数值模拟具有实验手段无法比拟的优点，因而比较适用于这一问题。然而如果直接进行数值模拟，将遇到“动边界问题”，这在计算流体力学中至今仍是一个棘手的问题，即使可以计算，计算的时间也将比较长，对计算机的软硬件的要求也会比较高。本文通过转换坐标系的方法，克服了这些问题，大大缩短了计算周期，使得仔细研究控制飞行的各个参数对飞行的影响成为可能。文中通过改变控制飞行的各个参数：振幅，平动和转动间的相位差等等，较为详细的研究了各种情况下升力和阻力的变化，而且发现不同的翅尖轨迹可以造成不同的升力和阻力等等一些现象。二基本假设和简化的振翅模型昆虫是一个十分复杂的生物系统，为了研究的方便将其进行了一些简化，从而可以采用数值模拟的方法进行分析：(1). 昆虫振翅飞行所扰动的流场本质上是一个三维流场，翼尖涡等三维效应是存在的。然而由于昆虫的翅本身的特征长度较小，这样三维效应有可能并不十分显著，另外近似的分析：取通过翅纵向方向中点的横截面，在这一截面两侧距此截面距离相等的翅表面上所受到的升力，阻力的平均值将与该中截面近似相等。因此，本文中仅研究二维的昆虫振翅飞行。(2). 虽然昆虫的躯体不具备良好的流线型，但由于其特征尺度与其到我们所研究的翅的中点横截面的距离相比较小，所以在本文中忽略虫体对流场的干扰作用。(3). 不同种类的昆虫振翅频率相差较大，本文中以振翅频率处于0300Hz的昆虫为主要研究对象。由此可以计算出相应的流场的特征速度，对应于振翅频率为300Hz的情形，根据已有的生物学资料可估计出特征速度远小于100m/s因此，本文采用空气的不可压缩模型进行计算。(4). 生物学已有的资料表明：真实昆虫的翅表面分布有刚毛，鳞毛等毛状器官，这些结构的存在将破坏翅表面的边界层，但是在本文中认为翅表面是光滑的。此外，由于翅的柔性使得翅在拍动过程中受到空气的作用后将发生变形（扭曲成螺旋状），不过在本文中这也被忽略掉了。综合以上假设，我们进行数值模拟的假定是：昆虫的翅被简化为光滑刚性平板，周围的空气采取不可压缩模型。三虫振翅模型的数学描述如下图所示，fig.4(a)是真实昆虫振翅飞行时的翅尖轨迹曲线（空间轨迹曲线和本体轨迹曲线）。由图上可以看出，翅尖的本体轨迹基本上是一个很扁平的曲线，这条曲线的长轴方向所在平面被称为平均拍翅平面。Fig.4(b)是简化的拍翅过程示意图。Fig.4 昆虫振翅飞行模型的数学描述。Fig.4(a) 通过高速摄像取得的昆虫振翅时翅尖的轨迹。图中所画的两条轨迹，一条是相对于地面的翅尖的空间轨迹，另一条是相对于昆虫自身的本体轨迹，轨迹上的小短线代表翅面。可以看出翅尖的空间轨迹相似于正弦曲线。 Fig.4(b) 简化的振翅模型图。XOY是固连于地面的参考系，xoy是固连于翅面上的坐标系，msn是相对于XOY静止的坐标系，在此坐标系下用（，）描述点的坐标。图中实短线代表翅下拍过程中翅面的位置，虚短线代表翅上抬过程中翅面的位置。翅的运动可以分解为平动和转动，翅的平动用翅尖的中心坐标（，）来描述，采用攻角来表征翅的转动。综合分析fig.4(a)并发展文献【16】中的模型，给出如下方程描述翅的运动：翅尖的运动方程： （1） （2）攻角变化规律： （3）其中，A0 为振幅，B0为横向偏移量，为振翅频率，为横向振翅频率，为启动时的相位，为平动和转动间的相位差。微分（1）（3）可以得到在msn坐标系下的速度和加速度： （4） （5） （6） （7）在msn坐标系下的角速度和加速度： （8） （9）由于坐标架msn与坐标架XOY相对静止，所以上式也是在参考系XOY下的昆虫的翅的角速度和角加速度的表达式。为了获得在参考系XOY下昆虫翅的速度和加速度，利用以下坐标变换公式： （10） （11）其中，是平均拍翅平面的倾角。 是常数。得到相对于参考系XOY 的速度： （12） （13）相对于参考系XOY 的加速度： （14） （15）四振翅飞行的昆虫周围流场所满足的方程首先，在固连于地面上的参考系XOY下考虑问题，流体运动满足连续性方程（不可压缩条件）： （16）其中，是速度场流体运动满足Navier-Stokes方程： （17）其中，是空气密度，由不可压缩性条件，是常数。是压力场，是动力学粘滞系数初始时，假定昆虫静止，在t=0时刻开始起飞，因此流场满足初条件： （18）其中，是无穷远处来流速度，如果昆虫原地悬停或处于吊飞状态，则此速度为零；如果昆虫处于自由飞行的状态，此速度与昆虫的飞行速度相反。 是坐标系XOY的单位向量在翅的表面上，满足无滑移壁面条件： （19）其中，C为翅面的边界曲线， 为翅面上边界点的瞬时速度。在无穷远处，流场满足： （20）方程（16）（20）即是在固连于地面上的惯性参考系XOY中此问题所满足的流体运动方程。原则上，只要求解该方程组，即可解决此问题。然而，仔细观察上述方程组可以发现，内边界（即翅面的边界）是随着翅的拍动而不断变化的，这样如果采取上述方程组直接进行数值模拟将十分困难。此外，为了要达到研究不同的控制飞行参数对于昆虫飞行的作用的目的，将导致计算所用的程序十分复杂，有可能计算时间也比较长。我们采取转移参考系的办法，来克服此困难。将参考系由固连于地面上的惯性参考系XOY 换成固连于昆虫翅上的非惯性参考系xoy, 这样动边界问题消失了。原来的动边界“变成”非定常的边界条件出现在方程的边条件中。采用以下变换： （21）（22）其中，是流体质点在非惯性参考系xoy下的速度和加速度，是翅尖在坐标系XOY下的速度和加速度（按照前面的假定翅是刚性不变形的，因此，翅尖的速度和加速度是已知的）。，是分别是速度和角速度。是矢径。代入上述方程组（16）（20），得到非惯性参考系xoy下流场所满足的方程组： （23） （24） （25） （26） （27）其中，惯性力项的表达式为： （28）方程（23）（28）就是非惯性参考系xoy下流场所满足的方程。此方程组的分量形式为： 控制方程： （29） （30） （31）初条件： (32) (33)其中，翅表面的边条件： （34） （35）无穷远处所满足的边界条件： （36） （37）其中，Fx，Fy的形式为： （38） （39）实际进行数值模拟时，采用方程组（29）（39）。五网格生成和数值求解方法1 网格的生成计算中假定翅为椭圆形长轴长为c1cm，短轴b=c/10=0.1cm。计算区域为15倍于弦长c的圆形区域，计算区域的特征尺度远远大于翅的弦长，保证了在外边界上近似地满足上述无穷远边界条件（36），（37）。网格采用结构化网格，共有12960个结点，12800个控制体，最大控制体积1.004903E-4m3，最小控制体积为1.302605E-9m3。在翅表面，网格的特征尺度很小，小于翅表面附近的湍流边界层的厚度。考虑到翅尖附近流场变化较为剧烈，采用Gambit生成的网格如下图所示： Fig.5 计算时所使用的网格。左图是整个计算区域，可以看出相对于整个流场来说昆虫的翅是很小的。因此，翅振动引起的扰动将由于空气的粘性的耗散作用，而不能够传到外边界，外边界上的无穷远边界条件（36），（37）可以得到满足。右图是翅尖附近加密的网格。 2数值求解方法昆虫的翅是往复的振动，不但平动，而且同时转动，因此，可以认为其扰动所形成的湍流是充分发展的湍流。采用标准的模型进行计算。将NS进行平均化处理得到雷诺应力方程： （40） （41）为了使上述方程组封闭，采用Boussinesq 假设： （42）其中，k为湍流动能是湍流粘滞系数，可以通过湍流动能,湍流耗散率计算： （43）湍流动能,湍流耗散率满足的传输方程： （44） （45）方程（40）（45）就组成了封闭的方程组，即是标准湍流模式，其中的常数经过实验确定为： （46）实际计算中，采用Fluent完成计算。计算中发现大约经过1015个振翅周期升力和阻力开始呈稳定的周期性变化，分析中着重于这一阶段。六计算结果与讨论在早期的研究中，没有考虑到非定常效应对昆虫获得足够的升力所起到的作用，因而产生了“大黄蜂飞不起来”的谬论。为了确切的了解在昆虫振翅飞行中非定常效应的作用，我首先采取定常的计算程序对昆虫的翅的空气动力学性能进行计算（主要是升力和阻力系数）。另一方面，也可通过此计算了解计算程序的适用性和稳定性，以及所采取的计算方法是否可行。通过计算不同攻角下翅面上（fig.6）的升力和阻力，并且采用以下公式计算得到升力和阻力系数： （47） （48） 其中，V是来流速度，S是特征面积。数值模拟的结果如下图所示（fig.7）：（a）(b)Fig.7 不同攻角下翅上的升力和阻力系数. fig.6(a)是采用上述方法计算出来的结果；fig.6(b)是引自： M.H.Dickinson 的论文中的实验结果。比较两者，数值计算与实验结果十分相似，但计算结果稍小于实验结果，可能是由于试验和计算所采用的昆虫翅的模型不是十分一致。比较试验和计算的结果可以发现：数值计算的升力和阻力系数的变化趋势与试验所得到的结果完全一致，但稍小于实验的结果。这可能是由于试验和数值计算所采用的昆虫翅的几何外形不是十分一致。 以上定常的昆虫翅的空气动力学性能的计算结果说明Fluent的程序有可能满足我们的需要。首先采用以下参数我们进行了计算，振翅的振幅和横向振幅分别为 ， 角度 ， ， 振翅的频率为 ，空气的运动学粘滞系数 空气的密度 计算大约1015个周期后，升力和阻力稳定的周期性变化，如下图所示（fig.8）：fig.8升力和阻力随时间变化的曲线。从开始计算大约经过1015个周期升力和阻力开始随时间稳定的变化。图中横坐标表示周期数，纵坐标表示单位面积上作用的升力和阻力，因此，单位为N/m.图上显示的结果和Z.Jane Wang 十分相似，表明计算程序的结果是可信的。正如前文提到的，计算是在固连于昆虫翅上的参考系下进行的，所以不需要变动网格，仅仅边界条件是随时间变化的。然而在非惯性参考系下许多物理量，以及他们之间的关系是不同于惯性系下的关系，但是可以验证在惯性参考系XOY下的涡量和在非惯性参考系xoy下的涡量仅仅相差一个翅的转动角速度涡量，也就是： （49）计算所得的涡场随时间变化的过程如下图（fig.9）所示。虽然不同时间的涡场图之间没有可比性，因为不同时刻翅的转动角速度不同，然而从图上可以看出同一时刻涡量的分布。由图上可以明显看出，在平动速度最大时，即t=0 T，t=0.5 T时，翅两端的涡量方向相反；在转动角速度最大时，即t=0.25 T，t=0.75 T时，翅两端的涡量方向相同。转动速度最大时平动速度最大时Fig.9昆虫振翅飞行时周围空气的涡场随时间变化的过程（在非惯性参考系xoy下观察）。惯性参考系XOY下的涡量与非惯性参考系下的涡量仅相差一个常数：翅的转动角速度，而且这一关系在任一时刻对全流场的每一空间点都成立。左一列的两张图表示转动的角速度最大时的涡场图，在翅的两端附近涡量具有相反的方向；右一列的两张图表示平动速度最大时的涡场图，在翅的两端附近涡量具有相同的方向正如在前言中提到的，在研究昆虫飞行的早期，人们并未认识到非定常效应对昆虫获得升力的作用，为了确切分析非定常效应这一作用我们定义如下两个量： 定常的升力： （50） 定常的阻力 （51） 其中，是升力和阻力系数，可以采用上一节计算出来的值。是平动速度，由于翅的运动已知这一值可由（12）（13）式计算而得。下图表示这样处理的结果（fig.10）。fig.10 一个多周期内升力随时间变化的曲线图。图中横坐标轴代表周期数，纵坐标轴代表单位截面上升力的大小。其中，红线是数值计算的结果（定常效应非定常效应），蓝线是按上述定义计算的定常效应产生的升力，黄线是两者之差，即为非定常效应。非定常效应对阻力也有较大的影响。 从图中可以明显看出，非定常效应使升力的峰值明显增大，通过实际的计算发现，非定常效应还使得平均升力增大20.07，这样大的差别足以解释为什么采用定常流理论无法解释昆虫升力之谜。类似于升力，昆虫翅上所受到的阻力也受到非定常效应的影响，虽然相对较小一些（8.21%）。 在实际设计微型飞行器的工程中，关心的另外一个问题就是振翅飞行过程中所消耗的功率。由于昆虫的翅的质量相对于昆虫的躯体的质量来说是一个小量，所以以下在忽略昆虫翅的质量的假设下，讨论飞行中功率的分配。昆虫飞行过程中消耗的功率由三部分组成：产生升力消耗的功率，产生推力消耗的功率，克服扭矩转动而消耗的功率。以上三部分功率可以表达成如下的形式： （52） （53） （54） 其中，分别是翅平动速度的X轴，Y轴分量，转动角速度。这样，昆虫飞行时消耗的总工率就可以表达为三者之和： （55）按公式（51）（54），（55）计算三部分的即时功率，所得功率随时间变化图（fig.11）如下：fig.11 飞行过程中消耗的即时功率随时间变化的过程。图中，横坐标代表时间单位是妙；纵坐标代表翅的单位截面上消耗的功率。红线代表升力消耗的功率，蓝线代表推力消耗的功率，黄线代表翅转动所消耗的功率，浅红色的曲线代表飞行过程中消耗的总工率。由图上可以清楚看到前半个振翅周期（即翅由最高点向下运动的过程）主要是升力在消耗功率，而推力消耗的功率很小；在后半个振翅周期（即翅由最低点向上运动的过程）主要是推力在消耗功率，而升力所消耗的功率很少。翅的转动所消耗的功率在整个过程中都很小。对这些曲线在一个周期内积分，可以得到平均功率： （56） （57） （58） （59） 在振翅过程中，推力和升力消耗的功率基本上分别消耗总功率的50%，而转动消耗的功率实际上可以忽略。对于飞行器设计来说这样的功率分配是十分不合理的。因此，从这个意义上来说振翅飞行并不是一种有效的飞行方式。正如在引言中提到的，计算是在固连于翅面上的参考系中进行的，所以不需要在计算过程中更改网格，动边界条件转化为随时间变化的边界条件，而且保持边界不动。这样对于不同的控制飞行的参数，仅仅需要改变若干参数即可实现对于不同条件的模拟。利用这一方法研究了若干控制飞行的参数。昆虫的振翅过程可以看作翅的平动和转动的叠加，因为平动最大速度发生于翅面处于水平位置，而对于的情形最大转动角速度发生于上拍的最高点和下拍的最低点，这样振幅其实控制了平动和转动之间的空间距离。此外，显然另一控制平动和转动的参数就是平动和转动间的相位差，它控制了平动和转动之间的时间距离。增大振幅，升力曲线的最大值显著提高，而在其余的位置上变化不大。利用下式（59）对平均升力进行无量纲化， （60）这样就得到了无量纲化的平均升力于振幅之间的关系（fig.12）：fig.12 无量纲化的平均升力与振幅之间的关系。平均升力采用（59）式进行无量纲化，横轴表示振幅，单位是cm。从图上明显可以看出当振幅增大到一定程度，无量纲化的平均升力趋于一个基本固定的常数（0.55左右）。 无量纲化的平均升力随着的振幅增大而增大，当振幅增大到一定值时，无量纲化的平均速度趋于一近似恒定的值（0.55）。这是因为，振幅增大使得转动和平动的最大值发生的位置的空间距离加大，两者所扰动的流场相互的干扰减小。另一方面，平动和转动间的相位差也表现出相似的作用。下图（fig.13）表示不同的相位差所对应的升力曲线（随时间变化）。ffig.13 不同平动和转动间的相位差下升力随时间变化的曲线，其中，相邻的曲线间相差。图中，时间轴是横轴，纵轴表示单位长度的翅的截面上所受到的升力。明显可以看到，除了在几个曲线的极大值和极小值的位置上，升力随相位差的不同而发生变化外，在其他的时间上各条曲线基本上是重合在一起的。这一现象是由于相位差改变了最大转动速度发生的位置（最大平动速度发生的位置不随相位差变化，参看（1），（2），（3）式）,通过平动和转动所扰动的流场的作用，升力和阻力发生变化。通过积分一个周期内的升力值，得到平均升力随平动与转动间相位差变化的曲线图（fig.14）fig.14 平均升力虽平动和转动间相位差变化的曲线图。图中，横轴表示相位差，纵轴表示单位长度上的升力。图中，红色的叉子表示计算所得到的平均升力，相邻两个计算值得相位差相差。当相位差减小时，平均升力也随之减小；反之平均升力随相位差的增大而增大。图上每一个红色的叉子代表一个计算值，相邻两点之间相位差相差。相位差增大时，平均升力也随之增大，反之，平均升力减小。相位差发生较小的变化就可以引起升力和阻力的较大变化，这一现象在试验中已经被发现，并被称为“升力和阻力对相位差的敏感性”。这一现象的发现为解释昆虫如何在飞行中十分灵活的控制飞行的方向，飞行的姿态提供了依据。一种简单的假设是：昆虫在飞行过程中可以通过控制左右两个翅的振动的平动和转动的相位差，使左右两个相位差不同，从而使作用于左右两个翅上的升力和阻力不同，这样就可以形成偏航力矩，俯仰力矩等等。由于控制这种相位差相对于控制振翅的振幅来说，需要的消耗能量少的多，而且敏感性比较强，所以在微型飞行体的设计中是一个可以考虑的控制飞行的方案。在第二节的“基本假设”中，昆虫的翅被简化为具有光滑表面的平板，而忽略了翅表面上的各种微观结构对周围流场可能的影响。可以通过另外一个途径，近似地了解这种微观结构的作用。由于微观结构主要是影响翅面附近流场的边界层，并且主要是通过于流体的粘性摩擦作用耗散流场的动能，而在边界层中粘性力和惯性力相当，粘性力与流体的粘滞系数成正比，因此微观结构对边界层的这种作用可以通过改变流体的粘滞系数加以模拟。流体具有较小的粘滞系数就相当于昆虫的翅表面比较光滑；反之，具有较大粘滞系数的流体的流场代表昆虫的翅表面各种微观结构比较多。利用数值计算的方法可以比较容易的实现各种情况，下图（fig.15）就是计算的结果：fig.15 平均升力随空气的动力学粘滞系数变化的曲线，红色的点是计算的平均升力的值，相邻的点对应的空气的动力学粘滞系数相差一个量级。值得注意的是，横坐标是空气的动力学粘滞系数，此轴为对数轴。在空气的动力学粘滞系数处于一定范围内时，平均升力变化很小。由图上可以看出，当空气的动力学粘性系数在正常空气的动力学粘滞系数左右变化时平均升力变化不大。当空气的动力学粘滞系数增大到的量级上时，平均升力急剧下降，说明如果昆虫翅上的微观结构的作用还是比较大的。前面的讨论都是在假定昆虫是静止不动（相当于昆虫处于吊飞状态）的前提下进行的。这种假定情形与处于飞行的昆虫的情形有一个十分明显的力学上的差别：处于水平飞行的昆虫，从固连于昆虫上的参考系来看，其无穷远处的边界条件上来流速度不为零。考虑到这种差别以后，翅上所获得的升力和阻力有相应的变化，实际计算的结果如下图所示（fig.16）：（假定昆虫处于水平匀速飞行状态，无穷远来流速度与飞行速度相反，参看（20），（36），（37）式）fig.16 平均升力随速度（用2m/s无量纲化）变化的曲线。除了上面提到的无穷远边界条件不同外，吊飞昆虫和自由飞行的昆虫在其他方面也存在着差异。根据已有的生物学资料：处于吊飞状态下的昆虫在振翅过程中其翅尖轨迹是近似于椭圆的一条十分扁平的曲线；而处于自由飞行的昆虫的翅尖轨迹是“八字形”或“双八字形”；或是更为复杂的曲线，如图所示（fig.17）。fig.17 昆虫翅尖轨迹图。(a)椭圆形轨迹代表昆虫处于吊飞状态下，(b),(c)“八字形”轨迹和“双八字形轨迹”对应于处于自由飞行的昆虫振翅模式。轨迹上的小箭头表示翅尖的运动方向。以上的讨论都是假定翅尖沿某一倾斜的直线往复运动，即横向振幅。关于这一问题的研究似乎未引起太多的注意，但计算的结果显示翅尖轨迹对昆虫振翅所获得升力和阻力有较大的影响。在前面提到的昆虫振翅的数学模型（1），（2），（3）中，取不同的的值就可以实现不同的翅尖轨迹。具体计算中，取该比值为1，2，3分别模拟了椭圆形轨迹，八字形轨迹，双八字形轨迹，横向振幅，。计算结果如下图所示（fig.18）: fig.18(a) 椭圆形翅尖轨迹时，升力随时间变化的过程。fig.18(b) “八字形”翅尖轨迹时，升力随时间变化的过程。fig.18(c) “双八字形”翅尖轨迹时，升力随时间变化的过程。Fig.8 不同的轨迹下，振翅飞行的昆虫所获得的瞬时升力随时间变化的过程。fig.18(a) 椭圆形翅尖轨迹时，升力随时间变化的过程; fig.18(b) “八字形”翅尖轨迹时，升力随时间变化的过程。fig.18(c) “双八字形”翅尖轨迹时，升力随时间变化的过程。横轴表示时间，纵轴表示瞬时升力。明显可以看出，在(a)中三条曲线基本上重合在一起，而在(b),(c)中在极大值的位置上三条曲线明显分开，尤其是在(c)中这一现象更加明显。由图上可以看到，当翅尖轨迹是椭圆形时，增大横向振幅，升力曲线基本上没有变化。但当翅尖轨迹是“八字形”时，增大横向振幅，升力曲线在原极大值的位置上继续增大，当翅尖轨迹是“双八字形”，这一效应更加明显了。同样地，对瞬时的升力和阻力在一个振翅周期内积分，可以得到平均的升力和阻力，这一计算结果如下表(tab.1)所示：tab.1 不同翅尖轨迹对应的平均升力和平均推力值，箭头表示横向增大时，平均推力或平均升力变化的方向。虽然目前为止还找不到一个合适的方式解释这一问题，但可以提供一个可能解释此问题的想法：仔细观察不同的翅尖轨迹，可以发现在椭圆形轨迹，“八字形”轨迹，“双八字形”轨迹，三者中椭圆形轨迹所围的区域最大，“双八字形”轨迹所围的区域（三个小块区域）最小。因此，椭圆形轨迹对流体的“限制作用”最小，“八字形轨迹”其次，“双八字形”轨迹对流体的“限制作用”最强，因为它围成的区域最小。“限制作用”越强所获得的升力的峰值也就越大，但同时升力变化得也就越剧烈。这一解释方式明显是不十分严谨的，需要进一步的研究。然而这一现象对于微型飞行器的设计却有比较重要的意义：一方面，设计中可以利用这一方式提高可以获得的升力大小，还可以作为一种产生控制力，控制力矩（俯仰力矩，偏航力矩等等）的方案。但另一方面，必须防止由这一现象产生的升力或推力剧烈变化，导致飞行器丧失飞行的稳定性。利用这一振翅模型和计算方法，我们还计算了起飞段的情形。简单的，起飞段的模型可以从（1），（2），（3）式稍加变化的来。考虑下式： （61）其中，是稳定时的振幅，n控制达到稳定所需要的时间，T为振翅周期。将（60）式带入（1），（2），（3）式即可得到起飞段昆虫振翅过程的数学描述。采用同样的计算方法可以得到升力随时间的变化如下图所示：fig.19 起飞段升力随时间的变化过程。振翅模式按照（61），（1），（2），（3）决定。七结论以上的计算结果表明转换参考系处理这一问题的方法是可行的，所得到的结果与已有的试验的结果可以比较好的符合。虽然从理论上可以有许多参数（如振幅，频率等等等，频率，最大攻角的影响比较复杂，但我们也做了相关的数值计算，但需要进行进一步的研究。）可以控制飞行得速度，方向，飞行的姿态等，但是对于生物体昆虫来说，控制某些参数（翅的几何外形等）是不可能的，或者控制某些参数不是“十分经济的”。这样根据以上的计算结果，可以确定一些参数可能被昆虫所使用，来控制飞行。例如平动和转动间的相位差，以及翅尖轨迹等等。最后，值得一提的是，三维情形下的振翅的数值模拟相对于二维有更多的困难。转换参考系的方法原则上可以同样地应用于三维情形，但是此时就必须考虑昆虫的躯体对流场的影响，因此可能需要其他的方法加以克服。致谢首先，我要感谢李政道和他夫人提供这一基金，使得我有机会与处于科学研究第一线的学者接触，直接聆听他们的教诲，学习科学研究的方法。本文是在陈滨教授的指导下完成的，在论文工作的各个阶段，陈滨教授总是给我很多的启示。正是在陈教授的指导，鼓励，帮助下，这篇论文才能得以完成。在我完成“政”学者计划的两年期间，陈教授不仅在专业上引导我进入学术研究的殿堂，教给我科学的研究方法，而且更重要的是交给我做人的道理。陈滨教授的渊博的学术知识，敏锐的抓住问题核心的能力，和对力学发展方向的洞察力都是我学习的榜样。陈滨教授两年间对我的学习，生活上的指导，帮助将使我终身受益。在我的后一年的研究时间里，我到陈耀松教授的实验室，完成数值计算的工作。在此期间，陈耀松教授不但无私的提供给我必要的设备外，还热忱地对我的科研工作给于了具体的指导和帮助。对我的具体的研究方向和论文的撰写进行了细致的关心。陈耀松教授渊博的学术知识，敏锐的抓住问题核心的能力，简单清晰的物理建立模型的能力，“活到老学到老”的学习精神和对力学发展方向的洞察力都是我学习的榜样。陈耀松教授对力学无尽的痴迷之情，对事业的执着追求，对他人无私的帮助，以及堂堂正正的生活态度给我以深刻的印象，值得我永远回味和学习。系里很多教授同样给我以热忱的帮助。与魏庆鼎教授，陆来教授，翟大熙教授，陈十一教授，陈国谦教授，吴介之教授，李植老师，谭文长教授，苏卫东教授，蔡庆东教授的讨论使我受益匪浅。特别是，李植老师一直关心我的研究进展，并且提供了许多好的建议。他还对我的文章进行了十分认真的修改，达到字斟句酌的程度，他的认真的态度令我感动。陈国谦教授，陈十一教授，吴介之教授在我的研究过程中一直关心我的研究进展情况，他们一直积极的引导我将研究推向更加深入的境界。同时，实验室里浓厚的学术气氛，轻松开放的学术氛围，都是我得以完成这一论文的必要条件。在这