基于风洞实验的蜜蜂飞行机理研究.pdf 179页

基于风洞实验的蜜蜂飞行机理研究赵红燕2015年1月 中图分类号：O3UDC分类号：531基于风洞实验的蜜蜂飞行机理研究作者姓名赵红燕学院名称机电学院指导教师马巍研究员答辩委员会主席秦太验教授申请学位工学博士研究方向计算力学与工程仿真学科专业力学学位授予单位北京理工大学论文答辩日期2015年1月 ReserchonFlightMechanismsofHoneybeesBasedonWindTunnelExperimentsCandidateName：HongyanZhaoSchoolorDepartment:MechatronicalEngineeringFacultyMentor:Prof.WeiMaChair,ThesisCommittee：Prof.TaiyanQinDegreeApplied:DoctorofEngineeringMajor：MechanicsResearchDirection：ComputationalMchanicsandEngineeringSimulationDegreeby:BeijingInstituteofTechnologyTheDateofDefence：Jan，2015 研究成果声明本人郑重声明：所提交的学位论文是我本人在指导教师的指导下进行的研究工作获得的研究成果。尽我所知，文中除特别标注和致谢的地方外，学位论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得北京理工大学或其它教育机构的学位或证书所使用过的材料。与我一同工作的合作者对此研究工作所做的任何贡献均已在学位论文中作了明确的说明并表示了谢意。特此申明。签名：日期： 北京理工大学博士学位论文摘要微型扑翼飞行器具有体积小、功耗低和机动灵活性高等一系列优点，能在复杂地形和恶劣环境下完成普通飞行器无法完成的任务，可广泛应用于军事侦察、航空摄影和环境监测等领域，因此对扑翼飞行机制的研究有重要的理论价值和应用前景。本文借鉴自然界中的飞行生物研究扑翼飞行的机理，选取蜜蜂为研究对象，采用风洞实验和数值计算相结合的方法，对蜜蜂的飞行机理进行了系统的理论研究。论文的主要研究工作如下：（1）设计并制作了用于小型昆虫流场显示实验的低速风洞。为保证风洞的流场品质，研究了风洞实验段、扩散段和收缩段的尺寸对流场品质规律的影响，确定了风洞各段的最佳尺寸，选择了维多辛斯基曲线作为收缩曲线线型；为降低湍流度，保证气流的平稳性，在稳定段和动力段添加了蜂窝器和阻尼网；设计并调试了发烟装置，对高速摄像机和灯光进行了合理的实验布局，开发了完整的流场显示实验系统。（2）在低速烟风洞中，采用两台高速摄像机同步拍摄蜜蜂，获得了蜜蜂的振翅规律；利用烟线法研究了蜜蜂振翅时的远场和近场流场的变化特征，研究了尾迹涡对诱导速度的影响，提出了蜜蜂振翅时涡的演变规律。翅下拍时，连续的翅尖涡尾迹形成涡管；翅外旋时，产生下拍停止涡与上拍启动涡；翅的弦向变形使得前、后翅尾缘处的涡量方向相反；下拍末期，前缘涡最明显，由翅根至翅尖方向前缘涡的尺寸是递增的；在蜜蜂振翅的下拍末期，左、右翅上各有一涡环，涡环由沿着翅展的前缘涡、翅尖涡尾迹涡管和翅根涡尾迹组成；根据涡环计算的升力大于蜜蜂的重力，这表明蜜蜂振翅产生的升力，不仅能够维持自身重力，还可携带一定的重物进行飞行。（3）建立了蜜蜂模型和坐标系，基于蜜蜂的振翅规律给出了翅的运动方程，利用动网格技术数值模拟计算了蜜蜂振翅时的流场。数值模拟结果验证了风洞实验的结论：翅在拍动时，翅上附着有前缘涡、翅尖涡、翅根涡和尾缘涡；翅下拍末期，由翅根至翅尖，前缘涡的尺寸增大，与翅尖涡连接在一起；下拍末期，尾缘处存在下拍停止涡，会降低升力，上拍初期，尾缘处存在上拍启动涡，能够提高升力；下拍末期，在左（右）翅上，翅尖涡尾迹、翅根涡尾迹和翅面附近的涡量能够组成完整的涡环，且涡环为圆形。此外，数值结果还表明：翅的旋转使得前缘涡在上（下）拍末期尺寸达到最大；在上拍末期，上拍期间形成的翅尖涡尾迹、翅根涡尾迹和翅边界处的涡量连接在一起形成耳状涡环；在上拍涡环和下拍涡环的中心，均有垂直于涡环平面向下的射流。I 北京理工大学博士学位论文（4）基于实体蜜蜂翅，建立了蜜蜂翅模型及翅的坐标系，定义了翅的升力、推力、气动功率和惯性功率；数值模拟计算对比了单翅模型、双翅模型和整体蜜蜂模型的力和功率，结果表明左、右翅的相互作用以及蜜蜂身体均对单个翅的气动影响可忽略不计；针对蜜蜂单翅模型，通过数值模拟计算分析了频率、最大拍动角、中期攻角、翻转时刻和旋转轴位置等不同影响因素对气动力和功率的变化规律。研究结果为将来研制微型仿生扑翼器有一定的指导意义。关键词：蜜蜂；风洞；涡；流场显示实验；数值模拟II 北京理工大学博士学位论文AbstractTheMicroAerialVehicle(MAV),whichcanfulfiltasksthatthenormalaircraftcannotperforminthecomplextopographyandbadenvironment,hasaseriesofadvantages,suchassmallervolume,lowerpowerdissipationandhighermaneuverability,andcanbewidelyappliedtomilitaryreconnaissance,aerialphotographandenvironmentalmonitoring.Consequently,flappingflighthasimportanttheoreticalvalueandapplicationprospect.Honeybeeischosenforthestudyobjectofthepaper,andtheflappingflightmechanismsofhoneybeeshavebeensystematicallystudiedwithwindtunnelexperimentsandnumericalsimulations.Mainworksareasfollows:(1)Low-speedwindtunnelusedforflowfieldvisualizationofinsectshadbeendesignedandmanufactured.Inordertoimprovetheflowfieldqualityofwindtunnel,theoptimumsizesoftestsection,diffusersectionandcontractionsectionofwindtunnelaccordingtocharacteristicsofflowfieldwereresearchedanddetermined,whilewitozinskycurveselectedascontractioncurve.Honeycombsandnylonscreenswereaddedinsettlingchamberandpowersectionofwindtunneltoreduceturbulenceintensityandensureflowsteady.Thedesignedsmokedevicewasinstalledinthetunnel,andhigh-speedcamerasandlightswerearrangedreasonably,thusdevelopinganintactexperimentalsystemofflowfieldvisualization.(2)Flowfieldvisualizationexperimentsofhoneybeesflappingwereconductedbythelow-speedsmokewindtunnel.Withthehelpofmodelestablishedandimagessynchronouslyshotbytwohigh-speedcameras,characteristicsofwingflappingwereobtained,includingflowfieldstructuresoffarfieldandnearfield.Vortexevolutionofhoneybees’flappingflightwasresearched,andinfluenceofwakevorticesoninducedvelocitywasqualitativelyanalyzed.Theflowvisualizationexperimentshowedthat:continuousvortexwakesofwingtipformedvortextubeduringdownstroke;downstrokestopvortexandupstrokestartingvortexweregeneratedduringsupination;vorticitydirectionsofforewingandhindwingwasoppositecausedbythedeformationalongthechord;andleadingedgevortexwasthemostobviousandthesizeofleadingedgevortexIII 北京理工大学博士学位论文wasincreasingfromwingroottotipattheendofdownstroke.Duringtheendofdownstroketherewasavortexwing,whichwascomposedoftheleadingvorticesalongthespan,wakevortextubeofwingtipvortices,andwingrootvorticeswake,onleftandrightwingrespectively.Vortexringliftestimatedwaslargerthantheweightofhoneybee,suggestingthatliftgeneratedbywingflappingcansupporthoneybeetoflywithsomeextraweight.(3)Ahoneybeemodelanditscoordinatesystemwereestablished,withkinematicequationsofthemodelwingprovidedbasedontheregularpatternsofhoneybees’flappingwing.Theflowfieldofhoneybeeflappingwassolvedbydynamicmeshtechnique.Resultsofpostprocessingofflowfieldindicatedthatleadingedgevortex,wingtipvortex,wingrootvortexandtrailingedgevortexappearedonwings.Numericalsimulationresultsverifiedtheconclusionsofwindtunnelexperimentthat:attheendofdownstroke,thesizeofleadingedgevortexwasincreasingfromwingroottotipandleadingedgevortexmetwingtipvortex;attheendofdownstroke,downstrokestopvortexappearedatthetrailingedge,whichwoulddecreasethelift,whileatthebeginningofupstroke,upstrokestartingvortexattrailingedgeincreaseedthelift;onleft(right)wing,wingtipvortexwakes,wingrootvortexwakesandwingsurfacevorticityconstitutedtheintactvortexringwhichwascircular.Thenumericalsimulationresultsstillindicatedthatwingrotationmadethesizeofleadingedgevortexmaximumattheendofupstroke(downstroke),andattheendofupstrokewingtipvortexwake,wingrootvortexwakeandwingboundaryvortexmettoformtheEar-likevortexring.Inthecentersoftheupstrokeanddownstrokevortexring,thedownwardjetverticaltothevortexringplaneappeared.(4)Accordingtothewingentitiesofhoneybee,thecoordinatesystemofhoneybeewingmodelwasestablishedtodefinelift,thrust,aerodynamicpowerandinertialpower.Aerodynamicsandpowersofsingle-wingmodel,double-wingmodelandintegralofthehoneybeemodelwerecomparedthroughnumericalsimulations,indicatingthattheinteractioneffectsofaerodynamiccharacteristicsonsinglewingcanbeneglected.Meanwhile,throughalteringworkingconditions,suchasfrequency,maximumflappingangle,mid-strokeattackangle，reversetimingandpositionofrotationaxis,thevaryingpatternofaerodynamicforcesandpowerofthesinglewingmodelwereobtainedbyIV 北京理工大学博士学位论文simulations.Currentworkaimstomakecontributiontoresearchandmanufacturingofthebio-inspiredmicroflapping-wingaerialvehicle.KeyWords:honeybee;tunnelwind;vortex;flowvisualizationexperiment;numericalsimulationV 北京理工大学博士学位论文目录第一章绪论.....................................................................................................11.1研究背景及意义.....................................................................................................11.2扑翼研究的发展现状.............................................................................................31.3扑翼飞行的气动机制.............................................................................................81.3.1准稳态分析................................................................................................81.3.2合拢打开和剥离打开................................................................................91.3.3前缘涡与展向流......................................................................................111.3.4尾迹效应..................................................................................................151.4涡的基本理论.......................................................................................................171.4.1环量与涡..................................................................................................181.4.2Helmholtz第一定理.................................................................................181.4.3Biot-Savart定理........................................................................................191.5蜜蜂结构与飞行特点...........................................................................................201.5.1蜜蜂的结构..............................................................................................201.5.2周期运动..................................................................................................221.5.3拍动平面..................................................................................................231.5.4翅的扭转..................................................................................................241.6论文的研究内容...................................................................................................25第二章风洞结构参数对流场品质的影响规律研究...................................272.1引言.......................................................................................................................272.2风洞的整体结构...................................................................................................282.3风洞流场品质的评价指标...................................................................................292.4实验段尺寸对流场品质规律的影响...................................................................302.4.1实验段尺寸、模型及边界条件..............................................................302.4.2实验段截面边长对流场品质规律的影响..............................................322.4.3实验段长度对流场品质规律的影响......................................................35VI 北京理工大学博士学位论文2.4.4实验段尺寸的确定..................................................................................382.5扩散段尺寸对流场品质规律的影响...................................................................392.5.1扩散段尺寸、模型及边界条件..............................................................392.5.2扩散段长度对扩散段出口流场品质规律的影响..................................402.5.3扩散段长度对扩散段入口流场品质规律的影响..................................422.5.4扩散段尺寸的确定..................................................................................432.6收缩曲线及收缩比对流场品质规律的影响.......................................................442.6.1收缩段模型及边界条件..........................................................................442.6.2收缩曲线对流场品质规律的影响..........................................................452.6.3收缩比对流场品质规律的影响..............................................................482.6.4收缩段尺寸的确定..................................................................................502.7风洞整体结构参数及流场品质验证...................................................................512.8本章小结...............................................................................................................52第三章流场显示实验系统开发.................................................................533.1引言.......................................................................................................................533.2烟线显示的意义...................................................................................................553.3整流装置...............................................................................................................573.3.1喇叭口和稳定段......................................................................................573.3.2蜂窝器和阻尼网......................................................................................573.3.3动力段......................................................................................................583.4发烟装置...............................................................................................................593.5建立实验平台.......................................................................................................603.6准备实验对象.......................................................................................................623.7流动显示实验.......................................................................................................633.7.1圆柱绕流实验..........................................................................................633.7.2蜜蜂流场显示实验..................................................................................643.8本章小结...............................................................................................................65第四章基于蜜蜂流场显示实验的涡演变规律...........................................66VII 北京理工大学博士学位论文4.1引言.......................................................................................................................664.2蜜蜂的振翅规律...................................................................................................684.2.1直线与投影关系......................................................................................684.2.2蜜蜂的坐标系..........................................................................................694.2.3蜜蜂的拍动角..........................................................................................714.3蜜蜂振翅的远场流场结构...................................................................................724.3.1蜜蜂周期振翅的流场显示实验分析......................................................724.3.2下拍停止涡与上拍启动涡......................................................................784.3.3翅扭转产生的尾缘涡对..........................................................................804.4蜜蜂振翅的近场流场结构...................................................................................834.4.1蜜蜂周期振翅的流场显示实验分析......................................................834.4.2前缘涡......................................................................................................854.4.3下拍停止涡与上拍启动涡......................................................................874.4.4翅尖涡......................................................................................................874.5蜜蜂振翅的涡.......................................................................................................894.5.1翅旋转阶段产生的尾缘涡......................................................................894.5.2下拍涡环与平均升力..............................................................................904.6本章小结...............................................................................................................92第五章蜜蜂飞行时涡演变规律的仿真验证...............................................935.1引言.......................................................................................................................935.2蜜蜂模型的坐标系...............................................................................................955.3翅的运动方程.......................................................................................................965.3.1简谐运动方程..........................................................................................965.3.2中期拍动速度和攻角不变的方程..........................................................975.3.3提前、对称和延迟翻转的运动方程....................................................1005.4边界条件和求解器.............................................................................................1025.5计算有效性验证.................................................................................................1055.6涡现象分析.........................................................................................................107VIII 北京理工大学博士学位论文5.6.1蜜蜂振翅产生的涡................................................................................1075.6.2前缘涡....................................................................................................1105.6.3下拍停止涡与上拍启动涡....................................................................1125.6.4涡环........................................................................................................1135.7本章小结..............................................................................................................116第六章蜜蜂翅模型的力和功率分析........................................................1186.1引言.....................................................................................................................1186.2气动力和功率分析.............................................................................................1196.2.1翅的气动力............................................................................................1196.2.2翅的功率................................................................................................1206.3蜜蜂身体、双翅和单翅对力和功率的影响.....................................................1226.3.1单翅与双翅对力和功率的影响............................................................1226.3.2蜜蜂身体对气动力和功率的影响........................................................1266.4典型算例.............................................................................................................1276.5力和功率的影响因素.........................................................................................1316.5.1频率........................................................................................................1316.5.2最大拍动角和中期攻角........................................................................1326.5.3T1..............................................................................................................1376.5.4T2..............................................................................................................1396.5.5翻转时刻................................................................................................1406.5.6旋转轴位置............................................................................................1426.6本章小结.............................................................................................................144结论与展望...................................................................................................146参考文献.......................................................................................................150攻读学位期间发表的论文与研究成果清单...............................................165致谢...............................................................................................................166IX 北京理工大学博士学位论文第一章绪论1.1研究背景及意义微型飞行器（MicroAirVehicle―MAV）是20世纪90年代中期源于军事目的而[1]发展起来的一种新型飞行器。MAV一般是指无人自动控制的飞行器，尺寸一般在[2-3]15厘米以下。由于尺寸小，可操作性强，能够独立作战，其在军事以及民用领域[4]均有其特殊用途，具有极好的发展前景。MAV可以看起来像是鸟类，也可以像是蝴蝶或者蜻蜓，蜜蜂等昆虫，该类飞行器最大的优势就是体积尺寸小，便于隐蔽，能够[5]在狭窄的空间内发挥作用，这是普通的飞行器无法与之相比的。凭借这些优点，MAV可以在各种复杂地形、恶劣环境（例如沙漠，丛林，海洋，高山和北极）下单独工作或者群体作战，也可以飞行到指定位置破坏电路，甚至能够深入到隧道或者深埋的目标执行普通飞行器无法完成的任务。现代MAV广泛应用于商业，科研，政府和军事等领域。MAV可以远程摇控，能[1]够到达人，地面车辆或其他设备无法进入的危险环境。与常规无人飞行器相比，MAV具有体积小、重量轻、成本低的飞行平台优势，操纵方便、机动灵活、噪音小、隐蔽性好，具有巨大应用潜力，因此许多发达国家己经开始对其进行重点研究。在军事方面，MAV可进行低空军事侦察、监视、战场损伤评估等；作为反辐射和微型攻击武器，摧毁敌方雷达等电子设施以及携带微型战斗部进行攻击；用于目标搜索和通信中继；进行生化探测并标定危险区域等。在民用方面，MAV可以用于交通监控、边境巡逻、森林及野生动植物勘测、航空摄影、输电线路检查、环境监测、气象监测、森[6]林防火监测等。MAV的特征尺度远小于传统的飞行器，不能照搬传统飞行器设计理论和方法，[7]飞行器的微型化将面临诸多的来自科学和技术上的挑战。MAV分为固定翼飞行器、[8]旋翼飞行器和扑翼飞行器三种。与固定翼和旋翼飞行方式相比，微型扑翼飞行器具有明显的优势，其具有较高的升阻比，抗干扰性好，可以快速的起飞、加速、悬停，[9]具有高度机动性和灵活性，飞行效率高，可以用很小的能量进行长距离飞行，因此[10-11]是最有发展空间的一类MAV。扑翼飞行被自然界中的鸟类和昆虫广泛所采用，被认为是生物进化的最优飞行方1 北京理工大学博士学位论文式。研究微型扑翼飞行器，可以基于仿生学原理，借鉴自然界中的飞行生物的飞行特点和原理。鸟类和昆虫能够在强风和复杂环境下稳定飞行，具有很高的飞行机动性和[1]稳定性。因为在飞行过程中，当外界条件变化时它们的翅膀和身体可随之产生自适应变形，以保持最有利的飞行姿态。飞机能够在高空中飞行，主要靠机翼产生升力。气流到机翼前缘，分成上、下两股气流，机翼上表面气流速度较快，压力低，而机翼下表面，气流受阻挡作用，压力[12]增大。于是机翼上、下表面出现了压力差，压差的总和就是机翼的升力。飞机就是借助机翼上获得的升力克服自身重力，从而翱翔于蓝天。传统空气动力学中，飞机的[13-14]飞行速度很高，雷诺数（Re）也高，机翼周围的流场是稳态的，固定翼飞机的升[15-16]力原理、推进方式及可操控性的研究已较为成熟。扑翼飞行的难点主要集中在低雷诺数下的非定常流场研究，因为随着翅的不断拍动，周围的流场也在不停变化，难以用解析方法讨论和解决扑翼飞行的气动性能。扑翼飞行是微尺度下低雷诺数的非[17-18]定常流场问题，该方面研究面临着复杂和困难的挑战。自然界中的昆虫具有体积小、动作灵活、能够感知外环境的变化并快速响应等优势。很多学者和研究人员从飞行的昆虫或鸟类中获得灵感，有力推动了微小型飞行器的发展。诸如昆虫扑翼运动的非稳态空气动力学，神经传感分布和整合，以及信息处理技术，都是在研究微小型飞行器中值得借鉴的。昆虫是自然生物进化中一个完美的例子，它们尺寸微小但仍可探知世界，它们可以飞行、跳跃、爬行、追逐及逃跑，甚至能够找回自已的巢穴。在飞行过程中，昆虫的神经系统不仅可以在高速飞行时协调昆虫运动，而且能够监测周围环境的变化，使其在很小的落地点精确着落，能在不稳定的气流中平衡姿态、探寻食物和栖息场所，甚至可以指引同伴回巢路径。不同于传统飞行器的飞行方式，昆虫在飞行过程中扇动翅膀，整个胸腔会快速收缩，翅膀就会被动地随之进行拍动，周围的气流会不断的产生变化，伴随着涡的产生和脱落，昆虫会根据环境的变化调整飞行姿态和方向，以适应复杂的外界条件。受昆虫启发，工程师们制造微小型扑翼飞行器（Flapping-wingMicroAirVehicle，简称FMAV），他们希望FMAV可以像昆虫一样机动灵活，能够实现起飞、降落、转弯、悬停、前飞、侧飞等动作，在完全湍流和低雷诺数环境下能够利用扑翼产生高升[19]力进行飞行。FMAV应该像昆虫一样，具备体积小、能耗小等优点，在复杂环境中2 北京理工大学博士学位论文有较好的飞行稳定性，有较高的机动灵活性。为了解决这些问题，借鉴昆虫的飞行，研究昆虫的飞行机理以及空气动力学原理，是非常必要的。1.2扑翼研究的发展现状上世纪90年代后，扑翼飞行已成为一个热门领域，各行各业的专家包括航空机械工程师、机器人专家、流体力学专家、数学家以及生物学家用不同的方法对其进行[20]研究。与固定翼飞行方式相比，扑翼飞行速度低，雷诺数低，空气的粘性较高。对扑翼飞行机制的研究主要包括实验研究和计算模型数值模拟研究。随着高速摄像机和传感器的发展，计算机性能的逐渐提升，学者们多通过研究机械仿生扑翼实验和数值计算模拟扑翼飞行流场来研究扑翼飞行的机理。北京航空航天大学孙茂教授团队基于CFD方法，对昆虫扑翼飞行的气动研究做了大量工作。他们基于拟压缩性思想，即在连续方程中加入压力对拟时间的偏导数项，每个物理时间步用拟时间进行内迭代，对流项采用Roe格式后进行通量差分裂，用线[21][22][23]性Gauss-Seidel进行迭代求解N-S方程，计算了不同种类昆虫（蝴蝶，蜻蜓，[24-25][26-27]果蝇2等）在悬停和前飞时的气动力和功率，研究了昆虫动态稳定的飞行控制[28-29]问题。孙茂等人通过研究悬飞的果蝇，提出了高升力的三个机制，快速加速机制，[30]延迟失速机制和快速旋转机制；利用三台高速摄像机拍摄了自由悬飞的蜻蜓和食蚜蝇，测量了翅和身体的形态参数，根据所测量的运动参数利用N-S方程研究了蜻蜓和[31-32][33]食蚜蝇的气动特性。此外，还研究计算了悬停昆虫翅变形对气动力的影响，变形翅比刚体平板翅产生的升力较大（高10%），需要的能耗较少（低5%）。西北工业大学设计并制作了专用于研究微扑翼飞行器特性的低湍流风洞，研制了扑翼飞行器样机，对其进行了多方面研究。刘岚等人通过对扑翼飞行器进行风洞实验，[34]测量了升力。张西金等人在仿蜜蜂类昆虫动力学方程基础上，采用分层控制策略研[35]究了悬停控制问题。杨智春等人通过对微型扑翼机进行空气动力特性试验，得出了[36]影响升力的因素，并提出柔性扑翼产生的升力和推力较高。谢辉等人应用有限体积法结合双时间推进技术求解了微型扑翼在低速、低雷诺数下的粘性绕流，研究了扑翼[37]参数对气动特性的影响。熊超和宋笔锋利用研究微型飞行器的专用风洞，研究了迎[38-39]角、扑动频率和风速对气动特性以及纵向力矩特性的影响。侯宇等人根据仿生学结果设计了飞行器的传动布局、动力方案以及总体结构，并对微扑翼飞行器进行了优[40]化设计。3 北京理工大学博士学位论文南京航空航天大学昂海松等人对自行研制的微型柔性结构扑翼机进行了风洞实[41]验，并成功进行了飞行试验。曾锐基于流固耦合方法，利用有限体积法离散任意拉格朗日欧拉N-S方程，使用Delaunay图映射方法实现网格变形，研究分析了仿鸟扑翼的流场涡结构和升力推力的产生原理，分析了仿鸟复合振动的扑翼气动特性，在原有匀速刚性模型的基础上，提出考虑了扑翼扑动速率变化和形状变化的柔性扑翼模[42]型，使之更接近鸟翼真实扑动情况。昂海松和肖天航教授等人模拟计算了柔性扑翼气动功率及扑动效率随着扑动角、来流速度等参数的变化，从气动角度解释了为何鸟[43-44]类在不同的飞行阶段扑翼规律各不相同的问题，并对蜻蜓柔性翼进行了气动计[45-46]算，研究了蜻蜓前后翼相位差对飞行姿态的控制。[47]清华大学吴子牛教授等人基于Rayner的扑翼涡理论模型，利用涡动量理论，计算了脱落的多个尾迹涡环对平均升力的影响，提出尾迹涡诱导拍动平面处的下洗速[48]度会降低升力，同时对昆虫身体对升力的影响进行了理论研究，研究表明身体会降[49]低升力，身体尺寸影响了升力降低的数值。清华大学的曾理江教授等人利用条纹投影法，从两个互相垂直的方向投射条纹，测量昆虫的拍打运动，重点进行了昆虫运动机理研究和应用以及有关昆虫运动参数的[50-51]测量和分析，建立了昆虫运动模型，研究了昆虫的运动规律，并测量了蜻蜓前飞[52]和转弯时翅的运动参数、飞行轨迹和身体姿态。上海交通大学的陈文元和张卫平教授等人从仿生角度研究了MAV的飞行机理与样机制作，通过研究昆虫的运动特征与流场结构，为设计仿生微型扑翼飞行器提供了[53]参考。左德参基于MEMS技术和仿生学原理，设计了微扑翼飞行器结构，对其数[54-55]学建模并进行了仿真研究。彭松林利用MATLAB和ADAMS建立了曲柄摇杆扑翼机构模型，仿真分析了减速比和曲柄摇杆参数变化对频率的影响，提供了优化设计[56]方案。中国科学技术大学童秉纲教授和陆夕云教授等人对飞行和游行的生物进行了运[57]动学研究。鲍麟和童秉纲数值求解了二维昆虫翼的动态变形对气动性能的影响，计算结果表明，影响气动特性的主要原因是附加惯性效应，动态变形对平均气动性能影[58]响很小，但是影响瞬态的气动性能，这有可能与飞行的稳定性有关。余永亮和童秉纲理论分析了昆虫翅的轨迹对前飞气动性能的影响，指出8字型轨迹的气动性能优于[59]O型轨迹的气动性能。国外对扑翼飞行的研究，至少可追溯到1919年，工程师Hoff利用昆虫学家Demoll4 北京理工大学博士学位论文[60][61]的数据，得出自然界生物不能使用固定翼飞机的气动理论进行飞行的结论。Demoll认为，利用Hoff的方程，熊蜂将需要至少两倍于固定翼飞机的升力系数才能[62]飞行。随着近些年来对扑翼空气动力学理论的研究，人们才发现，解释昆虫飞行的高升力需要利用非定常气动理论，这与固定翼的定常理论是有本质区别的。Weis-Fogh[63-64][65-66]研究了昆虫的非定常机制。他观察分析了台湾小黄蜂的气动性能，当小黄蜂上拍结束时，双翅绕前缘合拢在一起，然后绕着翅后缘打开，接着向下拍动，他把这种拍动运动称为“合拢、打开（clapandfling）”。这种方式加速了下拍过程中环量的[67][68-69]产生，在下拍过程中会提供很高的升力，这个气动效应被称为Weis-Fogh效应。虽然合拢并打开这种方式被某些昆虫所采用，但是后来更多研究表明，对于大多数自然界的飞行生物而言，这种气动机制并不具有代表性。[70]Ellington在1984年提出用“准稳态假设理论”分析扑翼运动的气动性能，并[71]对不同种类飞行生物的形态学参数总结分类，利用高速摄像机拍摄记录了昆虫的飞[72]行参数，系统地建立了昆虫的飞行动力学模型，针对扑翼运动时的空气动力学特性，[73]分析了几种气动机制，尤其提到了利用涡理论来解释昆虫和鸟类的悬停飞行，并且[74]给出了估算平均升力和能耗的公式。随后，Ellington团队对昆虫的飞行及气动性能进行了多方面的研究。Willmott将鹰蛾置于开口风洞中，对其进行流场显示实验，利用立体相机技术拍摄了不同风速和展向位置处的烟线图，分析并提出鹰蛾飞行时有一[75][76]系列的垂直和水平涡环。Gilmour估算了熊蜂飞行的能耗。Wakerling和Ellington[77]研究了蜻蜓在稳态下的气动力、飞行速度和运动参数，还估算了蜻蜓的升力和能耗[78-79]。VanDenBerg制作了机械昆虫flapper，模拟鹰蛾悬飞的翅的运动，在机械翼前[80-81]缘处释放烟线，能够观察到前缘涡尺寸从翅根至翅尖是增长的。Ellington提出对于给定质量的MAV，尺寸较大的MAV需要的能量较少，但是体积小的MAV能够利[82]用高频拍动产生更快的速度。Tyell基于鹰蛾的动力学和形态参数，利用活塞管产[83]生了类似鹰蛾悬停时的涡环，在层流状态时可估算尾迹涡的环量和冲量。Dickinson团队开发了油桶实验系统平台，将机械果蝇扑翼模型置于油桶中拍动，拍动过程分为平动和转动两部分，通过传感器测量扑翼在各种工况下的速度和压力，[84]并提出了延迟失速、旋转环流和尾迹捕获三种机制。实验结果表明，扑翼模型在上拍初期有一升力峰，平动过程升力近似保持不变，在下拍末期又出现一升力峰。下拍与上拍过程中升力变化相似。Dickinson提出，平动过程中产生的较大升力是由于延迟失速机制，半拍中第一个峰值是由于之前形成的尾迹气流影响产生的，称为尾迹捕5 北京理工大学博士学位论文获机制；半拍中第二个峰值是由旋转环流机制产生的。模型扑翼旋转阶段的升力平均值为一周期升力总和的35％，这说明尾迹捕获和旋转环流是昆虫获取飞行所需升力的重要机制。这两个部分升力不仅提供了昆虫飞行所需的升力，而且提供了昆虫控制其飞行方向的手段。Sane利用等比例的果蝇扑翼模型，研究了翅的运动学参数对非稳态[85][86]气动力的影响，研究了扑翼模型的旋转对气动特性的影响。Birch等人利用油桶[87]实验，研究了在高雷诺数和低雷诺数下扑翼的前缘涡结构和力的产生，利用DPIV[88]技术研究了尾迹效应对扑翼飞行气动特性的影响。Fry等人利用三台高速摄像机拍摄了自由飞行的果蝇，捕获了果蝇的翅和身体的运动参数，将运动参数用于模型机械扑翼，测量了产生的气动力，随时间变化的力矩表明控制昆虫飞行动力的主要是惯性力，而不是摩擦力，对果蝇悬飞时的气动力进行了研究，模型扑翼在相同运动下所需[89-90]的能量大于果蝇。Lehmann用果蝇的模型扑翼研究了双翅交互作用对气动性能的[91][92]影响。Lentink提出了翅在旋转时加速，可以稳定翅上的前缘涡。Liu基于有限体积法利用多块overset网格法求解N-S方程，计算了鹰蛾悬停时翅周围的流场，通过后处理流线显示分析了前缘涡和展向流，计算结果与Ellington机械[93]翅的实验对比基本吻合，由此确定了前缘涡和展向流对翅产生高升力的重要性。Liu还基于真实的昆虫形态学尺寸和运动参数建立模型，计算了鹰蛾，蜜蜂，果蝇三种昆[94]虫悬停飞行的流场。Aono计算了昆虫在悬飞时近场和远场的三维流场情况，计算[95]的平均力和一个周期内的升阻力变化情况与实验结果相近。Aono和Liu通过显示一周期内涡结构的变化，提出拍动过程中形成的涡环是持续保持高升力的重要原因，前缘涡和翅尖涡对升力具有重要贡献，此外，还基于气动力和惯性力，估算了鹰蛾悬[96]停飞行时的能耗。Nakata计算了鹰蛾翅在变形状态下的气动力、惯性力、力矩和功率，计算结果表明翅的变形可以提高昆虫悬停飞行的效率，他指出昆虫的气动性能由[97]很多因素决定，包括翅的强度、翅脉的流通能量以及翅的运动参数等。随着近代测量技术的发展，特别是PIV技术的广泛应用，更多研究学者倾向于通过实验研究扑翼飞行的流场特征来分析扑翼的气动特性。在定量的流场显示实验中，不仅能观察涡流随扑翼运动的变化规律，还能根据测量的速度场计算环量和力的大小。粒子图像测速（PIV）是一种光学技术，即在流场中释放小的示踪粒子，通过追踪粒子运动测量速度，从而测量流场的速度场。最初这种技术用于固体力学领域，通过跟踪斑点图案测量表面运动。Barker和Fourney提出，该技术可用于测量流体速度[98][99-100]。Adrian和Prasad等人进一步研究了PIV技术的发展及应用。6 北京理工大学博士学位论文Spedding利用PIV技术测量了鸽子低速飞行时尾迹提供的动量与能量，研究了多种鸟类飞行的尾迹结构对扑翼飞行的影响，提出了低成本，高精度的DPIV测量方法[101-102]专门用于研究鸟类的扑翼飞行。他基于涡动量方程，改进了PIV测量力的方法，[103]分析了雨燕和蝙蝠的尾迹结构和对力的影响。Murphy等人用PIV测量了平板、波纹蜻蜓翼和基于波纹翼轮廓的光滑翼型的流场，对比了3种翼形不同雷诺数下的升阻力系数。实验结果表明，光滑翼型的升力系数受雷诺数影响大；波纹翼形升力系数几乎不受雷诺数影响。Murphy的PIV测量数据表明，对于波纹蜻蜓翼而言，在褶皱处有非稳态涡填充，提供了更高的能量，有助[104]于边界层的附着，从而阻止了气流分离。Poelma等人最先发表了扑翼随时间变化的三维PIV数据。他将果蝇翅模型置于矿物油中，针对锁定的时间阶段，测量了扑翼运动的流场数据，通过整合每个时刻单个截面的PIV数据信息，重构了扑翼随时间变化的三维流场信息。Poelma同时研究了翅模型瞬间启动的过程，首先，围绕翅的惯性环量的增长几乎是线性的，没有展向流；在一段时间后，展向流将环量从翅尖带走，前缘处的环量达到平衡，前缘涡能够[105]稳定地附着于翅上。Lu等人利用PIV技术，将模型扑翼置于水中，研究了扑翼模型在不同雷诺数、不同展弦比和不同攻角下的前缘涡结构，实验结果表明，在高雷诺数、高攻角和小展[106]弦比时，前缘涡系包括一个主涡和一个次涡。Lu还研究了模型蜻蜓翼前缘涡的演变过程，分别测量了翅模型在加速、减速和最大速度时不同展向位置处的速度场，重构了三维流场结构，实验结果表明，模型蜻蜓翼的前缘涡系统包括一个主涡和三个小[107]涡。扑翼拍动产生的气动力与多个因素有关，包括运动参数、翼的平面形状和周围流体特征。扑翼表面与附近空气的流固耦合交互作用，也会对气动力有一定影响。但是至今为止，用数值模拟计算气动弹性的相互作用都是很困难的，在上述的大部分数值计算求解扑翼周围流场中，均是刚体翼模型。在多数机械扑翼测量实验中，采用的也是刚体翼模型。但是，在自然界中，昆虫翅在周期拍动时，伴随着翅的变形运动。受此启发，越来越多的学者开始研究翼的变形对气动特性的影响。Combes利用真实鹰蛾翅，将其按真实频率和拍动幅度运动，对比了翅在空气和氦气中运动时翅的弯曲程度，结果表明翅的弯曲程度相差不大。Combes认为影响翅弯曲的主要因素不是气动力，而是流体阻尼作用，因此他倾向于利用阻尼有限元模型7 北京理工大学博士学位论文[108]研究翅的变形特性。Combes还研究了翅的刚度和变形，他指出昆虫翅脉络的分布[109-110]可能影响翅的刚度，影响翅的瞬态形状，从而影响飞行性能。Combes的研究还[111]指出，当蜻蜓翅被损坏时，会影响飞行性能和捕食能力。Wootton等人利用高速摄影和机械测量技术对沙漠蝗虫进行了研究，蝗虫飞行时，后翅随着翅的拍打会变形，变形程度依赖于翅的结构。在快速前飞时，后翅在下拍初期通过打开机制产生高升力，在平动阶段，翅的弯曲变形也能产生高升力。实验结果[112]还表明，翅的脉络分布对气动力有重要作用。Hu等人利用三种不同的材料制作翅模型，分别是木头（完全刚性），尼龙（有一定柔性），乳胶（柔性最好）。乳胶翅的气动性能最佳，产生的推力最高；刚性翅能产[113]生最高升力；尼龙翅的气动性能最差。Nakata通过求解非线性流场，计算分析了昆虫柔性翅的三维流场。对鹰蛾悬飞的流场求解结果表明，柔性翼会延迟前缘涡的脱落时间，因此能产生较高升力；展向变形和弦向变形均能增强前缘涡，提高扑翼运动[97]的气动特性。Mountcastle的实验数据表明，翅的柔性变形可能影响鹰蛾的诱导气流。[114]他发现，与刚性翼相比，柔性翼的拍动幅度更大，力的方向更有助于提高升力。Du计算了昆虫悬停时翅的变形对气动力的影响，柔性翅比刚体翅产生的升力高10%，[33]能耗少5%。在对柔性翼的研究中，大部分结论倾向于柔性翼能提高扑翼运动的气动性能。但是在某些研究中，柔性翼提升的气动性能很有限，有时甚至不如刚性翼的气动性能好。这一方面可能是由于机械扑翼结构被简化，另一方面是因为CFD计算模型与真实昆虫翅有区别。1.3扑翼飞行的气动机制1.3.1准稳态分析准稳态分析方法通过多个瞬态的升阻力工况，估算流场中翼的平均升阻力。[115]Drzewiecki最先基于准稳态分析，提出了叶素理论，用于研究螺旋桨的气动性能。[116]Osborne将这一理论应用扩展到研究扑翼飞行机制。叶素理论将机翼划分为多个翼展方向的单元，其中每个单元的升力L′和形阻D′为：pro12L′=rcUC（1.1）airrL212D′=rcUC（1.2）proairrD,pro28 北京理工大学博士学位论文r是空气密度，c是弦长，U是拍动速度和诱导速度的总和，C和C，分airrLD,pro别为升力系数和形阻系数。忽略相对速度的展向分量对升力和形阻的影响，在准稳态假设理论中，C和C仅与雷诺数和攻角有关，需要利用实验的方法预先测量出CLD,proL和C。由于不能预估诱导速度，叶素理论模型更适用于研究快速前飞，因为在较D,pro高的速度下，诱导涡量能更快地对流到下游，比起较低雷诺数下的悬飞，对翅周围的气流影响很小。Whitney和Wood利用叶素理论估算了微型仿生昆虫机械翅的升力，结果表明，叶素理论的估算值与实验测量的结果非常吻合。但是，准稳态分析方法不能解释延迟失速和翅反转产生的脱落涡，在一些实验中，准稳态估算的升力与实验结果并不一致[117]。Hubel和Tropea通过模拟鹰蛾的飞行，测量了模型扑翼产生的力，并利用PIV[118]技术测量了尾迹，实验测得的力比准稳态估算的力要高。这表明在扑翼飞行中，存在非稳态机制，准稳态假设是不适用的，且当雷诺数越低时，准稳态假设的结果就越差。准稳态分析模型往往会低估阻力和功率能耗，但是准稳态的预估值可通过改进模型来完善。Sane和Dickinson测量了模型果蝇翅在不同拍动幅度和攻角下的平均升力、阻力和功率，将实验结果与准稳态估值对比，准稳态模型虽然较为准确地预估了升力，[86]但是对阻力的估值明显偏低。Sane和Dickinson通过隔离分析果蝇翅模型的平动和转动，改善了模型，结果表明，在一个周期中，除旋转阶段的其他时间，准稳态分析[85]方法均能较好地预估改进模型的升力和阻力。由于尾迹捕获机制，旋转阶段无法用准稳态分析模型准确预估翅模型的气动力。Ellington基于准稳态分析方法，把昆虫振翅的拍动过程分离成若干个时刻，通过测量计算每个时刻的稳态力，得到了力随时间的变化历程。在每个时刻，气动力等同于稳态下以相同速度和攻角产生的气动力。Ellington指出，如果这种模型准确，用一系列简化的方程即可计算翅的气动力。但是，根据大量的实验数据，准稳态模型预测的最高升力小于昆虫悬停时的平均升力。因此，他提出，准稳态模型是不正确，或者[70,119]至少是不充分的，需要修改或者用其他模型研究昆虫扑翼飞行的机制。Ellington认为，在翅的旋转阶段，准稳态理论需要修改，可能需要非稳态机制来研究。在此之后，更多的学者认同Ellington的观点，提出需用非稳态气动机制研究扑翼飞行的气动特性。1.3.2合拢打开和剥离打开9 北京理工大学博士学位论文合拢打开（clapandfling）和剥离打开（peelandfling）是相似的，均是描述飞行生物的拍翅运动，左翅和右翅上拍到最高处时，左、右翅距离非常近（clap），或者靠[120]在一起（peel），之后像拉拉链一样，从前缘绕后缘快速打开（fling）。昆虫在达到向上拍翅和向下拍翅最大限度时，也就是分别在拍动平面的最顶部和底部，此时左右翅距离非常的近，这称之为合拢。最早是由Weis-Fogh观察到台湾小黄蜂时提出。当小黄蜂上拍结束时，双翅绕前缘合拢在一起，然后绕着翅后缘打开，[66]接着向下拍动，他把这种拍动运动称为“合拢打开”，如图1.1所示。这种方式加速了下拍过程中环量的产生，在下拍过程中会提供很高的升力，这种气动效应被称为Weis-Fogh效应。图1.1合拢打开三维示意图剥离打开（peelandfling）是描述一些昆虫下拍的变形和运动，尤其是拍打频率[121]低且尺寸小的昆虫，例如蝴蝶，石蝇，草蜻蛉等。剥离打开指的是由于上拍结束时，左、右翅有部分是重叠在一起的，所以下拍开始时，重叠在一起的部分会分离，在两翅尖随翅其他部分分开之前，翅绕尾缘旋转时分离。剥离打开这种方式对多数昆虫的影响较小，而且无法用数值模拟方法评估这种机制的影响。这种拍动方式对气动性能的影响主要也在于涡的产生和气流分离，能够直接影响翅面的环量，形成附着涡。左、右翅的合拢打开和剥离打开机制主要影响涡的产生，翅在旋转（特别是提前旋转）时，加速了环量的产生。左、右翅的打开和剥离能有效地产生翅面的附着涡。当气流流经翅的上下表面时，打开会直接引起强的展向流动。Ellington认为，打开的[73]主要作用在于产生强的前缘涡。Grodnitsky对飞行昆虫周围的流场实验表明，前缘处的气流分离和形成的涡主要产生于较高雷诺数下，在自然界多数昆虫飞行的速度范[121]围内是不会发生的，这也可能只是翅翻转阶段短暂的机制效应。左右翅打开时有较强的展向特征，即气流会翅尖流向翅根。图1.2为合拢打开的二维示意图。在合拢过程中，左右翅的前缘先于尾缘合拢，直至双翅完全合拢，如图1.2-A和图1.2-B所示。当双翅闭合时，如图1.2-C所示，每个翅上的涡量相反彼此抵消。这使得下次拍动时，每个翅上脱落的尾缘涡不存在。根10 北京理工大学博士学位论文[122]据Wagner效应，脱落的尾缘涡会延迟环量的增长。因此Weis-Fogh提出，尾缘涡[67]的缺失或衰减，加速环量的增长，在连续的拍动中提高升力。双翅合拢时产生的射流能为昆虫提供飞行的推力，如图1.2-C所示。在合拢末期，翅旋转时，前缘先于尾缘打开，如图1.2-D和图1.2-E所示。双翅中间产生了低压区，周围的气流进入这个区域，提供了能量，从而形成了附着前缘涡。之后，双翅带着相反的边界环量，继续拍动并分开。虽然每个翅的附着环量均能提供升力，但两个翅的总环量为0，如图1.2-F所示，因此满足开尔文环量守恒定理。Lighthill提出，在完全粘性流体中，开尔文环[69]量守恒定理适用于打开（fling）机制。总体而言，合拢打开机制能显著提高升力。图1.2合拢打开二维示意图但是，合拢打开机制这种假设是基于昆虫能达到最大拍动幅度时，升力才能明显提高。一些研究表明，在鸟和昆虫飞行中，通过增大拍动角可以提高升力，拍动幅度最大达到180°。因此，昆虫双翅的合拢运动，可能只是试图增大拍动角。虽然合拢打开和剥离打开机制被一些飞行昆虫所采用，但是自然界中采用这种方式做拍打运动的昆虫种类仅有很小的一部分，大多数昆虫和鸟类不采用这种飞行方式，因此，合拢打开和剥离打开机制并不具有代表性。合拢打开和剥离打开机制的提出对扑翼的非定常研究起了先导性作用，引导了后继研究学者们对扑翼飞行的非定常机制进行探索。1.3.3前缘涡与展向流在传统空气动力学中，当机翼攻角大于某一临界值时，翼型上表面气流产生严重分离，机翼产生的升力急速下降，不能继续稳定飞行，这称为失速，此时的攻角称为[123]临界攻角。Weis-Fogh在研究合拢打开模型时，提出了延迟失速机制，首次提出11 北京理工大学博士学位论文扑翼产生较高的升力并保持稳定，是由于翅前缘处附着的前缘涡（Leading-edge[67]Vortex，简称LEV）。之后，众多学者对延迟失速和前缘涡理论展开了研究。在Wanger实验中，翼型突然改变工况，攻角大于临界攻角后，翼型能继续运动几个弦长[124]的位置，这种效应被称为延迟失速机制。与准稳态模型结果相比，在失速之前，翼型有更大的环量，这与前缘涡的增长有关。Usherwood对模型翼进行了流场显示和测力实验，LEV的存在同时提高了模型翼[125]的水平力与竖直力，证明了延迟失速的存在。延迟失速机制能显著增加升力，提高扑翼运动的气动性能，但不能像准稳态理论一样预估出升力数值。吴江浩和孙茂数值模拟了模型果蝇翅的雷诺数在13～1500区间的数据，根据数值结果，当雷诺数大[126]于100时，前缘涡和延迟失速机制是产生升力的主要原因。翅的旋转运动能建立向外的展向流动，将翅上的涡量带到翅尖，这能稳定LEV，有助于延迟失速机制的产生。Hong利用简化的铰链平板扑翼结构，研究流场中涡流的特性。Hong沿着平板翼的展向取n个截面，利用PIV技术，研究了涡流与升力的[127]关系，结果表明涡量场与扑翼的展弦比有关，翼尖附近的弦向涡量能提高升力。此外，他还对比了平板翼和弧面翼的实验结果，证明了展向流动能提高升力，且弧面[128]翼的升力更高。[68]Maxworthy对Weis-Fogh机制展开了研究，他在流场中释放有色示踪粒子，通过流场显示实验，得到了沿着平板翼展向的LEV结构。实验观察到机翼表面的单独旋涡，这有利于改变翼的有效弯度和厚度。翼展方向的压力梯度建立了由翼根至翼尖的展向流动，将前缘涡的涡量从最内的展向位置（翼根）带到最外的展向位置（翼尖），最终将其带入翅尖涡。[129]Luttges研究了悬吊的蜻蜓、鹰蛾以及机械仿生翅的LEV结构，与Maxworthy和Weis-Fogh的结论不同，Luttges结果显示，LEV沿着翼展的尺寸是相同的，没有观察到螺旋状或者锥形的LEV结构。Luttges认为，在LEV涡核内几乎没有展向流动，[130-131]沿着翼展，流动特征结构可看作是二维的。Srygley等人对蝴蝶进行了流场显示实验，结果与Luttges的结论相同，沿着翅的展向，LEV尺寸几乎不变的，这表明展[132]向流很弱。Srygley和Taylor等人认为，斯特劳哈尔数优化了翅的运动，LEV的稳定并不一定需要有展向流动。在拍动的大部分阶段，LEV增长不脱落，这可能是因为[133-134]拍动周期比LEV变得不稳定所需的时间要长。12 北京理工大学博士学位论文Ellington等人对仿鹰蛾机械扑翼模型进行了流场显示实验，实验观察到扑翼模型[80]有展向流动，这使得LEV在扑翼平动过程中保持稳定不脱落。Willmott和Ellington等人利用烟线实验研究了悬吊的鹰蛾，根据实验结果，Willmott提出，LEV是三维结[80]构，从翅根至翅尖，LEV尺寸逐渐增大，展向流是促使附着前缘涡稳定的必要条件。[68]与Maxworthy的结果不同之处在于，在Willmott和Ellington的实验中，在低于25%翼展位置时，未能观察到LEV，据此预测左右翅上的LEV并不通过昆虫的胸腔相连接。此外，Ellington等人观察到，由翅根至翅尖，LEV涡核内有较高的展向流速度，[135]展向流可能是由于翅的离心加速运动或是旋涡引起的诱导速度场。为研究恒定角速度的纯旋转运动，Lentink等人给出了非惯性坐标系下无量纲化的Navier-Stokes方程。通过理论分析，Lentink等人提出，翅根至翅尖的压力梯度导致了展向流动，向心科氏加速度通过调整平动方向的力能保持展向流动，展向流使[136]LEV附着于翅上。Aono和孙茂等人通过数值模拟，也认为向心科氏加速度是导致[68,137]展向流的主要原因。除理论分析外，Lentink还对模型果蝇翅进行了流场显示和气动力测量，实验结果表明，与旋转阶段相比，在平动阶段，LEV很快脱落且产生较小升力。翅的旋转能够稳定附着的LEV，旋转运动能提高升力，其升力系数为平动升[92]力系数的两倍。在纯旋转运动中，展向流被认为是稳定LEV的主要原因，但一些研究却否定了这个观点。Birch等人研究了绕中心轴旋转的模型果蝇翅。在研究中，为阻止展向流动，安置了沿展向的挡板，利用PIV技术进行了流场分析处理，他们观察到无展向流时LEV仍存在。有挡板时，LEV的环量一般会下降。但是，在翅尖半径处安装柱形壁面，以阻止来自翅尖处的展向流，LEV环量会增加14%。虽然没有观察到LEV涡核内的展向流，但是在LEV下方翅面处，能够观察到翅根流向翅尖的较高气流速度。因为展向流被阻挡时，仍能够观察到LEV，据此推测，翅尖涡稳定了附着的LEV，[138]诱导了下洗气流，将其推向翅表面。Shyy等人研究平板扑翼时也发现了这一现象，[139]在翅尖处，翅尖涡产生了低压区以稳定LEV。Birch等人研究了扑翼拍动时的流场结构，发现扑翼的气动性能与LEV相关，气动力的大小受LEV强度限制。翅尖涡诱导的下洗气流克制了LEV的增长，这预示着下洗气流能够有效抑制LEV的强度和力的大小。他们还利用PIV技术研究了尾迹效应对气动性能的影响，在模型翼根处安装应变片以测量力的大小。每个周期拍动中涡的产生和脱落均影响了下个周期的力。翅从静止开始启动，随着LEV的增长，力也13 北京理工大学博士学位论文增大；LEV达到稳定时，力保持不变。扑翼模型提前翻转，LEV增长也会引起力的增加。翼反转时，LEV和旋转启动涡脱落到尾迹中，初期能提高力，随后力有较小的[88]衰减。他们还研究了雷诺数对LEV和展向流的影响，当Re=120和Re=1400时，扑翼以固定攻角匀速拍动时，有稳定的前缘涡。当Re=1400时，涡核内有展向流，当Re=120时，没有发现展向流。研究结果还表明：通过轴向流的涡量输送并不是低雷诺下稳定LEV的必要条件；不同雷诺数下，由前缘向尾迹中涡量的输送方式可能不[87]同。Devoria等人研究了两个矩形平板（展弦比分别为2和4）的启动旋转。染色剂流场显示结果表明，这两个矩形平板，在LEV内部从翼根至翼尖均有较强的展向流动；[140]较小的展弦比时，LEV的稳定时间更久。Carr等人的PIV技术研究实验表明，平[141]板翼前缘处只有一个LEV。在几个平行弦截面的测量结果验证了LEV通过翅尖涡（TV），与尾缘涡（TEV）相连接，在平板翼表面形成了“马蹄状”的涡。对于展弦比为4的平板，他们观察到LEV在翼展50%处升高，倾斜流入气流方向。类似于Devoria[142]等人的染色剂流场显示结果，Carr等人观察到LEV附着于翼面，翅尖涡诱导的下洗气流稳定了LEV，更重要的是，他们观察到了在LEV内，存在着从翅根至翅尖和从翅尖至翅根的流动。Ansari等人在水箱中研究了展弦比为4的矩形平板翼在雷诺数500和15000之间的纯旋转运动，观察到了附着的LEV，PIV的速度测量结果表明由翼根至翼尖的展向[143]流动，展向流速度约为翼尖速度的80%。Ansari指出在低雷诺数下展向流动更明[144]显，这与Wilkins和Knowles的计算结果是一致的。Wilkins和Knowles二维模拟计算了矩形平板低雷诺数下的LEV特性，计算结果表明LEV与雷诺数和攻角相关。Re<5时，没有观察到LEV。当攻角为45°时，Re=25是LEV不稳定的边界值。在三维情况下，1200.03m时，速度趋于稳定，直至到壁面处，速度急速下降至0。当ah/=6、8、10、12、14时，如图2.6所示，当x值较小（x≈0.015m~0.03m）时，即在靠近实验对象边界处，速度有明显的变化；在靠近实验段壁面处，x值最大的附近处，速度变化剧烈，在x极小的变化范围内，壁面处速度降为0；在x相当长的范围内（0.03m<0.12m时，速度基本不再变化，因此只取0≤≤y0.12m时的速度分析，如图2.10所示。图2.10为实验对象边界处线段CD上的速度分布曲线图。图2.10-1是速度随y的变化曲线图，横坐标为y，纵坐标为速度V，不同颜色的曲线代表不同L下的速度曲线。图2.10-2是速度随L的变化曲线图，不同颜色的曲线代表11不同y下的速度曲线。如图2.10-1所示，不同L值时的速度曲线几乎重合，分别取y=0.004m、0.008m、10.018m、0.060m、0.100m作为特征点（图2-10-1中虚线所示），分析速度随L的变1化。如图2.10-2所示，在相同y下，L变化时，速度V几乎没有变化。这表明，当L11取0.30m、0.45m、0.60m、0.75m和0.90m时，均能够保证流场品质。因此，当L≥0.3m，1实验段长度能完全反应实验对象周围的绕流，流场品质能够得到保证。在风洞设计中，为了充分保证流场品质，将实验段长度取值为L=0.6m。1（1）（2）图2.10实验对象边界处CD上的速度分布曲线2.4.4实验段尺寸的确定利用数值模拟方法计算不同实验段长度L（L=0.45m、0.6m、0.75m、0.9m）和11不同实验段截面边长a（a=0.06m、0.12m、0.18m、0.24m、0.30m、0.36m、0.42m）时，风洞实验段流场的发展规律。通过监测实验段中分截面处的速度分布云图和实验对象边界处的速度分布特征，得到了风洞实验段尺寸对流场变化规律的影响。在经济条件许可的情况下，实验段尺寸越大，风洞内的流场品质越佳。但是考虑到场地、加工和经济性，应该在保证流场品质的前提下尽可能地减小风洞的尺寸。通过分析计算，确定了实验段的最佳尺寸，实验段截面边长a=0.3m，实验段长度L=0.6m。138 北京理工大学博士学位论文2.5扩散段尺寸对流场品质规律的影响2.5.1扩散段尺寸、模型及边界条件扩散段的作用就是把气流的动能变成压力能。由于风洞损失与流速的三次方成正比，所以经过实验段的气流应尽量降低速度，把动能转变为压力能。但是，减速必然伴随损失，主要就是摩擦损失和扩散损失。摩擦损失就是相当于气流流过管道遇到摩擦而产生的损失，因此扩散段长度增加，摩擦损失就会增加。扩散损失就是扩散段中的气流速度减小而压力增加，因而流动处于逆压状态，附面层内气流受到阻碍后，厚度增加很快，相当于壁面阻力增加，损失能量。由于设计的低速风洞主要用于流场显示实验，能效高低不是研究问题的重点，因此，在确定扩散段尺寸时，应重点考虑流场品质的问题。扩散段垂直于风洞中心轴线的截面均为正方形，扩散段入口的截面边长即为实验段的截面边长（a=0.3m）。扩散段出口的直径b由动力段风扇的直径确定，b=0.4m。因此主要通过流场品质的高低来确定扩散段的长度L。如图2.11所示，扩散段的长度为L，扩散角为β。扩散段22的长度L降低时，扩散角β增大，这时风洞边缘处容易发生气流分离，对流场的品质2有很大的影响。当扩散段长度L增大时，扩散角β降低，可以提升流场品质，但是风2洞的长度增加会提高风洞的造价成本。因此，需要确定合适的扩散段长度L。2图2.11扩散段示意图计算当扩散段长度L分别为0.2m、0.4m、0.6m、0.8m和1.0m时的流场。如图2.122所示，实验段入口给定速度入口边界条件，来流速度为2m/s，扩散段出口为压力出口，压力值为大气压101325Pa，风洞壁面采用无滑移固壁，采用标准k−ε湍流模型求解。为节省计算量采用对称边界条件，即只计算风洞的四分之一区域，网格均采用六面体网格，如图2.13所示。39 北京理工大学博士学位论文图2.12边界条件示意图图2.13计算域网格示意图如图2.14所示，选取截面p（实验段出口，扩散段入口）和截面q（扩散段出口）为参考截面。如图2.14-1所示，在截面p上，选取截面中心O到壁面端点E的线段OE作为衡量流场品质的参考线段。在扩散段长度S不同时，观测线段OE上的速度分布。同理，如图2.14-2所示，在截面q上，选取O"E"作为参考线段。（1）（2）图2.14截面p和截面q示意图2.5.2扩散段长度对扩散段出口流场品质规律的影响扩散段长度的改变会影响风洞内的气流品质，特别是在扩散段有可能出现气流分离现象。为此，选取扩散段出口截面，对比不同扩散长度时出口截面处的速度分布。如图2.15所示，随着L增大，风洞边缘处的气流分离逐渐减弱。当L=0.2m时，在扩22散段边缘处有明显的气流分离现象，当L>0.4m时，边缘处气流分离不明显。2为此，将截面处斜对角线O"E"上的速度分布进行定量分析。如图2.16-1所示，当L=0.2m时，沿着z坐标方向，在靠近壁面处的速度是不稳定的，速度曲线有振荡，22说明此处的气流有分离，而当L>0.4m时，沿着z方向的坐标的速度变化是平稳的，22即在风洞中心速度保持平稳，在靠近壁面处，速度会平稳地降至0。40 北京理工大学博士学位论文L=0.2mL=0.4m22L=0.6mL=0.8mL=1m222图2.15扩散段出口截面q的速度云图为衡量不同S的流场品质，取z的多个特征点，如图2.16-2所示，分别取z=0、220.100m、0.133m、0.180m时，监测速度V随L的变化曲线。在扩散段出口风洞中心2处（z=0），随着L增大，速度V在降低。当L=0.2m时，V>1.4m/s；当0.4m≤≤L1m2222时，速度V在1.2m~1.3m范围内，变化不大。z=0.1m时，V随L的变化趋势与z=0222时基本相同。z=0.133m时，不同L下的速度V基本相同，均在1.1m/s附近。当22z>0.133m时，在靠近风洞壁面处，速度V随z增加快速下降。z=0.18m时，速度V222随L增大而增大，这表明L越大，在壁面处速度的变化越平缓，边界处越不易产生气22流分离。z=0.18m，L≥0.6m时，速度V基本保持平稳。22综上所述，当L≥0.6m时，在扩散段出口截面处，风洞边界气流不会产生分离，2可以保持良好的流场品质。41 北京理工大学博士学位论文（1）（2）图2.16线段O’E’的速度曲线图2.5.3扩散段长度对扩散段入口流场品质规律的影响对扩散段入口截面处（即图2.12中的截面p）进行分析，图2.17为不同L下的截面2p的速度分布云图。如图2.17所示，速度范围为0~2m/s。当L=0.2m，风洞中心有较2大范围的低速区域；当L=0.4m时，低速区域的范围减小，当L≥0.6m时，低速区22域消失。对于风洞实验段的截面，风洞内大部分区域保持一致速度，只在壁面附近处有速度的变化，是最理想的情况，此时的流场品质应是最佳。因此，由图2.17，当L≥0.6m时，风洞的流场品质能够满足要求。2L=0.2mL=0.4m2242 北京理工大学博士学位论文L=0.6mL=0.8mL=1m222图2.17扩散段出口的速度云图图2.18为在不同扩散段长度L下，线段OE上的速度V随坐标z的变化曲线图。21如图2.18-1所示，在z最大值附近，也就是风洞壁面处，速度V变化范围大，急剧下1降为0。因此选择0≤0.6m时，V>2m/s。当L=0.2m和2222L=0.4m时，在0≤0.12m范2111围内，z值增大，速度V随之降低。当L=0.6m、0.8m、1.0m时，随着z值增大，速121度V会减小。在靠近风洞壁面处，应保证气流不会发生分离，才能够保证较好的流场品质，因此选取靠近壁面处的特征点，取z=0.12m、0.13m、0.14m时进行分析。如图2.18-3所1示，当z=0.12m和z=0.13m时，L≤0.4m时，速度V随S增大而增大，当L>0.4m时，1122速度V随L增大而降低。当z=0.14m时，速度V随L增大而降低。z=0.14m处的速2121度越低，在壁面处0.14m≤kk2LLm(1−m)33五次方曲线计算公式如式2.3所示。45 北京理工大学博士学位论文345kkkx"=110−+15−6(0.5c0.5a)0.5a−+（2.3）LLL3332.6.2.2收缩曲线对流场中分截面流场品质规律的影响收缩曲线影响着风洞内部的流场分布。不同收缩曲线下会产生不同的流场分布，图2.22为不同收缩曲线下风洞中分截面的速度云图。图2.22-1、图2.22-2、图2.22-3分别为收缩曲线选用维多辛斯基曲线、双三次曲线和五次方曲线时风洞中分截面的速度云图。如图2.22-1所示，收缩曲线采用维多辛斯基曲线时，在收缩段内，速度沿着风洞轴向方向均匀增加。如图2.22-2和图2.22-3所示，收缩曲线采用双三次曲线和五次方曲线时，在收缩段内，靠近收缩段入口处，速度沿轴向的变化较不均匀，在收缩段出口附近，速度沿轴向的变化较为均匀。在实验段内，三种曲线下的风速均能达到实验的稳定风速2m/s。在扩散段内，风速的变化也较为相似。（1）（2）（3）图2.22不同收缩段曲线下风洞中分截面速度云图2.6.2.3收缩曲线对收缩段出口截面流场品质规律的影响图2.23为不同收缩曲线下收缩段出口截面的速度云图。图2.23-1、图2.23-2和图2.23-3分别为收缩曲线选用维多辛斯基曲线、双三次曲线和五次方曲线时风洞中分截面的速度云图。如图2.23所示，在截面中心处的速度最高，越靠近壁面，速度越低。46 北京理工大学博士学位论文（1）（2）（3）图2.23不同收缩段曲线下收缩段出口截面速度云图图2.24为不同收缩曲线下收缩段出口截面边长的速度曲线图。如图2.24-1所示，当z<0.14m时，速度基本保持在2m/s附近，当z>0.14m时，由于靠近风洞壁面处，速33度急剧下降。图2.24-2中，z的范围是0≤≤z0.14m，在这个区域内，速度的变化范33围较小。当z=0时，双三次曲线的速度最高，维多辛斯基曲线速度最低。实验中的3风速为2m/s，当速度小于2m/s时，速度不稳定，开始下降。可以看出，维多辛斯基曲线到达2m/s时，z=0.13m；双三次曲线到达2m/s时，z=0.125m；五次方曲线到达332m/s时，z=0.12m。这表明，维多辛斯基曲线下，截面处的稳定速度范围区域最大，3流场品质最好。（1）（2）图2.24不同收缩段曲线下收缩段出口截面边长的速度分布2.6.2.4收缩曲线的确定分别数值模拟计算了维多辛斯基曲线、双三次曲线和五次方曲线这三种收缩曲线对风洞流场品质的影响，通过对比风洞中分截面及收缩出口截面的速度分布，可以得到使用维多辛斯基曲线对风洞的流场品质最佳。因此采用维多辛斯基曲线。47 北京理工大学博士学位论文2.6.3收缩比对流场品质规律的影响2.6.3.1收缩比对流场中分截面流场品质规律的影响22收缩比λ，是指收缩段进口面积与出口面积之比，即λ=ca/。收缩段出口边长为实验段截面的边长a，则收缩段进口截面边长c的尺寸由收缩比λ确定。一般情况下，收缩段的长度L与收缩入口边长c有关，收缩比越大，收缩段的长度就越长，取3Lc=。分别计算收缩比λ=4、5、6、7、8时的流场。3λ=4λ=5λ=6λ=7λ=8图2.25不同收缩段曲线下风洞中分截面速度云图图2.25为不同收缩比λ下风洞中分截面的速度云图，如图2.25所示，当λ=4、5、6、7、8时，实验段内的风速均在2m/s附近处，均能达到平稳的实验风速。在收缩段出口与实验段入口相接附近处，收缩比越大，风速变化越平缓，这代表流场品质越好。图2.26为实验段中心轴线O""O的速度曲线图。图2.26-1为不同收缩比λ下，速度随坐标y"的变化曲线，不同收缩比下，图中不同颜色曲线代表不同收缩比λ。如图2.26-1所示，当λ=4、5、6、7时，随着λ值增大，速度曲线上移，这表明在相同坐标y"下，λ值越大，速度越高。随着λ值增大，实验段入口和出口处的速度差值也在降低，这表明在这个范围内，λ值越大，流场品质越好。λ=7和λ=8速度曲线几乎重合。取y"=0、0.300m、0.465m、0.600m作为参考点，图2.26-2是参考点处速度随λ的变化曲线。如图2.26-2所示，不同颜色曲线代表不同坐标y"，在任一参考点处，速度均随λ增大而增大，当λ≥7时，速度基本保持不变，这表明此时已达到较佳的流场品质。48 北京理工大学博士学位论文（1）（2）图2.26线段O""O的速度曲线2.6.3.2收缩比对收缩段出口截面流场品质规律的影响图2.27为不同λ下的出口截面的速度云图。由图2.27，当λ=4时，风洞中心附近处的速度小于2m/s，由中心到壁面处，速度先增大后减小。当λ≥5时，在中心大部分区域，速度均能保持较为稳定的速度。这表明当λ≥5时，流场品质较好。图2.28为参考线段O""E""的速度曲线。图2.28-1是在不同λ下，速度随z3（z=0~0.15m）的变化曲线，如图2.28-1所示，在风洞中心的大部分区域，速度在32m/s附近，在靠近壁面处，速度急速下降为0。除去壁面附近的点，取z=0~0.14m，3图2.28-2为速度随z的变化曲线。如图2.28-2所示，在z=0处，λ=4时，速度V<2m/s；33λ≥5时，速度V>2m/s。这也是图2.27中λ=4时风洞中心附近处出现速度间断的原因。当λ=4、5、6时，在风洞壁面边界处（z>0.12m）速度梯度较大。当λ=7、83时，速度曲线几乎重合，且在风洞壁面处的速度变化较为光滑。取z=0.12m、0.13m、30.14m作为特征点进行分析，如图2.28-3所示。当z值一定时，速度V随收缩比λ增大3而降低，这表明，λ越高，壁面处速度梯度越小，速度变化越平缓，流场品质越好。当λ≥7时，速度随z的变化很小。因此，可以认为当λ≥7时，风洞能够达到较好的3流场品质。综上所述，考虑到流场品质的高低和风洞的制作成本，初步选择收缩比λ=7。49 北京理工大学博士学位论文λ=4λ=5λ=6λ=7λ=8图2.27不同λ下收缩段出口截面的速度云图（1）（2）（3）图2.28线段O""E""上的速度曲线2.6.4收缩段尺寸的确定通过对比维多辛斯基曲线、双三次曲线、五次方曲线的流场品质，选择维多辛斯基曲线作为收缩段曲线。根据收缩比λ=7，收缩段出口截面边长a=0.3m，可计算收缩段入口截面边长为0.79m，考虑到制作方便，对收缩段入口截面边长取值c=0.8m。22由a和c可最终确定收缩比λ=ca/=7.11。50 北京理工大学博士学位论文2.7风洞整体结构参数及流场品质验证综上所述，可以确定风洞实验段、扩散段、收缩段的尺寸。风洞其他结构部分的尺寸应根据具体情况给出，具体见第三章实验系统的介绍。风洞的整体尺寸如表2.1所示。根据风洞的尺寸不同，风洞洞体可以选择木结构、钢筋混凝土结构、钢结构或混合结构。这几种材料各有其优缺点，但就风洞气流品质来讲，这几种材料之间是没有差别的。木结构质量轻，容易加工和改造，吸收振动和噪声的效果比较好，价格也相对便宜，因此，制作的风洞洞体完全采用木结构。制作完工的小型低速风洞照片如图2.29所示。表2.1风洞的结构数据表喇叭口稳定段收缩段实验段扩散段动力段入口截面边长/(m)10.80.80.30.30.4出口截面边长/(m)0.80.80.30.30.40.42长度/(m)0.10.60.80.60.60.5图2.29风洞实体图将蜜蜂置于风洞中进行流场显示实验，如图2.30所示。烟线经过蜜蜂身体中间，烟线在未经过蜜蜂之前，烟线保持层流状态。烟线经过蜜蜂，在蜜蜂身后产生绕流的尾迹。未经过蜜蜂的烟线一直保持层流状态，这说明风洞具有良好的流场品质。图2.30中可以看到足够大的尾迹范围，这表明所设计的风洞用于小型低速昆虫的飞行是可靠的。因此，所设计的风洞完全可以用于小型昆虫的流场显示实验。51 北京理工大学博士学位论文图2.30蜜蜂振翅的烟线图2.8本章小结为专门研究蜜蜂等小型昆虫飞行的流场特征，对蜜蜂飞行进行流场显示实验，本章设计了小型低速风速的整体结构，风洞为开路闭口式，风洞的整体结构包括喇叭口、稳定段、收缩段、实验段、扩散段和动力段。为保证风洞内实验气流的流场品质，将气流速度作为评价指标，利用数值模拟计算方法，分别研究了实验段、扩散段和收缩段尺寸、以及收缩曲线线型对流场品质的影响，确定了实验段、扩散段和收缩段尺寸，并将维多辛斯基曲线作为收缩曲线线型，确定了风洞的整体结构参数。对蜜蜂的流场显示实验结果表明，所设计的结构及参数均能够保证风洞的流场品质。52 北京理工大学博士学位论文第三章流场显示实验系统开发3.1引言在流体力学领域，19世纪Reynolds在水平圆管内利用直接注入颜色水方法进行[188]了层流、湍流及其转捩实验，这个实验被追溯为流场显示技术的开端。流场显示实验在科学发展上起了很大的作用，是流体力学和空气动力学的一种重要实验手段。[189]例如湍流、涡流、边界层等现象都是首先被流场显示实验观察到的。流动显示实验的烟线图对提示发现流体力学的基本现象，建立理论模型均有较大的推动作用。随着流体力学的发展，流场显示技术也不断地得到发展。流场显示技术目前已被广泛应用于流体力学、空气动力学、爆炸力学、等离子体物理、燃烧学、传热学、大气物理、空间技术和化学反应工程等一系列领域中。特别是在风洞实验中，流场显示实验[190]是研究涡等现象可靠而有效的方法。在流体力学中，研究流场特征时，流场内的流体可能静止的，也可能是随时间变化的。通过跟踪流体随时间变化的规律，能够分析流场特征。但是，对于多数流体，例如水和空气，均是透明的介质，人的肉眼无法观察到这些透明流体的变化情况。在力求不改变流体力学性质的前提下，用图像显示流动现象的方法。在水洞或风洞中常用这类方法来显示流动分离、尾流、旋涡、边界层和激波等流动现象。流动显示方法可直接用于风洞流场校测和各类空气动力实验，为空气动力计算提供可靠的流动模[191]型，还能发现一些新的流动现象。在流体力学实验中，实现流场显示有三种方法：（1）表面流动显示：这种方法就是把示踪用的涂料涂布到物面上，或浸泡在物面上。当流体流过时，涂层便可显示出流动状态。这种方法主要用来显示物体表面的流动，适合于气体或液体的定常流，对于非定常流只能观察到时间的平均流态。这种方法主要是用于研究壁面处的边界层，特别是分离气流和附着气流。（2）示踪粒子法：在不影响流场结构前提下，需要向流场中引入有色示踪粒子，示踪粒子随流体运动。通过追踪粒子轨迹，不仅可以看到流动过程，而且可以用PIV方法通过测量粒子速度，求得流场的速度场，进行定量分析研究。（3）光学测量方法：通过改变光学折射率的方法进行流场显示，包括阴影、纹影、干涉等技术。光学方法的特点是不扰动流场，属于非接触性测量，适于可压缩流53 北京理工大学博士学位论文[192]动显示。烟线法就是示踪粒子法的应用。在流场中释放烟线，通过观察烟线绕流实验对象[193]边界的气流变化，分析流场的局部特征规律。研究昆虫飞行流场的烟线实验就是在上游释放烟线，烟线在上游保持层流状态，经过昆虫身体或翅时，层流状态被破坏，烟线可以有效地显示昆虫附近的绕流流场，由于昆虫的飞行是非定常的，在不同时间下烟迹也不相同，但是当烟线经过翅时，由于翅是周期性的拍动过程，烟线也呈周期性变化规律，这能有效地研究昆虫振翅飞行的机理。在烟线显示实验中，首先要保证风洞内的流场品质。第二章已经利用CFD数值模拟方法，根据流场品质的变化规律，确定了风洞结构和尺寸参数。但是在CFD计算中，给定的入口条件是假设均匀的。大气中的气流湍流度高，平稳性差，为保证进入风洞的气流是均匀平稳的，风洞内还需要设计整流装置，衰减气流中旋涡的尺度，降低湍流度。其次，保证良好的发烟质量。烟线的质量对流场显示实验至关重要。蜜蜂翅的尺寸很小，翅长约10mm左右，因此既要保证烟线的浓度，也要保证合适的烟线粗细程度。本章节中设计调试了发烟装置，通过对比发烟时间，发烟效果和发烟浓度，确定了电热丝的材质和形状，选取了烟油的品种。在实验中，为了确保烟线相对于翅的位置，需要两台高速摄像机拍摄两个角度，为保证高速摄像机的高帧率拍摄，需要保证充足的光线，在实验中，选取了LED光源，搭建了完整的实验平台。虽然在一些学者研究中，用悬吊方法研究昆虫的飞行，例如Willmott对悬吊的鹰[75][132][181]蛾，Srygley和Thomas研究了蝴蝶和蜻蜓的悬吊飞行和自由飞行。但是，[182]这些昆虫的尺寸较大，体积约为蜜蜂尺寸的四五倍。即使是Bomphrey研究的熊蜂，也是蜜蜂的三倍大小。发烟的时间短暂，如果蜜蜂自由飞行，高速摄像机难以抓拍到翅恰好经过烟截面的视频。蜜蜂的胸腔为硬壳，与蛾类和蝶类柔软的身体不同，不适合悬吊飞行。本文流场显示实验的重点在于研究蜜蜂翅的拍动对流场的影响，需要使得烟截面恰好经过翅上，最好的方法是将蜜蜂身体固定，这样便于调整烟线相对翅的位置，使烟经过翅上。本章根据搭建的实验平台和实验方法，研究了圆柱绕流和蜜蜂振翅的流场显示实验。实验结果表明，开发的流场显示实验系统能够用于研究蜜蜂振翅的流场特征。54 北京理工大学博士学位论文3.2烟线显示的意义描述流体运动有两种方法，拉格朗日方法和欧拉法。拉格朗日方法着眼于流体质点，目的是描述出每个流体质点的运动规律，即质点位置随时间变化的规律。欧拉方法着眼于空间点，目的是在空间每一点上描述出流体运动随时间的变化状况。迹线（Pathline）与拉格朗日方法相联系。迹线表示一段时间内同一个质点的经历，是流体质点的运动轨迹。同一时刻，质点的微元位移总是和它的速度同方向。因此，定常流动中，流线与迹线是重合的，即质点沿着流线运动。式3.1为迹线运动方程。ddrrxyzt=(,,,)（3.1）流线（Streamline）与欧拉方法相联系。流线是指示某一时刻流场中各点速度方向的假想曲线，曲线上任意一点的切线方向与该点速度方向一致。流线方程如式3.2所示。dddrV×=0（3.2）流向不能突然折转，因为流点有速度，就有惯性，因此流线只能是连续光滑的曲线。在流场中，流线是不相交的，流线的疏密程度反映了该时刻流场中各点的速度变化。在定常流动中，流线是稳定不变的曲线。脉线（Smokeline）是把相继经过流场中同一空间点的流体质点在某瞬时顺序连接起来得到的一条线。在流场固定点连续释放染色剂，每一瞬时，从固定点出发的由染[194]色剂组成的脉络线就是脉线（也称染色线）。脉线是瞬时线，但脉线与先后从出发点流出的众多流体质点的经历又密切相关，脉线兼有流线和迹线的特征。在定常流中，脉线、流线、迹线重合，而且形状都保持不变，因此常用脉线来表示流线。但如果流动是非定常的，脉线与流线不重合，就不能简单地把脉线当做流线。脉线很容易在实验中观察，因此脉线常被用于做流动显示[195]实验。根据亥姆霍兹旋涡定理，组成涡面的流体质点永远在涡面上，所以如果将这些染色点设置在绕流物体的分离点或分离线处，则那些染过色的带有涡量的流体微团都将[196]随自由剪切层进入旋涡中，从而使流场中的旋涡被显示出来。所以，脉线（染色线）流动显示技术是显示旋涡流动现象的有力工具。烟风洞就是利用了脉线流动显示技术，通过加热电热丝上的发烟油，产生可显示的有色烟线。电热丝两端加上可调节的电压，通过调整电压的大小调节电热丝的温度，55 北京理工大学博士学位论文当电热丝达到一定温度，附在电热丝上的油滴会气化为蒸汽。油蒸汽随气流离开电热丝，在常温的气流中又立即凝结成细小的油雾滴，这些雾滴会显示可见的绕流形态。研究小型昆虫飞行的空气动力学原理，难以直接测量翅上的压力和速度，虽然机械扑翼作为一种有效的研究方法，但是它只能反应刚体的运动，难以模拟生态翅的扭曲和弯曲，也就不能真实反映昆虫振翅时的流场特征。早在20世纪70年代，一些学者们提出了解释昆虫的高升力需要利用非定常空气动力学原理，非定常运动是昆虫飞行最重要的特征。在昆虫振翅飞行时，身体和翅周围均会产生较强的涡流，翅的运动会产生涡，涡也会对昆虫产生诱导速度和诱导升力。来流将烟引入到有涡量的区域，烟线会卷起，形成旋涡。同涡量一样，烟线符合伽利略变换，因此，烟线的形式是独一无二的，与参考系无关，这适用于研究昆虫复杂的飞行运动。利用烟线法研究昆虫附近产生的涡流以及涡流的演变，有助于分析昆虫在非定常运动下的空气动力学特性，是研究真实昆虫周围流场最直接有效的手段。利用烟线法显示流场时，烟线是脉线，与流线有本质的区别。在非定常运动中，烟线不同于流线。在翅的前缘处，可以近似看做定常流动，绕过前缘的烟线与流场中的流线是重合的，这与在来流中的机翼相似，来流绕过机翼的前缘，气流产生分离。翅在快速旋转（例如在上下拍动交替时）或加速时，脉线与流线是分离的，烟线不能代表当前流体的运动，而是运动在空间处时间历程的叠加。在每个单张烟线流场图中，脉线不代表真实的流线。高速摄像机拍摄了连续多帧的图片，由连续多张烟线流场图中，能够观测到烟线随时间历程的演变，分析流场时，脉线与流线的区别可以适当地降低或忽略。在单张烟线流场图中，烟线的运动可能影响瞬时的流场，烟线显示的某些离散特性（螺旋或者涡），可能表示的是时间叠加的流场，而非当前流场。因此，利用高速摄像机拍摄流场的连续变化是非常必要的，在分析流场结构时，需要细致地分析每一帧图片中烟线的变化，将多张连续帧图片组合分析，才能得到流场的变化规律。在烟线流场显示实验中，涡的发展是流场中最主要的特征。烟线是由油蒸汽粒子组成的，随着涡发展的粒子有可能被留在后面，不能标记出流场中重要的涡。但是，涡一般都是经过一段时间累积形成，这段时间的长度远高于昆虫飞行的频率。在实验中，烟应该具有足够的浓度用于流场显示实验，此外，烟线产生的时间应当远远高于蜜蜂的振翅周期，这样才能够保证能够利用烟线法对昆虫飞行的流场进行研究。在风洞中，在实验段上游释放烟线，如果没有干扰，烟线会一直保持层流状态，56 北京理工大学博士学位论文当烟线被来流输送至蜜蜂附近处，经过蜜蜂及翅时，层流状态被破坏，产生分离绕流。随着振翅运动及来流的速度，烟线被输送到下游，形成由于振翅产生的尾迹流场结构。由于蜜蜂振翅运动是周期性的，烟线在绕过翅形成的尾迹结构是规律性的，因此可以被用作分析流场。3.3整流装置3.3.1喇叭口和稳定段在稳定段前方，有一段光滑唇口，使风洞外的气流能从各个方向光滑地进入风洞。将此光滑唇口设计为喇叭形入口，位于稳定段的前部，保证气流平滑进入风洞，防止在风洞入口边缘形成旋涡，改善进入稳定段气流的流动。喇叭口弧形为1/4圆弧，半径取20cm。在风洞收缩段的上游，通常与等截面管道的稳定段相连接。因此稳定段的边长等于收缩段入口的边长。稳定段的功能在于使外面进来的紊乱不均匀的气流稳定下来，衰减旋涡，使气流速度平行于风洞主轴。如果来流不均匀，气流经过收缩段后也是不均匀的，实验段的气流就不均匀，无法保证风洞的流场品质。因此在收缩段的上游，必须要经过稳定段，使气流均匀，从而保证风洞的流场品质。稳定段内如果不设置任何整流装置，就必须有足够的长度，使气流在流动过程中有足够时间调整运动方向、速度分布并衰减紊流度。从加工和经济角度来看，应当尽可能减小风洞的尺寸，降低成本。因此，一般在稳定段内都装有整流装置（即蜂窝器和阻尼网）。这样既能够降低稳定段长度，又能提高气流均匀性和稳定性，保证流场品质。稳定段的长度首先要保证安装整流装置，其次还要有一段长度，使气流经过整流装置后逐渐稳定下来并衰减残存的小旋涡。对于收缩比λ≤5的风洞，稳定段长度一般为mc=~1.5c，对于收缩比λ>5的风洞，稳定段长度m=0.5~cc。由第二章已确定风洞的收缩比λ=7.1，给定稳定段长度mc=0.75，由于截面边长为c=0.8m，则m=0.6m。3.3.2蜂窝器和阻尼网蜂窝器常由方形、圆形、六角形等截面的小格子组成阵列，形同蜂窝。蜂窝器的作用在于导直气流，使其平行于风洞轴线，将气流中的大尺度旋涡分割成小旋涡，有利于加快旋涡的衰减。同时，蜂窝管道对气流的摩擦，会逐渐减弱小旋涡的能量，还57 北京理工大学博士学位论文有利于改善气流的速度分布。在开路风洞中，气流由四面八方流入风洞，蜂窝器是不可缺少的。气流经过蜂窝器后，旋涡的尺度必将小于蜂窝口径，气流紊流度下降。同时，蜂窝器限制或减小了气流的横侧运动，使气流方向更接近于平行风洞轴线。蜂窝格子的形状影响风洞的损失系数，但是对于小型低速风洞，可以忽略损失系数的影响。主要考虑蜂窝格子对气流品质的影响。一般来说，六角形蜂窝格子整流的效果最好。影响蜂窝器性能的主要参数是蜂窝长度和口径。长度越长，整流效果越好。[183]口径越小，降低紊流度的效果越显著。选用口径6.5mm、长度为210mm的六角形格子蜂窝器。阻尼网的作用是可以进一步破碎蜂窝器后的气流旋涡，降低气流的湍流度。阻尼网网眼及网线的直径都很小，作用与蜂窝器类似，将较大的旋涡分割成小旋涡。阻尼网越细，网的层数越多，整流效果越好。阻尼网与蜂窝器的基本区别在于阻尼网不能对气流起导向作用。阻尼网一般装在蜂窝器之后，与蜂窝器一起使用时整流效果非常[184]显著。在所设计的风洞中，在稳定段与动力段均加入了蜂窝器和阻尼网，既能将风洞内的气流进行导直，使气流平行于风洞中轴，又能降低旋涡的尺度，以降低风洞内的紊流度，从而改善风洞内的气流品质。图3.1风洞的蜂窝器示意图3.3.3动力段动力段主要用来产生人造风并调节风速，以保证气流以一定的速度平稳地在风洞中流动。在动力段中，一般在风扇前使用反扭导流片，用于消除风扇所造成的旋流，从而改善气流的状态，提高流场品质。在所设计的风洞中，利用导流蜂窝器替换反扭导流片，作用是将气流平行导出，防止因风扇旋转引起的气流周向流动和离心运动影响到扩散段进而对实验段造成干扰。通过调节电机的功率，控制风扇的转速，从而可以控制调整实验段气流的速度。装在风扇前的导流片和装在风扇后的止旋片都是用于58 北京理工大学博士学位论文消除风扇所造成的旋流，从而改善气流的状态，提高流场品质。该风洞的动力段直径0.42m，长度0.5m。3.4发烟装置发烟装置用于产生烟线，是流场显示实验的重要组成部分。通过在电阻丝上涂抹油质，电阻丝通电生热，上面的小油珠受热后变成蒸汽，随气流离开电阻丝后又立即凝结成油雾滴，因此产生烟雾。发烟电阻丝宜采用电阻率大、强度高、抗氧化性好的材料，通常选用镍铬或铁铬铝。根据实验段截面宽度，电阻丝取为长度35cm左右的镍铬（Cr20Ni80）丝，垂直安装在实验段的入口中间位置。电阻丝可以选取不同的直径和截面形状，从而得到不同的阻值和发烟效果。选用薄片状电阻丝和圆柱状电阻丝进行发烟实验，对比两者的发烟效果。薄片状电阻丝由于单位长度表面积较大，可以附着更多的油质，加热后产生的烟浓度也较高，但是烟雾连成均匀的烟片，没有明显的烟线结构（如图3.6所示）。柱状电阻丝由于相对表面积小，附着的油较少，产生的烟雾较淡。但由于油质与电阻丝之间不浸润以及油的表面张力，其表面上的涂抹的油每隔一小段汇聚成一颗小油滴，加热时从每颗油滴处拉出一条清晰的烟线，所以能够观察单根烟线的流动轨迹（如图3.7所示）。柱状电阻丝适合用来观察流场的细微结构，而使用片状电热丝则可以观察流场的整体扰动情况。对于昆虫振翅实验，为记录翅周围的精细流动特征，因而选择柱状电阻丝。在强度允许的条件下尽量选用较细的阻丝，经过实验对比烟线的密度和清晰度，选取直径0.1mm的镍铬丝Cr20Ni80做电阻丝，每米阻值138.8欧。烟油以发烟白、油雾颗粒小、无毒、无腐蚀性为好。为了使电阻丝上能黏有足够的油，通常选用粘度较大的。实验对比了甘油、橄榄油和婴儿油的发烟效果，强生婴儿油产生的烟线最为清晰，因此被选为实验烟油。烟线的浓度主要由改变电阻丝两端的电压来调节。电源使用220V的可调压交流变压器，调压范围为0～220V。实验时需要按照不同的风速调节合适的电压值，风速在1m/s～2.5m/s时，相应的电压在20V～40V之间。相同风速下，电压值越高，发烟浓度越大，烟线也越清晰。但发烟时间会越短，而且电阻丝越热，会导致附近的对流扰动越严重，烟线也会越乱。反之降低电压，则烟线浓度减小，清晰度变差，但发烟时间会延长且烟线趋于稳定。所以使用烟线法必须在清晰度和稳定性两者间进行平衡。图3.2-1是电压45V的烟线图，烟线浓度高但却很乱；图3.2-2是电压20V的烟线图，烟线稳定但不够清晰。59 北京理工大学博士学位论文（1）45V（2）20V图3.2不同电压值的烟线对比图在烟线流动显示实验中，从发烟到烟的扩散、消失的时间很短，为此需要设计专门的时间电路以确保发烟、延时、闪光和照相的协调进行。通常，从通电到发烟之间的时间约为1ms，闪光照相延时时间为10ms。对光源的基本要求是应具有足够强而均匀的照度，以便得到烟线较强的散射光。3.5建立实验平台采用奥林巴斯i-speed3高速摄像机，帧率选择10000fps，曝光时间95.938μs，图片分辨率528×396pixel。高速摄像对光源有极高的要求，蜜蜂的振翅频率每秒高达250Hz，所以摄像照明需要连续或者高频的高亮度光源。连续光源可以选择碘钨灯或LED灯，高频光源可以选择无极灯。碘钨灯属于热光源，色温偏低（小于3000K）且热功率大（一个灯管就有2000W），蜜蜂常被烘烤致死。LED灯和无极灯都属于冷光源，亮度高发热少（功率只有100W左右），色温（5300K~6000K）接近正午阳光，适合用来拍摄活体昆虫。实验中选用了两台LED灯，置于风洞实验段的顶部，透过天窗进行照明。实验时高速摄像机对烟雾面进行对焦，调整蜜蜂的位置使得拍动的翅刚好划过烟雾。对于弦长只有4mm左右的蜜蜂翅来说，实验选取尼康105mm微距镜头或者附加延伸管的85mm镜头作为高速摄像机镜头。光圈和拍摄距离都会影响景深。光圈决定进光量，光圈越大画面感光越充分，因此大光圈可以改善高速摄像经常发生的曝光不足问题。但是较大的光圈会导致较浅的景深，画面中清晰的部分会很浅。为拍摄清楚蜜蜂翅周围烟线的绕流细节，高速摄像机会离蜜蜂较近，以获得较高的放大率。然而拍摄距离越近，景深也会越浅（甚至只有两三毫米）。拍摄时需要调节合适的镜头光圈和拍摄距离，来实现最理想的放大率和景深。在烟线显示实验中，要求能够拍摄清楚翅附近的烟截面绕流，也要能够确定烟线60 北京理工大学博士学位论文相对翅的位置。这就需要两台高速摄像机拍摄两个视角，一台摄像机置于风洞顶部（正上方），用于拍摄烟线的位置，另一台置于风洞的正前方（水平放置），拍摄烟截面的流动显示变化。两台LED灯分别置于风洞实验段的对角方向，实验布局示意图如图3.3-1所示，实体图如图3.3-2和图3.3-3所示。（1）（2）（3）图3.3实验布局图在观察蜜蜂周期振翅的流场实验中，一般将烟线置于翅0.5R（R是翅长）处，即当翅拍动到水平位置时，烟线位于翅根和翅尖的中间位置，如图3.4所示。在分析展向流和前缘涡三维流场结构时，需令烟线流经不同的翅展位置（实验中使烟线分别位于0.25R、0.5R和0.75R处），观察并分析烟线的绕流特征，如图3.4所示。图3.4翅不同展向截面的示意图61 北京理工大学博士学位论文3.6准备实验对象首先进行了不同雷诺数下的圆柱绕流实验。选取了4种不同直径（分别为0.76cm、1cm、2cm和3cm）的圆柱。实验使用片状电阻丝发烟，将圆柱固定在风洞中进行绕流实验，观察形成的卡门涡街。主要工作是蜜蜂振翅的流动显示实验。蜜蜂是在北京植物园，通过人工收集方式得到。捕捉到的蜜蜂，装进带有呼吸孔并且事先添加了饲料的小瓶中，然后放进恒温箱里。返回实验室后立刻进行实验，以保证蜜蜂的活性。在实验过程时，室内温度在20℃～25℃，湿度在40%～50%。为了能够拍摄蜜蜂振翅时周围的流场，烟线必须准确流经蜜蜂身体和翅，因此我们采取固定蜜蜂身体的方式，也就是系飞的办法。研究者使用过各种不同的固定方法，[197]例如Willmott等使用细钢管和尼龙绳从下方套住鹰蛾颈部，Kitagawa等利用塑料管和细绳从上方套住甲虫颈部。我们也用类似的方法进行过实验，利用细塑料管和尼龙绳分别从上方和下方套住蜜蜂的颈部，然后置于风洞中，如图3.5-1和图3.5-2所示。但这种方法对小昆虫来说效果不好，小昆虫的头颈吃力较小，而套绳力度直接影响固定效果。使用套索系住蜜蜂的颈部，太紧则蜜蜂失去振翅能力，太松振翅时身体会发生转动。而且颈部的细绳导致蜜蜂的异物感强烈，经常会用脚爪去挠细绳，影响实验进行。（1）从上方套住颈部（2）从下方套住颈部（3）牙蜡黏接背部图3.5固定蜜蜂方式的对比图[198]最后使用牙蜡固定蜜蜂。首先在钢针上粘取少许牙蜡，接着将蜜蜂封入塑料袋。然后立即向袋中吹入纯氮气排出空气，很快（不过两三秒）就能使蜜蜂晕过去。立刻取出塑料袋中已麻醉的蜜蜂，同时用打火机加热融化牙蜡，并迅速将钢针粘在蜜蜂背部，如图3.5-3所示。麻醉的蜜蜂恢复活性，需要大概10秒左右，因此整套动作要求快速娴熟。使用牙蜡固定蜜蜂身体比较稳固，效果很好。而且牙蜡的熔点很低，不62 北京理工大学博士学位论文会对蜜蜂造成伤害。3.7流动显示实验实验记录了圆柱绕流和蜜蜂振翅的流场图片。在圆柱绕流实验中使用了宽度2mm～3mm的薄片状电阻丝，在蜜蜂振翅实验中使用的是0.1mm粗细的圆柱状电阻丝。可以观察到两者发烟效果的差别。薄片状阻丝发出均匀的烟雾面，圆柱状阻丝形成条纹明显的一条条烟线。3.7.1圆柱绕流实验圆柱绕流时，由于流动分离和旋涡破裂而在圆柱下游两侧不断交替脱落的卡门涡街，从而产生两排方向相反交错排布的旋涡尾迹。雷诺数Re决定了圆柱绕流的流动特5征，60。这表明蜜蜂能产生大于重力的升力。蜜蜂被系住时，振翅频率与最大拍动角均可能低于自由飞行时的数据。当蜜蜂自由飞行时，通过改变振翅频率与最大拍动角，可以产生更大的升力，这也是蜜蜂能够携带一定重物飞行的原因。4.6本章小结在蜜蜂的流场显示实验中，利用两台高速摄像机拍摄了蜜蜂两个角度的图像，建立了蜜蜂坐标系，根据直线与投影的关系，求得蜜蜂振翅时拍动角φ的变化规律，在一个周期内，拍动角φ随时间的变化曲线近似为半个周期的正弦曲线。对蜜蜂进行风洞流场显示实验，利用高速摄像机拍摄记录了远场和近场的流场结构。在远场流场显示实验中，拍摄观察到的尾迹结构包括下拍形成的翅尖涡涡管、翅外旋时产生的下拍停止涡和上拍启动涡、翅扭转产生的环量方向相反的前翅尾缘涡和后翅尾缘涡。定性地分析了尾迹涡对翅的诱导速度，在有来流的蜜蜂振翅运动中，只需考虑一个周期的尾迹效应。在近场流场显示实验中，对比了不同展向位置的前缘涡尺寸，沿着翅的展向，由翅根至翅尖，前缘涡尺寸逐渐增大，不停地向翅尖输送涡量，在翅尖处与翅尖涡相连接。下拍末期，下拍停止涡形成于翅上表面，对升力贡献为负；上拍初期，上拍启动涡形成于翅下表面，能够提高升力。翅尖涡形成后，在一个周期内，翅尖涡尾迹的位移量很小。最后，根据流场显示实验结果，提出了蜜蜂飞行时涡的演变规律：翅外旋时，产生下拍停止涡和上拍启动涡；翅的扭转产生了较强的弦向变形，使得前翅尾缘与后翅尾缘处的涡量方向相反；下拍末期，在蜜蜂左翅和右翅上，均形成了由前缘涡、翅尖涡尾迹涡管和翅根涡尾迹组成的闭合涡环。根据涡环估算了翅产生的平均升力，平均升力大于蜜蜂的重力，这表明，蜜蜂通过振翅产生的升力，不仅能够支撑自身的重力，而且能够携带一定的重物飞行。92 北京理工大学博士学位论文第五章蜜蜂飞行时涡演变规律的仿真验证5.1引言研究昆虫扑翼飞行的非定常流场，直接测量翅上的力面临非常大的挑战，首先，昆虫的尺寸较小，这要求测量力的装置也要足够小，鉴于直接测量昆虫翅的力非常困难，随着计算机性能的发展，越来越多的学者利用CFD方法对扑翼飞行的问题进行[148][156]数值模拟求解。特别是对于昆虫扑翼飞行时周围涡流的影响。Lehmann、Sane、[205][206][207][208][209]Wang、Iima、Minotti、Keigo和Lee等人将扑翼运动简化为二维的运动，通过求解N-S方程，分析了昆虫飞行时脱落的尾迹涡特征和扑翼周期运动的气动[210][211][212][213]力特性。Miller、Ramamurti、孙茂和Liu等人对三维扑翼运动进行了数值模拟。Aono等人基于有限体积法求解N-S方程，数值模拟了昆虫悬停时近场和远场的三维涡流变化规律。他指出在下拍和上拍初期，有马蹄形的涡成长为环形涡环，在涡环中心有强烈的射流，形成下洗气流；最后环形涡环分解为两个圆形涡环。上拍行程中有较大的升力，可能是存在较强的前缘涡和翅尖涡，与下拍行程相比有较强的下洗[95]气流。根据计算的气动力矩、惯性力矩和功率，估算了悬飞的能量。Ramamurti等人数值计算了果蝇悬飞时的三维非稳态气动力，计算结果与实验结果基本一致。转弯时，左翅和右翅的拍动角与偏差角有较小的差值，这导致了偏航（yaw）力矩，在转弯初期偏航力矩最大。根据左翅和右翅的压心位置与力臂，研究了偏航力矩，结果表明，向前的力和升力的分量共同产生了转弯力矩。从果蝇翅前缘[211]和翅尖脱落的涡量形成了环状结构，当翅旋转到下个半拍时，环状结构被破坏。于鑫和孙茂数值计算模拟了昆虫飞行的流场，结果表明，昆虫在悬停和前飞时，在下拍和上拍期间，左翅和右翅均会产生一个涡环，涡环中心有很强的射流，在涡环[214]外的气流则相对较弱。左右翅产生的两个涡环在昆虫身体两侧，并不相连接。在第四章中，对蜜蜂进行了流场显示实验，根据蜜蜂振翅时烟线的变化规律，分析了蜜蜂扑翼飞行时周围的流场结构，特别是涡流的变化情况，提出了蜜蜂振翅时产生的涡的演变规律。但是，蜜蜂振翅的烟线实验属于定性分析实验，在流场显示实验中，烟线截面是二维截面，虽然利用了2台高速摄像机在两个角度同步拍摄，但对于93 北京理工大学博士学位论文实验所提出的三维涡结构，仍需要借助流体计算动力学（CFD）方法进行验证。为此，需要建立蜜蜂动力学模型和计算模型，通过求解计算N-S方程，既能分析蜜蜂周围的流场变化是否与所提出的涡演变规律一致，又可求解蜜蜂振翅产生的力。在烟线流场显示实验中，蜜蜂振翅时既有附着涡（前缘涡、尾缘涡和翅尖涡），又有脱落的涡形成的三维尾迹结构（翅尖涡形成的涡管），利用数值模拟计算能够更好地研究涡的形成与脱落，特别是更有利于研究涡的三维结构，能够分析前缘涡由翅根至翅尖的分布情况，翅上拍和下拍时翅尖涡的尾迹，以及翅的边界涡与尾迹涡的关系。昆虫飞行时，身体的姿态随着运动方式的不同发生改变，例如蜜蜂在悬停、前飞、后飞等不同状态包括不同机动情形下，身体的方位角以及拍动平面（近似看成平面）的方位角都会发生变化。当昆虫处于稳定的飞行状态时，一般认为其双翅做对称的拍动；当昆虫做偏航（yaw）以及滚转（roll）等机动动作时，认为双翅的运动会有些微差别，经过高速摄影观察也可以发现这一点，目前研究昆虫飞行的机动动作比较困难，几乎所有研究都是关注稳定的飞行状态，如悬停、前飞等。当昆虫处于稳定的飞行状态时，双翅的运动对称，身体各部分的方位角和拍动平面的方位角共同构成了昆虫的飞行姿态。在计算中仅考虑蜜蜂翅在拍动平面内的拍动运动。在数值模拟过程中，将蜜蜂翅的运动近似为刚体的定点转动，拍动角的变化为翅绕着翅根点的公转运动，翅还有绕主翅脉的自转运动。公转运动称为翅的平动，对应拍动角φ的变化；自转运动称为翅的翻转，对应攻角α的变化。值得指出的是，真实蜜蜂翅的自转运动伴随着翅的柔性变形，从翅根至翅尖方向，不同截面的变形均不同，近似的刚体翅只能将其简化为攻角的运动。因此，攻角α只是为了描述翅的翻转引入的近似物理量。根据蜜蜂振翅拍动角的实验数据（第四章已给出）和攻角近似理论，给出了蜜蜂翅的运动方程，通过优化方程能够调节振翅频率、最大拍动角和中期攻角等值，便于研究扑翼运动的流场规律。本章选用标准k-e湍流模型，利用动网格和网格重生技术求解N-S方程。为验证数值计算和选用参数的可靠性，将计算结果与孙茂的计算数据和Dickinson的实验数据进行了对比。计算结果与孙茂和Dickinson的结果有较高的一致性，这表明选用的湍流模型与设置参数能用于研究蜜蜂振翅飞行的流场结构和气动特性。本章数值计算模拟了蜜蜂振翅时附近流场涡流的变化情况，根据流线图分析了附着于翅上的尾缘涡、前缘涡和翅尖涡；分析了翅外旋时，尾缘处的下拍停止涡和上拍启动涡；根据等涡量面三维云图，分析了蜜蜂振翅过程中三维涡环结构。对比数值模拟与流场显示实94 北京理工大学博士学位论文验结果，一方面验证了流场显示实验中提出的蜜蜂振翅时涡的演变规律，另一方面补充了实验中未能观察到的现象，例如翅根涡和翅上拍形成的涡环结构。5.2蜜蜂模型的坐标系蜜蜂飞行是多自由度的运动，为了研究蜜蜂飞行时周围的流场，需要简化蜜蜂模型，建立蜜蜂的运动学模型。在数值模拟计算中，忽略蜜蜂触角和腿的影响，将蜜蜂简化为身体与翅膀两部分，翅在运动过程中忽略柔性变形，视作一均匀厚度的刚体平板，如图5.1-1所示，红色的部分是蜜蜂身体，黄色部分是翅。如图5.1-1所示，以左翅根为原点建立坐标系，xyz为蜜蜂翅的随体坐标系，xy面与翅平面重合，y沿翅脉方向，x沿翅的弦向。XYZ为惯性坐标系，XY为水平面。翅的运动轨迹近似为一平面，称为平均拍动平面，如图5.1-1中蓝色椭圆线所示。平均拍动平面与水平面（图中黑色椭圆线）的夹角称为拍动平面角β。翅的位置由拍动角φ和攻角α来确定。翅尖在拍动平面上扫过的角度称为拍动角，用φ表示，翅面与拍动平面的夹角为攻角，用α表示。（1）（2）图5.1蜜蜂模型的坐标系和拍动平面图5.1-2为拍动平面俯视图。拍动分为下拍和上拍，下拍时，拍动角φ从0变化到φ；上拍时，拍动角由φ变化到0。下拍时，拍动角φ从0变化到φ；上拍时，maxmaxmax拍动角φ由φ变化到0。φ是翅的最大拍动角。翅在上下拍动交替时，翅会翻转，maxmax翻转分为内旋和外旋，由下拍变为上拍的翻转过程是外旋，反之是内旋。惯性坐标系XYZ与随体坐标系xyz的转换公式如式5.1所示。95 北京理工大学博士学位论文Xxcosφφ−sincosαsinαYy=sinφφcos1Zz1−sinααcos（5.1）coscosφα−sinφcossinφαx=sincosφαcosφsinsinφαy−sinααcosz5.3翅的运动方程5.3.1简谐运动方程翅的位置由拍动角φ和攻角α确定，要确定翅的运动方程，就要构造φ和α关于时间t的函数，拍动周期记为T。根据第四章蜜蜂拍动角的变化规律，拍动角φ关于0时间t的函数为半周期的正弦曲线。简单起见，假设蜜蜂翅的运动为简谐振动，满足三角函数关系。根据蜜蜂振翅规律，t=0的起始时刻，y轴与Y轴重合，即φ=0，α=π/2；当tT=/4时，φφ=/2，αα=；当tT=/2时，拍动角达到最大φφ=，α=π/2；0maxmin0max当tT=3/4时，φφ=/2，αα=−π。0maxmin经过待定系数法计算可得到以下运动方程：φφmaxmax2π3φ=++sin(tπ)22T02（5.2）αα=+−π(π)sin(2πt+π)min22T0πφmax2π3φ=sin(t+π)TT002（5.3）αα=2π(π−+)sin(2πtπ)minTT200式5.2为拍动角φ和攻角α随时间t的方程，式5.3为拍动角速度φ和攻角速度α随时间t的方程。在式5.2和式5.3中，可变参量包括最大拍动角φ、拍动周期T和最max0小攻角α，给定这三个变量，可确定翅的运动方程，图5.2为简谐运动的角度和角min速度曲线图。图5.2-1为拍动角φ和攻角α的周期曲线图（式5.2的曲线图），图5.2-2为拍动角速度φ和攻角速度α的周期曲线图（式5.3的曲线图）。如图5.2所示，在1个振翅周期内，蜜蜂翅的拍动角φ、攻角α、拍动角速度φ和96 北京理工大学博士学位论文攻角速度α均随时间在变化。控制翅的运动只能调节φ、α和T三个变量，为了maxmin0更有利于控制翅模型的运动，对翅的运动方程进一步优化。（1）（2）图5.2简谐运动的角度和角速度曲线5.3.2中期拍动速度和攻角不变的方程现假设振翅时，拍动速度φ先从0加速、然后匀速、再减速到0从而完成半个周期；而攻角α则与拍动速度φ满足类似的变化趋势，攻角α亦经历了变化和保持不变的阶段。在1个周期内，令拍动速度φ对时间t的加速（包括减速）曲线取简单正弦函数，在剩余时间内拍动速度φ保持不变，在1个周期内，拍动的变速阶段组合起来刚好是正弦函数的1个周期，记为T（TT≥>0），拍动速度φ保持不变的时间为101TT−。拍动速度φ在保持不变时取其最大值φ，其中φ为最大拍动角。01maxmax在式5.2中，攻角α一直随时间t变化。在蜜蜂振翅实验中观察到，翅在上下拍动交替时，翅的翻转运动最明显，在下拍和上拍中期，翅的翻转运动不明显。因此可以假设攻角α在下拍和上拍中期保持不变，这与孙茂和Dickinson的模型假设相同。攻角α的变化时间小于拍动角φ的变化周期。1个拍动周期内除去翻转（攻角α改变）时间外，在剩余时间内攻角α保持不变，这称为翅的平动。设翻转周期（攻角α的变化时间）为T，下拍中期攻角αα=，上拍中期攻角αα=−π。构造攻角速度α对2minmin时间t的函数，使其处处光滑。令攻角速度α对时间t的加速曲线为周期等于T的余弦2函数（TT≥>0）；攻角速度α的最大值为α；攻角在保持不变时取其最小值α。02maxmin令t+T1TT411tt=−tTint,∈−[，T−)（5.4）1010T44097 北京理工大学博士学位论文t+T2TT422tt=−tTint,∈−[，T−)（5.5）2020T440式5.4和式5.5中，int()x为取整函数，表示不超过x的最大整数。拍动速度φ和攻角速度α如式5.6所示。2πTTφsint,t∈−[11，)max11T441TTTφ,[t∈−11，0)max1424φ=2πT0TT00TT11φsin−tt,∈−[，+)max11T122424TTT−φ,[t∈+011，T−)max10（5.6）244ααmaxmax4πTT22−−costt22,∈−[，)22T442TTT2200,[t∈−，)2424α=ααmaxmax4πT0TT0TT220+costt−,∈−[，+)2222T222424TTT0220,[t∈+，T−)20244当t=0时，拍动角φ=0；当t=T/2时，φφ=。对拍动速度φ积分，可得待110max定参数φ，如式5.7所示。max2πφφ=max（5.7）max2TTT+−π()101当t=0时，攻角α=π/2；当t=T/4时，攻角αα=，得到待定参数α，222minmax如式5.8所示。4π−8αminα=（5.8）maxT2拍动角和攻角的表达式如式5.9所示。式5.9为式5.2优化后的运动方程。在式5.9中，有5个变量，分别为T、T、T、α和φ。通过改变这5个变量，可以调012minmax整翅的运动。在式5.9中，当TTT==时，式5.9等同于式5.2。12098 北京理工大学博士学位论文Tφ2πTT1max111cos−tt11,[,)∈−2TTT1+−π(01)T1442πφmaxTT11T1T0T1tt11+−,[,∈−)2TTT1+−π()201π4424φ=TTT+−π()Tφ2πTTTTT101φ+1maxcostt−0,∈−[011,0+)max112TTT1+−π()201TTT1+−π(01)T1224242πφTTTTTmax11011TT01+−−tt,[10∈+,−)（5.9）2TTT+−π()2π4244101πtπ−2α4πTT2min22−(2π−4αmin)−sintt2,2∈−[,)2TT222π44TTT220αt,[,∈−)min2424α=πt−−T0π2α4πTTTTT22min00022+(2π−+4α)sintt−,∈−[,+)min222T22πT222424TTT022π−αt,[∈+,T−)min20244图5.3为中期拍动速度φ和攻角α不变的角度和角速度曲线图。图5.3-1为拍动角φ和攻角α的周期曲线图（式5.9曲线图），图5.3-2为拍动角速度φ和攻角速度α的周期曲线图（式5.6曲线图）。如图5.3所示，在半拍中期，攻角α保持不变，拍动角速度φ保持不变，可以通过T和T调整φ和α保持稳定的时间。T衡量拍动速度φ变121化快慢；T衡量翻转速度α变化快慢。2（1）（2）图5.3中期φ和α不变的角度和角速度曲线如图5.3所示，tT=0~/2时，是翅的下拍阶段；tT=/2~T时，是翅的上拍阶00099 北京理工大学博士学位论文段。翅的平动阶段，拍动角速度φ达到最大值φ并保持恒定，攻角速度α为0。在翅max的翻转阶段，拍动角速度φ和攻角速度α分别按正弦和余弦曲线变化。外旋时（下拍变为上拍的翻转过程），攻角α由α变到πα−；内旋时（上拍变为下拍的翻转过minmin程），攻角α由πα−变到α。下拍中期，αα=；上拍中期，απα=−；其他minminminmin阶段攻角α按正弦函数的变化。5.3.3提前、对称和延迟翻转的运动方程蜜蜂主动飞行时，能通过改变翅翻转的时机从而改变气动力的大小来控制飞行。式5.2和式5.9中，翅的翻转在拍动前后是对称的，所以被称为对称翻转；在此基础上，翻转相对于拍动相位提前，就是提前翻转；翻转相对于拍动相位落后，就是延迟翻转。如图5.3-2所示，将翻转速度曲线向左平移即是提前翻转，反之则是延迟翻转。平移量由T决定，将翻转过程全部平移到原点左边，即曲线向左平移T/4，达到提22前翻转的临界，否则会有负的平动升力产生；反之，将翻转过程全部平移到原点右边，即曲线向右平移T/4，达到延迟翻转的临界。修改式5.6中的攻角速度α，可得式5.10。2αα4πTTmaxmax22−−costt33∈−[，)22T442TTT2200t∈−[)，3424α=αα4πTTTTTmaxmax00220+costt33−∈−[，+)（5.10）22T222424TTT0220t∈+[)，T−30244Tt+−(1ϑ)T224t=−tϑ−Tintϑ∈−[11]，304T0对时间t进行平移，定义平移因子ϑ，在-1到1之间变化。ϑ<0时为提前翻转；ϑ>0为延迟翻转；ϑ=0时为对称翻转；ϑ=1为翻转时机的临界点；ϑ>1时，翅产生负的平动升力，对于描述昆虫飞行来说脱离现实意义，但是为了特殊目的的人工模型可以做到ϑ>1，暂不考虑这种情形。100 北京理工大学博士学位论文图5.4提前、对称、延迟翻转的速度变化示意图Tφ2πTT1max111cos−tt11,[,)∈−2TTT1+−π(01)T1442πφmaxtt+−TT11,[,∈T1T0−T1)112TTT+−π()2π4424101φ=TTT10+−π(1)T1maxφ2πT0TT0TT110φ+costt−,∈−[,+)max112TTT+−π()2TTT+−π()T2242410110112πφTTTTTmaxTT+−−11tt,[∈+01,−1)（5.11）01102TTT+−π()2π4244101πtπ−2α4πTT−3(2π−4α)−minsintt,∈−[2，2)min332TT2π4422TTT220αtmin,[3∈−，)424α=πt−T0π−2α4πTTTTT32min00022+(2π−+4αmin)sintt33−,∈−[，+)22T2πT222424TTT022π−αt,[∈+，T−)min30244t+T14t=−tTint10T0Tt+−(1ϑ)T224t=−tϑ−Tint,ϑ∈−[11]，304T0对于速度曲线需要给出初始条件，即t=0时翅的位置。前文中已经提到拍动角和对称翻转时初始条件（当ϑ=0，t=0时，α=π/2）。对于图5.4中给出的提前翻转，翅在初始时刻为最小攻角，即当ϑ=−1，t=0时，αα=；图5.4中的延迟翻转（ϑ=1），min当t=0时，αα=−π。对于一般情形，将式5.9中攻角α的时间做与式5.6相同的min101 北京理工大学博士学位论文处理，得到通过平移因子ϑ改变相位的攻角α与时间t的关系，如式5.11所示。图5.5为不同翻转时刻下的角度变化示意图。式5.11中代入t=0，求得初始条件，如式5.12所示。πϑπ−2αminααϑt=0=+−+(2π4min)sin(π)（5.12）242π图5.5提前、对称、延迟翻转的角度变化示意图5.4边界条件和求解器蜜蜂的身体和翅均采用无滑移壁面边界条件（wall），并通过UDF函数给定翅的运动，计算域边界条件为自由出流（outflow），假设蜜蜂左右翅的振翅运动是对称的，在蜜蜂身体中心对称面采用对称边界条件（symmetry）。边界条件示意图如图5.6所示。图5.6边界条件示意图由于扑翼模型运动复杂，需要捕捉三维流场特征，因此要对整个流场划分三维网格。扑翼模型拍动频率高，快速翻转时，翅上的网格要随之运动，网格需要不断刷新和重构，选择具有良好适应性的四面体非结构网格划分整个计算域。非结构网格是指[215]网格单元和节点彼此没有固定的规律，节点公布完全是任意的。任何三维空间区102 北京理工大学博士学位论文域都可以被四面体单元的网格所划分，能够方便地生成复杂外形的网格，能够通过流场中的大梯度区域自适应来提高对间断的分辨率，分区及并行计算比起结构网格更直接。但是在相同的体积下，划分四面体非结构网格比六面体结构网格需要的网格数量[216]要多，而且在相同网格数量前提下，计算非结构网格所需内存更大，周期更长。非结构网格由于舍去了网格节点的结构性限制，易于控制网格单元的大小、形状以及网格点的位置，因此比结构网格具有更大的灵活性，对复杂外形的适应性更强。另外对于结构网格，在计算域内网格线和平面都应保持连续，并正交于物体边界和相邻的网格线和面，而非结构网格无此限制，这样就消除了网格生成中的一个主要障碍，且其网格中一个点周围的点数和单元数都不是固定的，可以方便地做自适应计算，合理分配网格的疏密，提高计算精度。利用Gambit划分网格，在Gambit中生成网格的过程是先划分翅上特征线上的网格，由线网格再生成面上的网格，然后根据面网格再生成体网格。对于动网格计算的特殊性，要求网格全部为非结构网格，即面网格全部为三角形网格，体网格全部为四面体网格。如图5.7所示，图5.7-1为蜜蜂身体的面网格，图5.7-1为翅上的网格。（1）（2）图5.7网格示意图模拟扑翼模型运动的流场是不可压的，因此，选择压力求解器，采用隐式求解方法。扑翼模型运动是周期性的，是非定常过程，流场在每一瞬时都会发生变化，选择非稳态模型。在计算过程中采用了动网格技术，通过编写UDF函数定义扑翼模型的运动轨迹，在每个时刻调用UDF赋值每个时刻扑翼模型上网格节点的位移。利用弹[217]簧系数法和网格重生技术更新每个时刻的网格。弹簧系数（SpringConstantFactor）选为0.8。动网格技术主要用于模拟边界在不同时刻下的运动情况。基于新的边界位置，在[218]每一时间步自动更新网格位置，从而对流场进行计算。在随体坐标系生成网格计103 北京理工大学博士学位论文算的过程中，随着时间的推移，翅模型需按照指定的轨迹进行运动，因此需要在计算中识别翅膀的运动方程。计算中采用了DMM（Dynamic-MeshModal）功能函数，通过DMM模块读入翅在每一时刻的位置参数。翅模型采用的是刚体运动，采用CG_MOTION函数来控制刚体运动。翅膀在运动过程中，翅膀周围的网格时时刻刻在发生变化，为了使在任意时刻的网格能够满足解算流场的条件及迭代过程的收敛性，采用弹簧系数法和网格重生技术更新每个步长的网格。弹簧系数算法是指将网格节点之间的边等效为一根弹簧，通过调整计算域内部的网格节点通过边界网格节点的刷新，弹性系数法在更新体网格的过程中不会改变网格的连续性。当边界结点的位移相对局部网格的尺寸很大时，网格的质量将变得很差。为避免这一问题，采用局部重构算法对质量差的网格进行合并或拆分。不合要求的网格定义为超过某一给定体，低于某一指定体积或者网格倾斜率大于某一数值的那些网格。流场计算基于有限体积法求解N-S方程，有限体积法的求解策略就是用边界面或线上的通量计算出控制点上的变量。将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积[219]分，必须假定值在网格点之间的变化规律，即假设值的分段的分布剖面。蜜蜂的飞行速度较低，整个流场中，空气密度变化很小，可以近似的看做常数，因此，把流场中的空气当做不可压缩气体处理。在扑翼模型的整个运动过程中，流场是低雷诺数下发展的完全湍流，通过求解雷诺平均N-S方程计算扑翼模型在完整周期内的气流运动。雷诺平均N-S方程是流场平均变量的控制方程，其相关的模拟理论被称为湍流模式理论。湍流模式理论假定湍流中的流场变量由一个时均量和一个脉动[220-221]量组成，以此观点可以得到雷诺平均N-S方程。N-S方程具体形式如下：∂uj=0∂xj（5.13）∂uupii∂∂∂()tr∂+ruj∂∂∂=−+()ttij+ijtxxxjij104 北京理工大学博士学位论文()t2在式5.13中，分子应力tij=2µlSij，雷诺应力tij=2µlSkij−rδij，31∂∂uuijSij=()+2∂∂xxji。采用标准k−ε模型做湍流模型。通过求解湍流动能k方程和湍流耗散率ε方程，得到k和ε的解。∂∂∂µt∂k(rk)+(rku)=[(µ+)]++−−+GGrεYSikbMk∂∂tx∂xσ∂xiikj∂∂∂µ∂εεε2t(rε)+(rεu)=++[(µ)]C(GCG+−+)CrSi1εkb32εεε∂∂tx∂xσ∂xkkiijε（5.14）2kµ=rCtµεω∂∂uu""iiε=()()r∂∂xxkk∂µ""jG表示由层流速度梯度而产生的湍流动能，G=−rµµ，G是由浮力产生kkijb∂xiµt∂T2的湍流动能，Gg=β，Y是由过渡的扩散产生的波动，YM=2rε，根据biMMtPr∂xti经验值，对模型参数取值C=1.14、C=1.92、C=0.09、σ=1.0、σ=1.3。1ε2εµkε再根据Boussinesq假设，如式5.15，()t∂∂uuii∂uj2tµ=(+−+)(rkµδ)（5.15）ijttij∂∂uu3∂xjii联合式5.13、式5.14、式5.15，即可求得雷诺应力的解。5.5计算有效性验证[30][84]为验证计算的可靠性，与孙茂的计算数据和Dickinson的实验数据进行对比。验证模型选用果蝇翅模型（与孙茂的计算模型一致，如图5.8所示）。翅形状参数：翅2长R=3mm，翅厚δ=0.1mm，单翅面积S=3mm。翅的运动参数：最大拍动角φ=150°，拍动中期攻角α=40°，拍动频率fT=1/=240H，Tf=0.36/，maxmin0Z13Tf=0.72/。参考速度U=2φfr=2.187m/s，rR=0.58，空气密度r=1.223kg/m，222105 北京理工大学博士学位论文2升力系数C=F/0.5rUS。LL图5.8果蝇翅模型的平面形状根据孙茂的数据，分别计算果蝇翅模型在提前翻转（ϑ=−0.5），对称翻转（ϑ=0）和延迟翻转（ϑ=0.5）时的升力系数C，与孙茂和Dickinson的结果进行对比。图L5.9为计算升力系数对比图。图5.9-1、图5.9-3和图5.9-5分别为本文计算的提前翻转、对称翻转和延迟翻转时的升力系数曲线图，图5.9-2、图5.9-4和图5.9-6分别为孙茂和Dickinson提前翻转、对称翻转和延迟翻转时的升力系数曲线图。在图5.9-2、图5.9-4和图5.9-6中，实线表示孙茂的计算升力系数曲线，虚线表示Dickinson的实验结果曲线。如图5.9所示，本文的计算结果与孙茂和Dickinson的数据基本吻合，这表明计算所选用的湍流模型和设置参数可用于研究扑翼模型的计算。（1）（2）（3）（4）106 北京理工大学博士学位论文（5）（6）图5.9果蝇扑翼模型升力系数对比图5.6涡现象分析5.6.1蜜蜂振翅产生的涡根据蜜蜂振翅的数据，翅的运动参数：最大拍动角φ=120°，振翅频率maxf=190Hz，TT=。目测蜜蜂中期攻角大约为30°，取下拍中期攻角α=30°，10minTT=/2，ϑ=0。参考时间t=cU/，参考速度U=2φfr，参考长度cSR=/，20max22rR=0.61，根据1.5.1节翅的参数，翅长R=9.81mm，翅的面积S=26.89mm。翅的2运动方程如式5.11所示，速度曲线如图5.10所示，1个周期中，翅的运动包括下拍（Downstroke）和上拍（Upstroke）；每个半拍中期，攻角速度为0，攻角保持不变，这个阶段为平动（Translation）；攻角速度改变时，为翅的翻转阶段，下拍转变为上拍的过程，称为翅的外旋（Supination），上拍转变为下拍的过程，称为翅的内旋（Pronation）。由于昆虫翅的面积小，速度快，在实验中捕捉翅周围的流场运动困难且复杂。根据5.4节的边界条件和设置参数数值计算模拟研究蜜蜂振翅时周围的流场。在计算过程中，蜜蜂身体边界网格没有位移，组成翅面的网格随翅的运动随时更新。图5.10蜜蜂翅的运动速度曲线图107 北京理工大学博士学位论文在1个拍动周期内，根据蜜蜂翅的流线图，根据涡的位置分类，翅上附着有前缘涡（LEV）、翅尖涡（TV）、翅根涡（RV）和尾缘涡（TEV），如图5.11所示。在流场显示实验中，当烟线截面靠近翅根处时，蜜蜂身体对烟线干扰作用较强，因此实验中高速摄像机拍摄的烟线图主要观察分析LEV、TV和TEV，需要借助CFD方法补充研究RV。（1）翅尖涡（2）前缘涡（3）翅根涡（4）尾缘涡图5.11涡的流线图计算蜜蜂振翅的流场，主要分析涡的演变过程。翅的流线图最能直观显示翅上的附着涡。在计算周期振翅流场时，需要保存每个时间步长的流场结果文件，在1个周期中选取10个时刻，通过对每个时刻的流场文件后处理，显示翅的流线图，分析涡的演变过程，如图5.12所示。图5.12是蜜蜂在1个拍动周期的流线图，不同颜色代表不同的流线，TT=。翅0在下拍初期内旋阶段（tT=0.1）和上拍初期外旋阶段（tT=0.6），有TEV附于翅面上，且在翅的最宽截面处。翅内旋时，TEV附着于翅的上表面，翅外旋时，TEV附着于翅的下表面。在翅下拍末期（tT=0.5）和上拍末期（tT=），由翅根沿尾缘向下有RV。在1个周期中，在每个时刻，均存在翅尖涡（TV）。当LEV存在时，沿着翅的展向，TV与LEV连接在一起。在下拍阶段，TV附着于翅的上表面，如图5.12-1～图5.12-5所示，在上拍阶段，TV附着于翅的下表面，如图5.12-6～图5.12-10所示。108 北京理工大学博士学位论文在下拍初期（tT=0.1），LEV不明显，随着翅向下拍动，LEV尺寸增长（tT=0.2），附着于翅上表面前缘处，在下拍中期（tT=0.3），由翅根至翅尖方向，LEV尺寸逐渐增大。翅开始旋转（tT=0.4）时，LEV的尺寸逐渐增大，在下拍末期（tT=0.5），LEV尺寸达到最大，如图5.12-5所示。翅旋转由下拍转变为上拍后（tT=0.6），LEV逐渐脱落，如图5.12-6所示。在翅上拍中期，LEV附着于翅下表面前缘处，如图5.12-7～图5.12-9所示。翅开始内旋时，LEV尺寸增长，在上拍末期（tT=），LEV尺寸在整个上拍阶段最大，如图5.12-10所示。根据1个拍动周期中的流线图，LEV在下拍和上拍中期分别稳定地附着于翅上表面和翅下表面，翅的旋转使得LEV在上拍末期和下拍末期的尺寸最大，到下个半拍（下拍变为上拍或上拍变为下拍）后，LEV脱落，在平动阶段，LEV又附着于前缘处。（1）t=0.1T（2）t=0.2T（3）t=0.3T（4）t=0.4T（5）t=0.5T（6）t=0.6T109 北京理工大学博士学位论文（7）t=0.7T（8）t=0.8T（9）t=0.9T（10）t=T图5.12一个周期的流线图5.6.2前缘涡图5.13-1～图5.13-3是蜜蜂翅在下拍中期（tT=0.3）的不同展向截面等涡量图。黑色圆点表示翅的前缘，红色表示高涡量值，蓝色表示低涡量值，无量纲涡量值（Non-dimensionVorticity）范围为Vor=−5~5。图5.13-1、图5.13-2和图5.13-3分别为沿着翅的展向截面为0.25R、0.5R和0.75R的涡量云图。图5.13-4、图5.13-5和图5.13-6分别烟线经过翅的展向截面为0.25R、0.5R和0.75R的流场显示实验图。如图5.13-1～图5.13-3所示，翅的上表面处，涡量值为正；翅的下表面处，涡量值为负。在翅前缘处，涡量明显高于其他区域，且涡量为正，这就是前缘涡存在的区域。在图5.13的二维截面视图中，LEV的环量为顺时针。截面越靠近翅尖，LEV越明显，LEV尺寸也越大。在图5.13-2和图5.13-3中，在翅前缘处，前缘涡绕过前缘附着于翅的上表面处，前缘涡尺寸均明显大于图5.13-1。在流场显示实验中，从翅根至翅尖，LEV尺寸也在逐渐增长，如图5.13-4～图5.13-6所示。数值模拟计算结果与实验结果是一致的。110 北京理工大学博士学位论文（1）0.25R（2）0.5R（3）0.75R（4）0.25R（5）0.5R（6）0.75R图5.13不同展向位置的前缘涡（1）（2）（3）（4）图5.14前缘涡的三维结构图5.14是翅下拍末期（tT=0.4）的数值模拟的前缘涡流线图和示意图。图5.14-1和图5.14-2是翅两个视角的流线图，在沿着翅前缘的展向方向，翅根至翅尖LEV的尺寸逐渐增大，在翅尖处LEV与TV相连接，黑线标识出了LEV的三维结构。图5.14-3是111 北京理工大学博士学位论文蜜蜂双翅上的流线示意图，翅前缘处存在环量，形成了翅上的附着涡主要为LEV和TV，提供了飞行所需的升力。在蜜蜂流场显示实验中，高速摄像机拍摄了在不同展向位置的LEV结构，实验结果表明，越靠近翅尖方向，LEV尺寸越大。数值模拟的结果与流场显示实验结果完全一致，如图5.14-4所示，蜜蜂飞行时，在下拍末期，翅前缘处，由翅根至翅尖均有LEV，翅根处LEV尺寸最小，越靠近翅尖，LEV尺寸越大，在翅尖处LEV与TV连接在一起。5.6.3下拍停止涡与上拍启动涡图5.15-1和图5.15-2分别为下拍末期（tT=0.4）和上拍初期（tT=0.6）时翅0.5R截面处的涡量云图，图中黑色圆点表示翅的前缘。在下拍末期，如图5.15-1所示，翅尾缘处环量方向为逆时针，即下拍停止涡（DSV）；在上拍初期，如图5.15-2所示，翅尾缘处环量方向为顺时针，即上拍启动涡（USV）。图5.15-3和图5.15-4分别是DSV和USV的示意图。图5.15-5和图5.15-6分别是DSV和USV的实验图片。如图5.15所示，DSV和USV的数值模拟结果与风洞流场显示实验结果一致。即在翅由下拍转为上拍的外旋阶段，下拍末期，尾缘处有DSV，速度环量为逆时针，上拍初期，尾缘处附着有USV，速度环量为顺时针。根据儒可夫斯基升力定理，DSV降低升力，USV提高升力。（1）（2）（3）（4）112 北京理工大学博士学位论文（5）（6）图5.15DSV和USV的数值模拟与实验对比图图5.16-1和图5.16-2分别为下拍末期（tT=0.4）和上拍初期（tT=0.6）时翅的流场图。如图5.16-1所示，流线图中并未见到尾缘涡，这表明此时翅尾缘处虽有环量，但是环量值较小。在图5.16-2中，根据流线所示，尾缘涡附着于翅最宽截面的尾缘处，这表明，在上拍初期，USV的环量较高。因此，USV对升力的提升大于DSV对升力的降低。（1）（2）图5.16尾缘涡的流线图5.6.4涡环图5.17为蜜蜂在一个周期内的等涡量面图。图5.17-1～图5.17-4分别为翅下拍末期、上拍初期、上拍末期、下拍初期的等涡量面云图，图中灰色表示蜜蜂身体，蓝色表示无量纲涡量Vor=1的等涡量面。如图5.17-1所示，在下拍末期，翅尖涡尾迹涡量、翅根涡尾迹涡量和翅的边界涡量连接在一起形成了涡环。翅尖轨迹决定了涡环的边缘及位置。翅运动到上拍初期，如图5.17-2所示，由于能量的耗散，翅尖涡尾迹与翅根涡尾迹涡量降低，涡环中心部分变空，形成了圆形的涡环。涡环处的气流均向涡环中心聚拢，如黑色箭头所示。如图5.17-3所示，在翅上拍末期，上拍过程中形成的翅尖涡尾迹涡量、翅根涡尾迹涡量和翅的边界涡量连接在一起，形成了上拍的尾迹涡环，涡环的形状像耳朵，称113 北京理工大学博士学位论文为耳状涡环。翅运动至下拍初期，如图5.17-4所示，旋转运动使得翅边界处的涡量增加，空气的粘性使得尾迹涡量耗散，耳状涡环中心部分变空。（1）t=0.4T（2）t=0.6T（3）t=0.9T（4）t=1.1T图5.17等涡量面图为分析流场的速度特征，取上拍初期的截面1和下拍初期截面2分析，如图5.18-1和图5.18-2所示。图5.18-3和图5.18-4分别为截面1和截面2处的速度矢量云图。如图5.18-3和图5.18-4所示，在涡环中心区域，有向下的速度射流，射流左右两边速度环量方向相反。（1）t=0.6T（2）t=1.1T114 北京理工大学博士学位论文（3）截面1（4）截面2图5.18速度矢量云图（1）Vor=1（2）Vor=2（3）Vor=3（4）Vor=4（5）Vor=5图5.19t=0.4T时的等涡量面云图（1）Vor=1（2）Vor=2（3）Vor=3（4）Vor=4（5）Vor=5图5.20t=0.9T时的等涡量面云图在流场显示实验中，根据烟线图分析，在下拍末期，前缘涡与翅尖涡尾迹连接在一起，如图4-19所示。数值模拟中，在下拍末期和上拍末期的等涡量云图中，如图5.17所示，翅尖涡尾迹和翅的附着涡量连接在一起，LEV并不明显。取tT=0.4和tT=0.9时刻分析，显示不同的Vor面。图5.19-1～图5.19-5为tT=0.4时Vor=1、2、3、4、5的等涡量面云图。图5.20-1～图5.19-5分别为tT=0.9时Vor=1、2、3、4、5的等涡量面云图。如图5.19和图5.20所示，随着Vor增大，等涡量面值的区域就越少，翅面的附着涡越明显。当Vor=5时，LEV结构较为明显，如图5.19-5和图5.20-5所示。翅附近的涡量值最高，距离蜜蜂身体和翅越远，尾迹的涡量值越小。与涡环（下拍圆形涡环115 北京理工大学博士学位论文和上拍耳状涡环）的涡量值相比，前缘涡的涡量高很多。在第四章的流场显示实验分析中，翅在下拍末期有前缘涡、翅尖涡尾迹涡管和翅根涡尾迹形成的涡环。数值模拟的结果与流场显示实验的分析结果是较为一致的。在数值模拟中，在上拍和下拍末期，翅尖涡尾迹、翅根涡尾迹和翅的附着涡量能够组成完整的涡环。下拍末期的涡环为圆形涡环，上拍末期的涡环为耳状涡环。在涡环中心，有垂直于涡环平面的较强射流速度。在数值模拟过程中，对翅上的压力积分，输出每个步长的力数值，计算得到在下−4拍过程中，单个翅上的平均升力为Ldown=7.09210N×，这与第四章估算的下拍升力−4−4Ldown=5.21810N×大小相近。蜜蜂重G=8.0710N×，2LGdown>，这表明蜜蜂被系住时，翅产生的升力不仅能抵消重力，还能产生额外的力以携带重物。5.7本章小结本章基于真实蜜蜂的飞行特点，建立了蜜蜂模型坐标系，给出了翅的运动方程。为了能研究振翅在不同运动状态下的气动特性，引入了多个变量，优化了运动方程，采用标准k−ε湍流模型，利用动网格技术求解了N-S方程，计算得到了蜜蜂振翅的流场。通过与孙茂和Dickinson的数据对比，验证了数值计算的可靠性。在一个拍动周期内，翅上附着有翅尖涡、翅根涡、前缘涡和尾缘涡。在一周期的任何时刻均存在翅尖涡。在上、下拍动初期，尾缘涡附着于翅的最宽截面尾缘处；在上、下拍动末期，有翅根涡。在下拍和上拍中期，前缘涡分别附着于翅上表面和翅下表面，翅的旋转使得前缘涡在下拍末期和上拍末期尺寸达到最大。数值仿真模拟结果表明：（1）下拍末期，前缘处有前缘涡，翅尖处有翅尖涡，由翅根至翅尖，前缘涡的尺寸增大，与翅尖涡连接在一起；（2）下拍末期，尾缘处有下拍停止涡，降低升力；上拍初期，尾缘处附着有上拍启动涡，提高升力；（3）下拍和上拍过程中翅尖涡尾迹涡管、翅根涡尾迹涡量和翅边界的涡量均能够组成完整的涡环，下拍末期的涡环为圆形涡环，上拍末期的涡环为耳状涡环，在涡环中心，有垂直涡环平面向下的射流。数值仿真模拟结果与流场显示实验结果的大部分结果是一致的：由翅根至翅尖，前缘涡尺寸增大；下拍过程中翅尖涡尾迹形成涡管；翅外旋时，尾缘处有下拍停止涡和上拍启动涡；翅下拍末期，翅的边界涡量与尾迹涡量形成了涡环。在流场显示实验中，未能观察到翅上拍过程中的前缘涡和翅尖涡，数值模拟结果对其进行了补充。116 北京理工大学博士学位论文综上所述，一方面数值模拟能验证流场显示实验中提出的三维涡结构和涡的演变规律是正确的，另一方面也表明所建立的蜜蜂数值仿真模型能够用于研究扑翼飞行机理和气动性能。117 北京理工大学博士学位论文第六章蜜蜂翅模型的力和功率分析6.1引言针对流场显示实验中所提出的三维涡结构与涡的演变规律，第五章已经验证并进行了补充。烟线实验研究蜜蜂飞行的流场属于定性研究，流场显示实验中能观察分析涡流的变化，特别是翅上附着涡的形成与脱落，这对提示扑翼飞行机理有重要作用。但是，要定量研究扑翼飞行产生的气动力和功率，还需要借助其他方法，例如PIV技术、机械扑翼的实验测量和数值模拟方法。近年来，PIV技术趋于发展成熟，在研究扑翼飞行时，已经有Spedding等人利用PIV技术重构随时间变化的三维流场结构，通过积分速度场可以求得扑翼飞行时的气动力。Dickinson团队对果蝇翅模型进行了油桶实验，测量了果蝇翅模型在不同运动参数下的气动性能。国内西北工业大学对自制的扑翼机进行了多种测量实验以研究扑翼飞行的气动性能。但是PIV测量实验和机械扑翼测量实验，均需要完整成熟的实验平台，对实验条件要求较高，成本也较高。随着计算机性能的发展，更多的学者利用数值模拟计算求[222]解N-S方程的方法研究了扑翼飞行的流场结构和气动性能。例如，国外Ramamurti、[223][94][95][224]Persson、Liu、Anno和Yu等人计算模拟了昆虫飞行时三维流场。国内孙[24-25]茂团队对昆虫悬停、前飞的产生的高升力机理进行了研究，南京航空航天大学昂[44-45]海松团队研究了柔性翼模型的气动特性和功率。本章利用CFD方法计算研究了蜜蜂翅模型在拍动时产生的气动力和功率。在计算过程中，假设蜜蜂的左翅和右翅的运动是完全对称的。因此，蜜蜂飞行时的力有两个方向的力，升力和推力。翅的功率分为气动功率和惯性功率两部分。蜜蜂飞行时，左翅和右翅的相互作用是否对翅的气动力和功率产生影响，以及蜜蜂身体是否影响翅的气动性能，均是值得研究的问题。本章分别建立了单翅模型、双翅模型以及完整蜜蜂模型，对比研究了不同模型条件下翅的气动力和功率，结果表明，左、右翅的相互作用和蜜蜂身体对翅的气动性能影响可忽略不计。在蜜蜂系飞的多组实验中能够观察到，蜜蜂振翅时，翅的运动参数会产生变化，例如最大拍动角，振翅频率等参数均会产生变化，因此有必要这些研究不同影响因素118 北京理工大学博士学位论文对翅气动性能和功率的影响。在第五章给出的翅运动方程中，可通过调节频率f、最大拍动角φ、攻角α、T、T、翻转时刻和旋转轴位置，调整翅的运动状态，以研max12究翅的气动性能和能耗，为未来研制微型扑翼飞行器提供指导。6.2气动力和功率分析6.2.1翅的气动力图6.1为翅坐标系。如图6.1-1所示，以翅根为原点，建立惯性坐标系XYZ，与翅固连在一起的随体坐标系xyz。XY面为翅的拍动平面，翅脉在XY面上扫过的角度为拍动角φ。主翅脉沿着y轴方向，翅平面与xy面重合，翅面xy与拍动平面XY的夹角为攻角α。图6.1-2是翅的拍动平面示意图。当蜜蜂左右翅上拍末期合拢时，沿着左翅脉为Y轴，如图6.1-2所示。如图6.1-1所示，在惯性坐标系XYZ中，作用于翅上的力可以分解为分别沿XYZ、、三个方向的力。在数值模拟研究中，假设左右翅的运动是对称的，左右翅在X方向的分力相互抵消，因此，可以忽略X方向的力。沿Z轴方向的力，即垂直于拍动平面方向的力，定义为翅的升力F，且规定沿Z轴正向的力为正。沿Y轴方向的力，L为翅的推力F，规定沿Y轴正向的力为正。T（1）（2）图6.1翅的坐标系在蜜蜂飞行过程中，可以调整拍动平面的位置，一般情况下拍动平面与水平面有一定的夹角，作用于翅上的气动力可以分解为水平力和竖直力。如图6.2所示，翅的拍动平面和水平面的夹角为β，翅的升力F和推力F沿着水平方向和竖直方向分解LT为水平力F和竖直力F。HV119 北京理工大学博士学位论文图6.2力分解示意图如图6.2所示，翅的拍动平面和水平面的夹角为β，将升力F和推力F沿着水平LT方向和竖直方向分解为水平力F和竖直力F。HVFF=sinββ−Fcos（6.1）HLTFF=cosββ+Fsin（6.2）VLT当蜜蜂悬停飞行时，平均竖直力F与蜜蜂重力大小相等，方向相反。平均水平V力F=0。因此平均升力F和平均推力F的合力为平均竖直力F，则有HLTVFTβ=（6.3）FL在本章研究内容中，只研究翅模型在不同影响因素下的气动性能，并不涉及蜜蜂的飞行状态，因此，在分析气动力时只研究升力F和推力F。LT6.2.2翅的功率翅的功率P（TotalPower）分为两部分，气动力功率P（AerodynamicPower）和ta惯性力功率P（InertialPower），即iPPP=+（6.4）tai翅的气动力功率为克服气动力矩所需的功率，等于气动力矩和角速度的矢量点乘，如式6.5。因此，一个拍动周期中的每个时间步长，需要计算对XYZ、、三个轴的气动力矩和XYZ、、三个轴的角速度。ddP==++MMMMωωωω（6.5）aaXXYYZZdd翅的惯性力功率P等于惯性力矩M与其对应角速度ω的点乘，即iiddPM=ω（6.6）iid对于翅的惯性力矩M，有：i120 北京理工大学博士学位论文dddHddMH=+×ω（6.7）idtxyzd对于翅的角速度ω，有：ddddωφααφα=−sinij++cosk（6.8）d对于翅的角动量H，有：dddHI=•ω（6.9）如图6.3所示，翅平面在xy面上，dm为翅的单位质量，翅的厚度很薄，翅模型[24]的厚度为0.05mm，因此在计算翅的质量矩时，可近似假设z=0。图6.3惯性矩示意图ddI是翅的质量矩，对于I的分量有=+≈222I()yzdmydmxx∫∫I=+≈()x22zdmxdm2yy∫∫I=()x22+ydmzz∫（6.10）I=xydmxy∫I=xzdm≈0xz∫I=yzdm≈0yz∫根据式6.9和式6.10，则有：II−−0φαsinxxxydddHI=•=−ωαII0xyyy00Iφαcos（6.11）zzddd=−(Iφααsin+++IiI)(αφαIsin)jI+φαcoskxxxyyyxyzz121 北京理工大学博士学位论文根据式6.5和式6.11，则有：PI=φφsin22α+−I(αααφsinαcos)αixxyy（6.12）+IIφφcos22α+(φαsinα++φαsinαφαcos)αzzxy[71]蜜蜂重0.0807g，蜜蜂一对翅的重量约为体重的0.5%，蜜蜂单翅的面积2−32S=26.89mm，假设翅的质量均匀分布，则蜜蜂单位面积翅质量dm=7.510kg/m×。当翅绕前缘旋转时，翅的惯性质量矩分别为−62I=610kgm×⋅xx−72I=×⋅5.4810kgmyy（6.13）−62I=×⋅6.5510kgmzz−62I=×⋅1.3410kgmxy6.3蜜蜂身体、双翅和单翅对力和功率的影响6.3.1单翅与双翅对力和功率的影响6.3.1.1不同最大拍动角下单、双翅模型的力和功率在本文研究内容中，均假设蜜蜂左翅和右翅的拍动是对称的。蜜蜂飞行时，研究左、右翅的相互作用是否影响气动力，需要分别建立单翅模型和双翅模型进行计算对比。图6.4为单翅和双翅的模型对比图，图6.4-1为单翅模型，图6.4-2为双翅模型，翅为无滑移固壁（wall），出口为自由出流（outflow）。如图6.4-2所示，假设左右翅的拍动完全对称，采用对称边界条件（symmetry）。选取标准k−ε湍流模型求解计算N-S方程。计算翅拍动5个周期的流场，选择第5个周期的气动力和功率曲线分析，−5时间步长为2.510s×。（1）（2）图6.4单、双翅模型的边界条件图6.5为不同最大拍动角φ下分别计算单翅模型和双翅模型得到的周期气动力max122 北京理工大学博士学位论文曲线。图6.5-1，图6.5-2，图6.5-3分别是φ为60°、90°和120°时的升力F和推力FmaxLT变化曲线图。计算中，翅的运动参数分别为：下拍中期攻角α=30°，拍动频率minfT=1/=200Hz，TTT==/2，ϑ=0。如图6.5所示，黑色线表示只计算单翅模0120型时的气动力，红色线表示计算双翅模型时一个翅上的气动力。在下拍第一个峰值处，双翅模型的升力略高于单翅模型的升力。在一周期内，单翅模型和双翅模型的升力FL和推力F周期曲线几乎重合。这表明，左、右翅的相互作用对翅的气动力影响不大。T（1）（2）（3）图6.5不同φ下的单、双翅模型的气动力曲线max表6.1为不同φ下单翅模型和双翅模型的平均升力和推力的对比结果。F和maxLsF表示计算单翅模型时的平均升力和推力，F和F表示计算双翅模型时的平均升TsLdTd力和推力，∆=−FFF，∆=−FFF。在不同φ下，单翅模型与双翅模型的LLsLdTTsTdmax平均升力差值不超过6%，平均推力差值均小于1%。这表明，双翅模型与单翅模型对气动力的影响是一致的，左右翅的交互作用对气动力的影响很小，几乎可以忽略。表6.1不同φ下单、双翅模型的平均气动力（mN）maxφmax/()°FLsFLdFTsFTd∆FL∆FT600.2370.2510.1720.1740.0140.002900.3880.3970.1980.1980.00901200.5900.5950.1830.1840.0050.001在计算单翅模型和双翅模型时，翅的质量惯性矩和运动参数均是相同的。如式6.12所示，惯性功率P只与翅的运动参数和惯性矩参数相关，与模型的边界条件无关，i因此单翅模型和双翅模型的惯性功率P是相同的，只需对比气动功率P曲线。图6.6-1，ia图6.6-2，图6.6-3分别是φ为60°、90°和120°时的气动功率变化曲线图。如图6.6max所示，在一周期内，单翅模型和双翅模型的气动功率曲线基本重合。已知单翅模型和123 北京理工大学博士学位论文双翅模型的气动力基本一致，因此对X、Y、Z轴的气动力矩也是一致的，根据式6.5，ddd气动功率只同气动力矩M和角速度ω有关，两个模型的对X、Y、Z轴的角速度ω也a相同，因此气动功率曲线几乎重合。（1）（2）（3）图6.6不同φ下的单、双翅模型的气动功率曲线max表6.2为不同φ下单翅模型和双翅模型的平均气动功率的对比结果。P表示计maxas算单翅模型时的平均气动功率，P表示计算双翅模型时的平均气动功率，ad∆=−PPP。单翅模型和双翅模型的平均气动功率差值不超过1%。因此，单翅模aasad型与双翅模型的气动功率是一致的。表6.2不同φ下单、双翅模型的平均气动功率（mW）maxφmax/()°PasPad∆Pa602.9102.9110.001903.5873.5730.0141204.8964.8580.0386.3.1.2不同中期攻角下单、双翅模型的力和功率分别计算当下拍中期攻角α为30°、40°和50°时，单翅模型和双翅模型的力和min功率。翅的其他运动参数为：最大拍动角φ=120°，拍动频率fT=1/=200Hz，max0TTT==/2，ϑ=0。图6.7为不同α下的一周期的力曲线。图6.7-1，图6.7-2，120min图6.7-3分别为α为30°、40°和50°时的升力和推力曲线图，黑色线表示计算单翅模min型的力，红色线表示计算双翅模型的力。如图6.7所示，在不同α下，单、双翅模min型的升力和推力曲线基本重合。这表明，左右翅的交互作用对翅气动力的影响可忽略不计。124 北京理工大学博士学位论文（1）（2）（3）图6.7不同α下单、双翅模型的气动力曲线min表6.3为不同α下单翅模型和双翅模型的平均气动力对比结果。如表6.3所示，min在相同α下，平均升力和平均推力差值均不超过2%。这表明，单翅模型与双翅模min型的气动力性能是一致的。表6.3不同α下单、双翅模型的平均气动力（mN）minαmin/()°FLsFLdFTsFTd∆FL∆FT300.5900.5950.1830.1840.0050.001400.6380.6460.0980.0990.0080.001500.7200.7220.0500.0590.0020.009单、双翅模型的运动参数和质量惯性矩均相同，因此单、双翅模型的惯性功率Pi在每个时刻均是一致的。图6.8为不同攻角α下单、双翅模型的气动功率P曲线图。mina图6.8-1，图6.8-2，图6.8-3分别是α为30°、40°和50°时的气动功率变化曲线图。min如图6.8所示，在一周期内，单翅模型和双翅模型的气动功率曲线基本重合。（1）（2）（3）图6.8不同α下的单、双翅模型的气动功率曲线表6.4为不同α下单翅和双翅的平均气动功率的对比结果。平均气动功率差值min不超过2%。因此，单翅模型与双翅模型的气动功率是一致的。125 北京理工大学博士学位论文表6.4不同α下单双翅的平均气动功率（mW）minαmin/()°PasPad∆Pa304.8964.8580.038406.1216.2310.110506.1316.0980.033综上所述，蜜蜂左右翅的交互作用对气动力的影响很小，求解单、双翅模型的气动力和功率是一致的，没有差别。6.3.2蜜蜂身体对气动力和功率的影响为研究蜜蜂身体对翅上的气动力的影响，建立完整蜜蜂模型，并与双翅模型的计算结果对比。图6.9为双翅和蜜蜂完整模型的对比图，图6.9-1为双翅模型（没有蜜蜂身体），图6.4-2为蜜蜂完整模型。如图6.9所示，均采用对称边界条件（symmetry），出口为自由出流（outflow），蜜蜂的翅和身体为无滑移固壁（wall）。计算翅拍动5个周期的流场，选择第5个周期的气动力和功率曲线分析。翅的运动参数为：中期攻角α=30°，最大拍动角φ=120°，拍动频率fT=1/=200Hz，TTT==/2，υ=0。minmax0120（1）（2）图6.9双翅模型和完整蜜蜂模型的边界条件如图6.10所示，黑线表示计算双翅模型时单个翅上的气动力和功率，红线表示计算完整蜜蜂模型时单个翅上的气动力和功率。图6.10-1为升力曲线，图6.10-2为推力曲线，图6.10-3为气动功率曲线。由于翅的质量惯性矩与运动参数相同，惯性功率是相同的，只需对比气动功率曲线。比较有蜜蜂身体和无蜜蜂身体两种情况，气动力曲线和气动功率曲线基本吻合，差别很小，可以忽略不计。这表明身体对于翅上的气动[214]力影响并不大，与于鑫的计算结果是一致的。126 北京理工大学博士学位论文（1）（2）（3）图6.10有身体和无身体的气动力和气动功率曲线图表6.5为计算双翅模型和蜜蜂整体模型的单个翅上的平均气动力和平均气动功率的对比结果。如表6.5所示，两种情况下，平均升力F、平均推力F和平均气动功LT率P的差别很小，可以忽略不计。因此在计算翅模型的力时，不需要建立蜜蜂身体模a型。表6.5有无身体的平均气动力和平均气动功率F/mNF/mNP/mWLTawithoutbody0.5950.1844.858withbody0.5910.1834.847蜜蜂身体对翅的气动性能影响可忽略不计，在6.3.1节中，已知单、双翅模型的气动性能基本一致，因此，只需建立单翅模型计算研究扑翼模型的气动力和功率。在以下研究翅模型的计算中，均只对单翅模型（如图6.4-1所示）进行计算。6.4典型算例分析翅模型的气动力和功率能耗，选取典型算例进行分析。在算例中，翅的运动参数为：最大拍动角φ=120°，下拍中期攻角α=30°，频率fT=1/=200H，maxmin0ZTTT==/2，ϑ=0。这些参数接近于实验中蜜蜂的振翅参数。计算振翅5个周期，120选取第5个周期的力和功率曲线分析。图6.11为翅模型一个周期的力和功率曲线图。图6.11-1为升力F和推力F周期LT曲线图，图6.11-2为惯性功率P和气动功率P的周期曲线图。ia如图6.11-1所示，在升力F曲线中，下拍和上拍的升力曲线变化趋势相似。下（上）L拍过程中均有两个升力F峰值，且下（上）拍末期的峰值高于下（上）拍初期的峰值。L下拍和上拍中期，F曲线较为平稳，在这个阶段，翅的攻角α保持不变。这表明，在L翅的平动阶段，能够稳定保持较高的升力值，Dickinson将这种现象称之为延迟失速127 北京理工大学博士学位论文[84]机制。如图6.11-1所示，在推力F曲线中，翅外旋（下拍末期和上拍初期）阶段，翅开T始加速外旋时，推力F增大，在下拍末期有F的峰值，随后推力F降低，当翅转变TTT为上拍后，推力F为负。下拍末期，攻角α要从30°变为90°，攻角的变化范围为60°；T上拍初期，攻角从90°变为120°，变化范围为30°，因此下拍末期推力F数值要大于T上拍初期推力F的数值。在一周期的其他阶段，下拍与上拍产生的推力F相互抵消。TT一周期的平均推力值主要由外旋阶段产生的推力贡献。在外旋阶段，推力F变化较为T明显，在内旋阶段，推力F变化并不明显。这是因为，推力沿着Y轴方向，翅内旋时，T旋转轴前缘几乎与Y轴平行，翅绕着Y轴旋转，沿着Y轴方向的力并不明显；翅外旋时，翅前缘与Y轴成一定角度，旋转时沿着Y轴方向有较大的分力，推力有明显变化。如图6.11-2所示，在惯性功率P曲线中，由于下拍运动和上拍运动完全对称，下i拍和上拍的惯性功率P曲线变化趋势是一致的。惯性功率P在下（上）拍初期和末期ii均有峰值，下（上）拍初期峰值大于下（上）拍末期的峰值，翅在平动阶段（攻角α保持不变），惯性功率P=0。在气动功率曲线中，气动功率P在一个周期中的值均为ia正，在下拍末期和上拍末期均会出现较高的峰值，这是由于翅突然加速旋转，此时升力有较高的峰值，翅的角速度也在增加，根据式6.6，气动功率P的值会剧增。当翅a由下（上）拍转为上（下）拍后，气动功率P降低。a（1）（2）图6.11一周期的力和功率曲线图图6.12为翅展向0.5R截面处的一周期涡量云图。图6.12为从翅尖向翅根的方向观察得到的云图，翅垂直纸面向外。如图6.12所示，无量纲涡量范围为Vor=−3~3。最大涡量值用红色表示，最小涡量值用蓝色表示，颜色代表涡量值的大小。128 北京理工大学博士学位论文在下拍开始阶段（t=0~0.05），翅绕主翅脉快速内旋，翅尾缘向后摆动，留下的尾迹携带正的涡量，翅被负涡量包围，此时对应升力为负。在t=0.075~0.125时，翅拍动速度增加开始平动，翅背面负的涡量向尾缘运动，积累产生启动涡，前缘涡开始生成并增大，如图6.12所示。此阶段升力F和推力F都增大。在t=0.1时刻附近出LT现F和F的第一个峰值，F的峰值更为明显，此时F=0.66mN、F=0.18mN，如LTLLT图6.11-1所示。在下拍中期阶段（t=0.15~0.375），翅保持攻角α=30°不变，翅只有平动没有转动。前缘涡逐渐增大并保持，从尾迹中可以观察到，前缘和后缘都在脱落二阶的小涡，主涡保持稳定，如图6.12所示。在此阶段，升力F和推力F平稳上升，如图6.11-1LT所示。这也验证了延迟失速机制。在t=0.4~0.425时，翅的前缘涡保持，由于拍动减速，尾迹中的负涡量在积累，如图6.12所示，升力F和推力F都增大。在下拍末期阶段（t=0.45~0.5），翅绕主LT翅脉快速后旋（外旋），翅后缘向前摆动，造成大面积负涡量尾迹区域，升力F和推L力F下降。在下拍结束的第二个升力F峰和推力F峰，时间上与翅的运动状态吻合，TLT此时F=1.79mN、F=2.14mN，如图6.11-1所示。LT在下拍过程中，升力F有两个峰值，第二个峰值明显高于第一个峰值。考虑翅的L转动方向，下拍开始时绕主翅脉内旋，对升力的贡献为负，下拍结束时绕主翅脉外旋，对升力贡献为正。在翅拍动过程中，上拍与下拍是完全对称的。如图6.11-1所示，上拍时的升力系数F曲线与下拍的F曲线基本一致。LLt=0.025t=0.05t=0.075t=0.1129 北京理工大学博士学位论文t=0.125t=0.15t=0.175t=0.2t=0.225t=0.25t=0.275t=0.3t=0.325t=0.35t=0.375t=0.4t=0.425t=0.45t=0.475t=0.5t=0.525t=0.55t=0.575t=0.6t=0.625t=0.65t=0.675t=0.7130 北京理工大学博士学位论文t=0.725t=0.75t=0.775t=0.8t=0.825t=0.85t=0.875t=0.9t=0.925t=0.95t=0.975t=1图6.12一周期内的涡量云图6.5力和功率的影响因素6.5.1频率蜜蜂在飞行时，能够根据外界环境改变翅的拍动频率。在风洞实验中，拍摄到蜜蜂的最高振翅频率约为190Hz。蜜蜂在自由飞行时，拍摄到的最高振翅频率约为250Hz。计算翅模型的拍动频率f分别为100Hz、110Hz、．．．（每间隔10Hz取值），300Hz时的气动力和功率。翅的其他运动参数为：最大拍动角φ=120°，下拍行程max的中期攻角α=30°，TTT==/2，ϑ=0。min120图6.13-1为平均升力F随频率f的变化曲线图。频率f增大，翅的振翅速度越L快，翅在一周期内的平均升力F越高。如图6.13-1所示，频率f越高，平均升力F的LL值越大。对其进行二阶拟合，平均升力F与频率f的关系为：L−−33−52F=×+×+×8.38101.0910ff1.5310（6.14）L图6.13-2为平均推力F随频率f的变化曲线图。如图6.13-2所示，频率f越高，T平均推力F的值越大。因为频率f增加，翅加速拍动，相同时间内，气流作用于翅上T131 北京理工大学博士学位论文的力越高，平均推力F越高（同平均升力增大的原理相同），对其进行二阶拟合，平T均推力F与频率f的关系为：T−−462F=−×+×0.045.8210ff8.3810（6.15）T图6.13-3为平均气动功率P和平均惯性功率P随频率f的变化曲线图。如图ai6.13-3所示，平均惯性功率P几乎为0，与频率f的改变无关。因为，在计算中，翅i模型的下拍运动与上拍运动是对称的，在一个周期中，惯性功率的均值为0。如图6.13-3所示，频率f增高，平均气动功率P越高。频率f增大，翅拍动速度加快，平均升力aF和平均推力F均有增加，气动力增大必然会消耗更多的气动功率，因此平均气动LT功率P随之增大。对平均气动功率P和频率f进行二阶拟合，则有：aa−42Pff=−0.410.034+×3.9810（6.16）a综上所述，频率f增高，平均升力F、平均推力F和平均气动功率P均增大，LTa且与频率f均符合二阶多项式拟合曲线。（1）（2）（3）图6.13平均力和功率随频率f的变化曲线6.5.2最大拍动角和中期攻角蜜蜂在振翅时，能够控制振翅的幅度，即能够通过改变最大拍动角调整自身的姿态和飞行状态。当最大拍动角φ分别为60°、80°、100°、120°、140°、160°和180°max时，计算翅的气动力和功率。翅的其他运动参数为：α=30°，频率f=200Hz，minTTT==/2，ϑ=0。图6.14为气动力随最大拍动角φ的变化曲线图。图6.14-1120max为在不同φ下，升力F和推力F的周期变化曲线，不同颜色的线代表不同φ，图maxLTmax6.14-2为平均升力F和平均推力F随φ的变化曲线，黑色线代表平均升力F值，LTmaxL红色线代表平均推力F值。T如图6.14-1所示，在升力F曲线中，相同时刻下，最大拍动角φ增加，升力FLmaxL增加。在相同振翅频率下，最大拍动角增大，翅在相同时间内拍动的角度增加，拍动速度增大，翅面的压力增大，因此最大拍动角变大，升力曲线整体上移。132 北京理工大学博士学位论文如图6.14-1所示，在推力F曲线图中，下拍初期和中期，相同时刻下，φ越大，Tmax推力F数值越高；上拍中（末）期的推力F数值变化与下拍中（初）期趋势相同，TT但是推力F的方向相反。在下拍末期，推力F有峰值，当φ=80°和φ=100°时，TTmaxmaxF峰值最高；当φ>°100时，φ越大，下拍末期的F峰值越低，出现峰值的时刻TmaxmaxT越早。在6.4章节典型算例的分析中曾提到，下拍末期推力F的峰值主要原因是翅的T旋转轴（翅前缘）与Y轴成一定夹角，当翅前旋与Y轴平行时，沿Y轴的分力（推力F）T最小，这也解释了当最大拍动角φ=180°时，下拍末期推力F峰值最低，在图6.14-1maxT中，当φ=80°和φ=100°时，推力F峰值最高，由此推论，当最大拍动角maxmaxTφ=90°，推力F的峰值应是最高的。maxT如图6.14-2所示，最大拍动角φ增大，平均升力F增加。在平均推力F曲线maxLT中，当φ=80°和φ=100°时F达到最大，这与之前分析的推力F周期曲线图的maxmaxTT结果是一致的。这表明，提高最大拍动角φ，能够有效地增加平均升力F。最大拍maxL动角φ=90°时，平均推力F最高。当最大拍动角φ=180°（趋于和Y轴平行）maxTmax时，平均推力F最低，几乎为0。因为当最大拍动角φ=180°时，上拍行程和下拍Tmax行程完全关于X轴对称，沿Y轴方向的分力几乎为0。（1）（2）图6.14不同最大拍动角下的气动力曲线图图6.15为功率随最大拍动角φ的变化曲线图。图6.15-1为在不同φ下，惯性maxmax功率P和气动功率P的周期变化曲线，不同颜色的线代表不同φ，图6.15-2为平均iamax惯性功率P和平均气动功率P随φ的变化曲线，黑色线代表平均气动功率P值，红iamaxa色线代表平均惯性功率P值。i如图6.15-1所示，在惯性功率P曲线中，在相同φ下，下拍阶段有两个P峰值，imaxi均在翅旋转阶段（第1个P峰值在下拍初期，第2个P峰值在下拍末期），在翅平动ii133 北京理工大学博士学位论文阶段，翅的惯性功率P=0。上拍的惯性功率曲线与下拍阶段的变化趋势相同。在相i同时刻下，最大拍动角φ增大，下拍初期惯性功率P峰值增高，下拍末期的惯性功maxi率P峰值降低，上拍阶段的惯性功率P曲线变化与下拍阶段相同。ii如图6.15-1所示，在气动功率P曲线中，在相同时刻下，最大拍动角φ增大，amax翅平动阶段的气动功率P增加，下拍初期和下拍末期的气动功率P峰值增高，上拍阶aa段的气动功率P变化与下拍阶段相同。a如图6.15-2所示，在不同最大拍动角φ下，平均惯性功率P均为0。这是由于maxi翅的上拍和下拍运动是对称的，攻角α与拍动角φ的变化在下拍和上拍阶段是对称的。根据惯性功率P计算公式，如式6.12所示，在1个周期内，平均惯性功率P=0。在ii平均气动功率P曲线中，最大拍动角φ增大，平均气动功率P增大。当φ增大，amaxamax平均气动力中，平均升力F明显增大，平均推力F的变化并不明显，如图6.14-2所LT示，随着φ越高，平均升力F较平均推力F越高，对平均气动功率P的影响越大。maxLTa因此，φ增大，翅的平均气动功率P增大。maxa（1）（2）图6.15不同最大拍动角下的功率曲线图蜜蜂振翅时，翅的扭转变形难以完全用数值计算方法模拟。在简化翅为单纯刚体运动时，为描述翅的翻转，引入攻角α。根据Dickinson油桶实验的机械扑翼模型和孙茂计算的仿生扑翼的运动特点，刚体翅在下拍和上拍中期阶段，攻角α保持不变，在下拍中期，攻角为α。min计算α分别为0、10°、20°、30°、……、90°时的气动力。翅的其他运动参数为：min最大拍动角φ=120°，频率f=200H，TTT==/2，ϑ=0。图6.16为不同α下maxZ120min的气动力曲线图。134 北京理工大学博士学位论文如图6.16-1所示，在升力曲线图中，翅上拍阶段和下拍阶段的升力F曲线变化相L同，只分析下拍阶段的升力F曲线。L当α≤°40时，在相同α下，下拍中期，升力F较为稳定，下拍初期和下拍minminL末期各有1个F峰值。在相同时刻下，下拍中期，α增大，升力F增大；在下拍的LminL第一个F峰值，α越大，升力F峰值越高；在下拍的第二个F峰值，α越小，LminLLmin升力F峰值越高。在翅旋转阶段，α值越低，升力F曲线振荡越明显。LminL当40°<α≤80°时，在相同α下，升力F较为稳定，翅旋转阶段并不见有明minminL显的F峰值。在相同时刻下，α增大，升力F降低；当α=90°时，升力F在一LminLminL周期内的数值基本为0。如图6.16所示，在推力F曲线图中，下拍初期和中期，相同时刻下，α越大，Tmin推力F越高；下拍末期，相同时刻下，α越大，推力F越低；上拍阶段，α越大，TminTmin推力F越低；相同α下，上拍中期和下拍中期的推力F由于方向相反，在一周期中，TminT几乎相互抵消，一周期内的平均推力F主要取决于翅外旋阶段的推力F。TT图6.16-2为平均升力F和平均推力F随α的变化曲线。当α=50°时，平均LTminmin升力F最大。α由0增大到50°时，平均升力F随之增大；α由50°增大到90°LminLmin时，平均升力F逐渐减小。α≤°50时，α增大，平均推力F降低，当α>°50LminminTmin时，平均推力F方向相反，且随着α增大，平均推力数值增大。Tmin（1）（2）图6.16不同α下的气动力曲线图min图6.17为不同α下的功率曲线图。图6.17-1为不同α下的惯性功率P和气动minmini功率P的周期曲线图，图6.17-2为随α变化的平均惯性功率P和平均气动功率P曲aminia线。在相同α下，P在下拍和上拍的变化是一致的，P也是一致的，因此只对翅模minia135 北京理工大学博士学位论文型在下拍过程的功率进行分析。如图6.17-1所示，在惯性功率P曲线中，在相同α下，惯性功率P在下拍和上imini拍的变化是一致的，只对下拍过程进行分析。下拍中期，由于攻角保持α不变，惯min性功率P=0。当α≤°40时，相同α下，下拍初期和下拍末期各有1个惯性功率Piminmini峰值；在相同时刻下，α越大，下拍初期的惯性功率P峰值越高，下拍末期惯性功mini率P峰值越小，P峰值出现的时刻基本一致。当α>°40时，同一α下，下拍初期iiminmin惯性功率P有峰值，下拍末期惯性功率P有谷值；在相同时刻下，α越大，下拍初iimin期P峰值和下拍末期P谷值出现的时刻越早，P峰值和P谷值数值变化不明显。iiii如图6.17-1所示，在气动功率P曲线中，相同时刻下，α越大，下拍中期的气amin动功率P越高，下拍初期和下拍末期的气动功率P越低。当α≤°40时，下拍末期aamin气动功率P有明显峰值，这是由于当α≤°40时，升力F和推力F在下拍末期均有aminLT较高峰值，因此气动功率P的值较高。a如图6.17-2所示，一周期的平均惯性功率P均为0。当α=30°时，平均气动功imin率P最低，当α>°30时，α增大，平均气动功率P增大。因为相同时刻下，αaminminamin增大，下拍和上拍中期的气动功率P增大。当α≤°30时，α减小，平均气动功率aminminP增大。因为相同时刻下，当α≤°30时，α越小，下拍和上拍末期气动功率的峰aminmin值越高。（1）（2）图6.17不同α下的功率曲线图min计算最大拍动角φ分别为60°、80°、100°、120°、140°、160°和180°时，α为maxmin0°、10°、20°、30°、40°、50°、60°、70°、80°和90°时翅模型的气动力和功率。翅的其他运动参数为：拍动频率fT=1/=200Hz，TTT==/2，ϑ=0。图6.18为不同0120136 北京理工大学博士学位论文φ和α下的平均气动力和平均气动功率示意图。当翅下拍和上拍运动完全对称时，maxmin平均惯性功率P=0，因此只给出平均气动功率P云图。图6.18-1为平均升力F云图，iaL图6.18-2为平均推力F云图，图6.18-3为平均气动功率P云图。如图6.18所示，每Ta幅图最上方为平均力或平均功率值的范围，最小数值用蓝色表示，最大数值用红色表示，纵坐标为φ，变化范围为60~180°°，横坐标为α，变化范围为0~90°°。图maxmin上每个圆圈代表对应于φ和α下的平均气动力或平均气动力的值。maxmin如图6.18-1所示，在相同α下，最大拍动角φ增大，平均升力F增加；在相minmaxL同最大拍动角φ下，α在50°附近处平均升力F最大，在0°和90°时平均升力F最maxminLL小。如图6.18-2所示，当最大拍动角φ≤°140，平均推力F受α的影响较大；在maxTmin相同φ下，α越小，平均推力F越大；当最大拍动角φ>°140，α对平均推maxminTmaxmin力F的影响减弱；平均推力F随α减小而增大。TTmin如图6.18-3所示，在相同α下，最大拍动角φ越大，平均气动功率P越高。minmaxa相同φ下，当α在20~30°°时，平均气动功率最低；低于这个范围时，α越低，maxminmin平均气动功率P越高；高于这个范围时，α越大，平均气动功率P越高。因此，当aminaα=90°，φ=180°时，平均气动功率P最高。降低平均气动功率P，可以降低φ，minmaxaamax将α调整在20~30°°左右。但是降低φ，会使平均升力F降低。minmaxL（1）（2）（3）图6.18平均气动力和平均气动功率随φ和α变化的云图maxmin6.5.3T1计算当T分别为0.2T0、0.4T0、0.6T0、0.8T0和T0时的气动力和功率。翅的其他运1动参数为：最大拍动角φ=120°，下拍中期攻角α=30°，TT=/2，拍动频率maxmin20fT=1/=200Hz，ϑ=0。0137 北京理工大学博士学位论文图6.19为不同T下的气动力曲线图。图6.19-1为不同T下的升力F和推力F周11LT期曲线图。图6.26-2为平均升力F和平均推力F随T的变化曲线图。LT1如图6.19-1所示，在升力F曲线中，T值越大，下拍中期相同时刻下的升力F越L1L高，但是下拍末期升力F的峰值越低。在推力F曲线中，T主要影响下拍末期的推LT1力F。T越小，下拍末期推力F峰值越高。如图6.19-2所示，T增大，平均升力F增T1T1L大，平均推力F降低。但是T对平均升力F和平均推力F的改变较小。T1LT（1）（2）图6.19不同T下的气动力曲线1图6.20为不同T下的功率曲线图。图6.20-1为不同T下的惯性功率P和气动功率11iP周期曲线图。图6.20-2为平均惯性功率P和平均气动功率P随T的变化曲线图。aia1如图6.20-1所示，在惯性功率P曲线中，T值越大，拍动速度变化的时间所占周i1期比重越高，惯性功率P=0的时间就越少。当TT=时，惯性功率P在1周期内均变i10i化。当TT=/2时，在旋转阶段，惯性功率P的变化剧烈，有较高的峰值和较低的谷10i值。如图6.20-1所示，在气动功率P曲线中，相同时刻下，T值越大，下拍中期的气a1动功率P越大，下拍末期的气动功率P越小。上拍时气动功率P的变化趋势与下拍阶aaa段相同。如图6.20-2所示，平均惯性功率P为0。当TT<0.6时，T值越小，平均气动功i101率P越高，当TT≥0.6时，T值越大，平均气动功率P越高。T对平均气动功率P的a101a1a影响较小。138 北京理工大学博士学位论文（1）（2）图6.20不同T下的气动力曲线16.5.4T2翅在旋转阶段所需时间为T，通过改变T的值，可以控制翅模型的旋转时间。计22算当T分别为0.3T0、0.4T0、0.5T0、0.6T0、0.7T0和0.8T0时的气动力和功率。翅的其2他运动参数为：最大拍动角φ=120°，下拍中期攻角α=30°，TT=，ϑ=0，maxmin10拍动频率fT=1/=200Hz。0图6.21为气动力随T的变化曲线图。图6.21-1为在不同T下，升力F和推力F22LT的周期变化曲线，图6.21-2为平均升力F和平均推力F随T的变化曲线。LT2（1）（2）图6.21不同T下的气动力曲线图2如图6.21-1所示，在升力F曲线中，T越小，翅模型旋转所用时间越短，旋转L2阶段翅的转动速度越快，产生的升力F峰值越高。T越大，旋转时间越长，翅在旋转L2阶段产生的升力F越趋于缓和。L139 北京理工大学博士学位论文如图6.21-1所示，在推力F曲线中，T只对外旋阶段（下拍末期和上拍初期）T2的推力F影响较大。T越小，旋转时间越短，下拍末期F曲线振荡越剧烈，F的峰T2TT值越高。如图6.21-2所示，平均升力F随T增大而减小，当TT/≥0.6时，平均升力F的L220L变化较小。在计算算例中，平均推力F均为负，即推力均沿着Y轴负向，这表明推力T能够推动蜜蜂向前飞行（如图6.1所示）。T增大，平均推力F降低。这表明当翅模2T型所用旋转时间越少，平均升力F和平均推力F均会增大。LT图6.22为功率随T的变化曲线图。图6.22-1为在不同T下，惯性功率P和气动22i功率P的周期变化曲线，图6.22-2为平均惯性功率P和平均气动功率P随T的变化aia2曲线。如图6.22-1所示，在惯性功率P和气动功率P曲线中，T只影响翅旋转阶段的功ia2率曲线。在惯性功率P曲线中，T越小，翅模型旋转速度越快，旋转阶段翅的惯性功i2率振荡越剧烈。在外（内）旋时，下（上）拍末期P峰值越高，上（下）拍初期P谷ii值越低。在气动功率P曲线中，T越小，翅模型的转动速度越快，下拍和上拍末期的a2气动功率P峰值越高。这表明翅模型旋转速度越快，能耗越高。a（1）（2）图6.22不同T下的功率曲线图26.5.5翻转时刻计算不同翻转时刻（ϑ=-1、0和1）情况下，力和功率的变化情况。翅的其他运动参数为：最大拍动角φ=120°，下拍中期攻角α=30°，频率f=200H，TT=，maxminZ10TT=/2。20140 北京理工大学博士学位论文图6.23为不同翻转时刻下力和功率的周期变化曲线图。图6.23-1为不同翻转时刻下升力F和推力F的周期曲线图，图6.23-2为不同翻转时刻下惯性功率P和气动功LTi率P的周期曲线图。a如图6.23-1所示，在升力F曲线中，提前翻转（ϑ=−1）和对称翻转（ϑ=0）L时，在外旋和内旋阶段，升力F均有明显峰值。提前翻转时，外旋阶段F峰值约为LL3.8mN；对称翻转时，外旋阶段F峰值约为2mN。即提前翻转的F峰值明显高于对LL称翻转的F峰值。延迟翻转（ϑ=1）没有升力峰值，周期升力F曲线较为平缓。由LL升力曲线推断，提前翻转能够提高平均升力F。L如图6.23-1所示，在推力F曲线中，提前翻转（ϑ=−1）和对称翻转（ϑ=0）T时，在外旋阶段，推力F均有明显峰值。提前翻转的推力F峰值高于对称翻转时的TT推力F峰值。延迟翻转（ϑ=1）时，推力F曲线较为平缓。由推力曲线推断，提前TT翻转能够提高平均推力F。T如图6.23-2所示，在惯性功率P曲线中，延迟翻转（ϑ=1）的P曲线较为平缓。ii提前翻转（ϑ=−1）的P曲线在旋转阶段振荡剧烈，P先上升到峰值，再下降到谷值。ii旋转阶段惯性功率P振荡的正负值相反而抵消，因此对平均惯性功率P的影响很小。ii如图6.23-2所示，在气动功率P曲线中，翅模型翻转时刻越提前，旋转阶段气动a功率P的峰值越高。当ϑ=1时，气动功率P曲线在一周期内较为平缓。由此可推断，aa当翅提前翻转时，所需平均气动功率P值较高。a（1）（2）图6.23不同翻转时刻下的力和功率曲线表6.6为在不同ϑ下的平均力和平均功率数值。如表6.6所示，提前翻转（ϑ=−1）时的平均升力F（0.85mN）最大，比对称翻转（ϑ=0）时平均升力F（0.71mN）LL141 北京理工大学博士学位论文大20%。延迟翻转（ϑ=1）的平均升力F（0.34mN）最小，比对称翻转时平均升力LF小52%。ϑ能明显改变翅模型的平均推力F。ϑ=−1时的平均推力F（0.44mN）LTT约为ϑ=0时平均推力F（0.05mN）的8.8倍。延迟翻转时，翅的平均推力F会降低，TT甚至改变平均推力F的方向。翻转时刻对平均惯性功率P影响较小，平均惯性功率PTii几乎为0。翻转时刻对平均气动功率P的影响很大，提前翻转会增加平均气动功率P，aa延迟翻转会降低平均气动功率P。ϑ=−1时的平均气动功率P（13.07mW）比ϑ=0时aa平均气动功率P（6.82mW）大90%，ϑ=1时的平均气动功率P（2.82mW）比ϑ=0aa时平均气动功率P小59%。综上所述，提前翻转时刻，能够增加平均升力F、平均aL推力F和平均气动功率P。Ta表6.6不同翻转时刻下的平均力和平均功率值ϑF/(mN)F/(mN)P/(mW)P/(mW)LTia-10.850.440.00113.0700.710.05-0.0596.8310.34-0.030.0022.826.5.6旋转轴位置在前面章节的数值模拟计算中，未加说明情况下，翅在旋转时均绕翅前缘主翅脉旋转，如图6.24所示，即绕k=0的轴旋转。计算旋转轴距前缘k（kc/值取0、0.1、0.2、．．．、0.9）时的气动力和功率，如图6.24所示，c为翅的弦长。翅的运动参数为：最大拍动角φ=120°，中期攻角maxα=30°，拍动频率fT=1/=200Hz，TTT==/2，ϑ=0。min0120图6.24旋转轴位置示意图图6.25为翅在不同旋转轴下的气动力曲线，图6.25-1为不同旋转轴下的升力F和L推力F周期曲线，图6.25-2为平均升力F和平均推力F随旋转轴位置的变化曲线。TLT142 北京理工大学博士学位论文（1）（2）图6.25不同旋转轴下的气动力曲线图如图6.25-1所示，同一旋转轴下，升力F曲线在下拍和上拍的变化趋势是相同的，L只对下拍的升力F曲线进行分析。当kc/≥0.4时，下拍初期升力F先上升达到1个LL峰值，后下降到1个谷值，距离前缘位置越远（k越大），升力峰值越高，谷峰也越低。当kc/<0.4时，下拍初期升力F先小幅下降再上升。在下拍中期，相同时刻下，不L同旋转轴下的升力相差不大。在下拍末期，相同时刻下，旋转轴距离前缘位置越近（k越小），升力F的数值越高。L如图6.25-1所示，在推力F曲线中，在下拍初期和中期，上拍中期和末期，翅的T旋转轴不同，相同时刻下，翅的推力F相差不大。在外旋阶段（下拍末期和上拍初期），T旋转轴对推力F的影响较为明显；在相同时刻下，旋转轴越靠近翅的前缘（kc/越小），T推力F越高。T如图6.25-2所示，旋转轴距离前缘位置越近（kc/越小），平均升力F和平均推L力F的值越高，当旋转距离前缘位置越远（kc/越大），平均升力F和平均推力F的TLT值越低。在Dickinson对果蝇翅模型的油桶实验中，果蝇翅模型的旋转轴位置越靠近[84]前缘，旋转环量值越高，产生的力越大。因此旋转轴越靠近前缘，平均气动力越高。图6.26为翅在不同旋转轴下的功率曲线，图6.26-1为不同旋转轴下的惯性功率Pi和气动功率P的周期曲线，图6.26-2为平均惯性功率P和平均气动功率P随旋转轴aia位置的变化曲线。如图6.26-1所示，旋转轴的位置只对气动功率P有影响，旋转轴的位置对惯性功a率P没有影响，不同旋转轴的惯性功率P曲线重合。在气动功率曲线中，旋转轴的位ii置对旋转阶段的气动功率P影响较大，对平动阶段的P几乎没有影响。在旋转阶段，aa143 北京理工大学博士学位论文当kc/<0.4时，旋转轴越靠近前缘（kc/越小），相同时刻下的气动功率P越高；当akc/≥0.4时，旋转轴距前缘越远（kc/越大），相同时刻下的气动功率P越高。a如图6.26-2所示，平均惯性功率P=0，旋转轴的位置主要影响平均气动功率P。ia当kc/=0.4时，平均气动功率P最低。当kc/<0.4时，旋转轴距翅前缘越远，平均a气动功率P越低；当kc/≥0.4时，旋转轴距翅前缘越近，平均气动功率P越低。aa（1）（2）图6.26不同旋转轴位置下的平均气动力系数图6.6本章小结本章基于实体蜜蜂翅，建立了翅的坐标系，分析了翅的升力、推力、气动功率和惯性功率。通过数值模拟单翅、双翅和整体蜜蜂模型的气动特性，确定了左、右翅的相互作用和蜜蜂身体均对翅的气动性能影响很小，因此，只需利用单翅模型研究气动特性和功率能耗。通过升力、推力、气动功率、惯性功率曲线和涡量云图，分析了翅运动一周期的气动特性。计算并分析了频率f、最大拍动角φ、中期攻角α、T、maxmin1T、翻转时刻和旋转轴位置对气动力和功率的影响。2平均升力F主要受f、φ、α、翻转时刻和旋转轴位置的影响，T和T对FLmaxmin12L的影响可忽略；增大f和φ，F升高；α在50°附近处时，F最高；提前翻转时maxLminL刻，F升高；旋转轴位置越靠近前缘，F越高。φ和T对平均推力F的影响较小，LLmax1TF主要受f、α、T、翻转时刻和旋转轴位置的影响。增大f、降低α和T、延Tmin2min2迟翻转、将旋转轴位置靠近前缘，能提高F。T不同影响因素对翅的平均惯性功率P的影响很小，可忽略不计。平均气动功率Pia主要受f、φ、α、T、翻转时刻和旋转轴位置的影响，T对P影响可忽略不计。maxmin21a144 北京理工大学博士学位论文增大f和φ、降低T、将翻转时刻提前，均能提高P；当α=°°20~30时，P最max2amina低，当α=90°时，P最高；当旋转轴位于前缘与尾缘中间时，P最低，旋转轴在minaa前缘或尾缘处时，P最高。a145 北京理工大学博士学位论文结论与展望本文结论扑翼飞行被自然界中的鸟类和昆虫广泛所采用，被认为是生物进化的最优飞行方式。微型扑翼飞行器具有体积小、重量轻、成本低等优势，有良好的应用前景。基于仿生学原理的扑翼飞行研究是当前的热点研究问题，国内的研究主要集中在仿生数值模拟计算和扑翼样机的研制两个方面，然而对于自然界中真实昆虫的飞行机理还尚未展开系统的研究。由于蜜蜂的身体尺寸和飞行效率都与理想的微型扑翼飞行器的特征相似，本文以真实蜜蜂做为研究对象，结合流场显示实验和数值仿真计算对蜜蜂振翅的流场进行了研究。本文的主要工作可总结如下：（1）设计制作了专门用于昆虫流场显示的低速烟风洞，将风洞内的速度分布作为流场品质的评价指标，利用数值计算方法研究了风洞实验段、扩散段和收缩段对风洞流场品质影响的规律。风洞为开路闭口式，整体结构包括喇叭口、稳定段、收缩段、实验段、扩散段和动力段；在实验段，研究了实验段截面边长和实验段长度对中分截面处速度分布的影响；在扩散段，研究了扩散段长度对扩散段出口截面和入口截面速度分布的影响；在收缩段，分别研究了维多辛斯基曲线、双三次曲线和五次方曲线三种收缩曲线线型对风洞流场品质的影响，通过对比分析风洞中分截面及收缩出口截面的速度分布，将维多辛斯基曲线做为收缩曲线线型；模拟计算研究了不同收缩比下风洞中分截面和收缩段出口截面的速度分布，确定了最佳收缩比。根据数值模拟结果，确定了风洞各段的最佳尺寸和整体结构参数。（2）为提高风洞内的流场品质，保证气流能够平滑进入风洞，降低湍流度，使气流均匀平稳，设计了风洞的整流装置，在稳定段前部设计了弧形喇叭口，在稳定段和动力段内添加了蜂窝器和阻尼网；由于风洞主要用于流场显示实验，烟的质量和效果均对实验的成败至关重要，因此通过对比不同材质、形状的电阻丝和发烟油的发烟效果，专门设计和调试了发烟装置，通过调整合适的电压值控制烟的浓度和清晰度；为从两个角度拍摄蜜蜂，对高速摄像机和灯光进行了合理的实验布局，开发了用于研究小型昆虫扑翼飞行的完整的流场显示实验系统。（3）利用实验研究了蜜蜂的振翅规律和涡的演变规律。首先，利用两台高速摄146 北京理工大学博士学位论文像机同步拍摄蜜蜂振翅运动，通过建模和图像分析，得到了蜜蜂翅拍动角的变化规律，拍动角周期变化曲线近似为半个周期的正弦曲线。然后，对蜜蜂进行了远场和近场流场显示实验，通过分析烟线实验图，提出了蜜蜂振翅时涡的演变规律：在外旋阶段，尾缘处有下拍停止涡和上拍启动涡，下拍停止涡对升力的贡献为负，上拍启动涡对升力的贡献为正；在旋转阶段，翅有较强的扭转变形，使得前翅与后翅尾缘处的涡量方向相反；下拍末期，存在前缘涡，根据不同展向位置的烟线实验图片，由翅根至翅尖，前缘涡的尺寸逐渐增大，并与翅尖涡连接在一起，形成三维锥形的前缘涡结构；蜜蜂振翅的下拍末期，在每个翅上，前缘涡与翅尖涡涡管和翅根涡尾迹共同组成了涡环，且左、右翅的涡环互相独立。（4）根据蜜蜂的振翅规律给出了翅的运动方程，为了数值模拟研究翅气动性能的影响因素，对运动方程进行了优化；利用动网格技术求解了N-S方程，验证了本文数值模型的可靠性，并求解了蜜蜂振翅的流场。数值仿真的流场特征与烟线实验所提出的涡的演变规律相吻合，包括外旋产生的下拍停止涡与启动涡、前缘涡的三维结构和下拍末期的涡环。此外，根据数值计算结果，得出以下结论：在上、下拍的中期和末期，均存在前缘涡，翅的旋转使得前缘涡在上、下拍末期尺寸达到最大；上、下拍动初期，尾缘涡附着于翅的最宽截面尾缘处；在上、下拍动末期，存在翅根涡；在上、下拍末期，翅尖涡尾迹涡管、翅根涡尾迹涡量和翅的边界涡量均能形成涡环，下拍末期的涡环呈圆形，上拍末期的涡环为耳型，在上、下拍末期的涡环中部均会产生垂直于涡环平面的射流。（5）基于实体蜜蜂翅模型及翅的坐标系，数值计算研究了左、右翅的相互作用和蜜蜂身体对翅的气动力和功率的影响，结果表明其影响均可忽略不计；针对蜜蜂单翅模型，数值计算研究了影响翅的气动力和功率的因素，包括频率、最大拍动角、中期攻角、翻转时刻以及旋转轴位置等。数值结果表明，这些因素对平均惯性功率的影响较小，主要影响翅的平均升力、平均推力和平均气动功率；为提高平均升力，可以增大频率和最大拍动角、调整下拍中期攻角在50°附近、将翻转时刻提前、将旋转轴靠近前缘；为提高平均推力，可以增大频率、降低下拍中期攻角、降低翅的旋转时间（降低T）、延迟翻转时刻、将旋转轴靠近前缘；为提高平均气动功率，可以增大频2率和最大拍动角、调整下拍中期攻角远离30°、降低翅的旋转时间、提前翻转时刻、调整旋转轴靠近前缘或尾缘。147 北京理工大学博士学位论文创新点本文的创新点主要有以下几方面：（1）为了进行小型昆虫的流场显示实验，设计并制作了低速烟风洞，通过数值计算研究了风洞各段尺寸对流场品质的影响规律，并确定了风洞的最佳尺寸和整体结构。在此基础上，设计了风洞的整流装置和发烟装置，搭建了完整的实验平台，开发了流场显示实验系统。（2）对实体蜜蜂进行了风洞流场显示实验，利用高速摄像机从两个角度拍摄了蜜蜂振翅过程和烟线的变化情况。通过对图像进行三维空间分析，得到了拍动角的变化规律；分析图像中烟线的变化规律，得到了蜜蜂周期振翅的远场和近场的流场特征，提出了蜜蜂振翅时涡的演变规律。（3）基于蜜蜂的振翅规律给出了翅的运动方程，并对运动方程进行了优化；利用动网格技术求解了N-S方程，并验证了数值计算的可靠性；通过数值计算求解了蜜蜂模型振翅时的三维流场特征，计算结果一方面验证了风洞实验所提出的涡的演变规律，另一方面表明所建立的蜜蜂数值仿真模型能够用于研究扑翼飞行的机理和气动性能。（4）将理论分析和数值模拟相结合，分析了影响翅气动力和功率的多种因素，包括振翅频率、最大拍动角、中期攻角、T1、T2、翻转时刻和旋转位置。分析结果对微型扑翼机的研制提供了指导意义。展望本文对蜜蜂的扑翼问题进行了研究，在风洞流场显示实验基础上，分析了蜜蜂的振翅规律及流场特征进行了分析，提出了蜜蜂振翅的涡的演变规律；利用数值仿真模拟了蜜蜂飞行的流场，验证了实验提出的三维涡结构和涡的演变规律；计算并分析了影响翅气动力和功率的因素，得到了具有一定价值的结论。但由于实验技术和实验条件的限制，仍有许多问题有待于进一步研究和深化，需要在以后的工作中逐渐完善。以下有待进一步研究的几个问题：（1）本文将蜜蜂身体固定，在烟风洞中研究了蜜蜂的周期振翅规律和流场特征，分析结果表明蜜蜂通过振翅可以产生较高升力，左、右翅的涡环结构可能是蜜蜂飞行操纵灵活性的重要原因。由于蜜蜂的飞行具有灵活多变性，还需研究蜜蜂在多种自由飞行状态下（例如起飞、降落和拐弯等）的流场结构。148 北京理工大学博士学位论文（2）本文基于烟线法研究了蜜蜂振翅的流场特征，实验成本低，经济性好，且流场显示效果良好，能直观地显示并分析振翅引起的涡的演变规律。为能够定量地测量真实蜜蜂振翅时周围的流场，可采用PIV技术，通过测量流场的速度变化，进一步研究蜜蜂振翅时周围流场涡强度的变化规律。（3）蜜蜂翅在拍动过程中，除了刚性转动运动，还伴随着复杂的变形运动，振翅周期过程中，不同时刻下，翅面上不同区域的变形程度均不相同。本文分析了流场显示实验结果中蜜蜂翅扭转变形对前、后翅尾缘涡涡量方向的变化规律，在数值模拟计算中，研究了刚体翅模型对流场结构的影响。下一步工作可以根据蜜蜂翅的可变形性，引入生物材料的粘弹性本构模型，采用流固耦合方法来研究蜜蜂扑翼飞行的气动性能。149 北京理工大学博士学位论文参考文献[1]陈国栋,贾培发,刘艳.微型飞行器的研究与发展[J].机器人技术与应用,2006,2:34-44.[2]StevenA.Palm-sizeSpyPlane[J].MechanicalEngineering,1998,11(3):74-78.[3]孙瑜,张杰,刘虎,武哲.MAV微型飞行器研究进展与总体设计[J].飞机设计,2010,30(6):11-16.[4]曲东才.微型无人机研制的关键技术及军事应用[J].飞机设计,2007,3:46-51.[5]DaltonS.BorneontheWind:TheExtraordinaryWorldofInsectsinFlight[M].Reader’sDigestPress,1975.[6]李占科,宋笔锋,宋海龙.微型飞行器的研究现状及其关键技术[J].飞行力学,2003,21(4):1-4.[7]杨淑利.适用于突风环境的微型扑翼飞行器柔性翼气动特性研究[D].西北工业大学硕士论文.2007.[8]袁昌盛,付金华.国际上微型飞行器的研究进展与关键问题[J].航空兵器,2005,6:50-53.[9]WoodsMI,HendersonJF,LockGD.EnergyRequirementsfortheFlightofMicro-airVehicles[J].AeronautJ,2001,105(2546):135-149.[10]Pornsin-sirirakTN,TaiYC,NassefH,HoCM.TitaniumalloyMEMSWingTechnologyforaMicroaerialVehicleApplication[J].SensorsandActuatorsA:Physical,2001,89:95-103.[11]Pornsin-sirirakTN,TaiYC,NassefH,HoCM.WingTechnologyforaBattery-poweredOrnithopter[C].MEMS,Miyazaki,2000.[12]杨俊.飞行原理[M].西南交通大学出版社,2012.th[13]AndersonJD.FundamentalsofAerodynamics4Ed[M].McGraw-Hill,2007.[14]ClancyLJ.Aerodynamics[M].CambridgeUniversityPress,1975.nd[15]KatzJ,PlotkinA.Low-speedAerodynamics2Ed[M].CambridgeUniversityPress,2001.[16]WhiteFM.ViscousFluidFlow[M].McGraw-Hill,1991.[17]OlMV,BernalL,KangCK,ShyyW.ShallowandDeepDynamicStallforFlappingLowReynoldsNumberAirfoils[J].ExpFluids,2009,46:883-901.[18]ShyyW,LianY,TangJ,ViieruD,LiuH.AerodynamicsofLowReynoldsNumberFlyers[M].CambridgeUniversityPress,2008.150 北京理工大学博士学位论文[19]MuellerTJ,DeLaurierJD.AerodynamicsofSmallVehicles[J].AnnualReviewofFluidMechanics2003,35:89-111.[20]崔尔杰.生物运动仿生力学与智能微型飞行器[J].力学与实践,2004,26(2):1-8.[21]兰世隆,孙茂.模型昆虫翼作非定常运动时的气动力特性[J].力学学报,2001,33(2):173-182.[22]孙茂,黄华.微型飞行器的仿生力学--蝴蝶飞行的气动力特性[J].北京航空航天大学学报,2006,32(10):1146-1151.[23]SunM,LanSL.AComputationalStudyoftheAerodynamicForcesandPowerRequirementsofDragonfly(Aeshnajuncea)Hovering[J].JExpBiol,2004,207:1887-1901.[24]SunM,TangJ.UnsteadyAerodynamicForceGenerationbyaModelFruitFlyWinginFlappingMotion[J].JExpBiol,2002,205:55-70.[25]SunM,WuJH.AerodynamicForceGenerationandPowerRequirementsinForwardFlightinaFruitFlywithModeledWingMotion[J].JExpBiol,2003,206:3065-3083.[26]孙茂,熊燕.微型飞行器的仿生力学-蜜蜂悬停飞行的动稳定性研[J].航空学报,2005,26(4):385-391.[27]SunM,WangJK.FlightStabilizationControlofaHoveringModelInsect[J].JExpBiol,2007,210:2714-2722.[28]吴江浩,王济康,孙茂.蜜蜂悬停飞行控制的仿生力学研究[J].航空学报,2007,28(4):776-782.[29]SunM,XiongY.DynamicFlightStabilityofaHoveringBumblebee[J].JExpBiol,2005,208:447-459.[30]SunM,TangJ.LiftandPowerRequirementsofHoveringFlightinDrosophilaVirilis[J].JExpBiol,2002,205:2413-2427.[31]LiuY,SunM.WingKinematicsMeasurementandAerodynamicsofHoveringDroneflies[J].JExpBiol,2008,211:2014-2025.[32]MouXL,LiuYP,SunM.WingMotionMeasurementandAerodynamicsofHoveringTrueHoverflies[J].JExpBiol,2011,214:2832-2844.[33]DuG,SunM.EffectsofWingDeformationonAerodynamicForcesinHoveringHoverflies[J].JExpBiol,2010,213:2273-2283.[34]刘岚,方宗德,张西金.微型扑翼飞行器的升力风洞试验[J].航空动力学报,2007,22(8)：1315-1319.151 北京理工大学博士学位论文[35]张西金,方宗德,杨小辉.一种仿蜜蜂类昆虫扑翼悬停控制的仿真估算研究[J].航空动力学报,2008,23(9):1561-1567.[36]杨智春,李思政,舒忠平,白存儒.一种柔性微型扑翼设计及其气动力特性的试验研究[J].机械科学与技术,2008,25(1):12-14.[37]谢辉,宋文萍,宋笔锋.微型扑翼绕流的N-S方程数值模拟[J].西北工业大学学报,2008,26(1):104-109.[38]熊超,宋笔锋.微型扑翼飞行器机理风洞试验研究[J].科学技术与工程,2007,7(11):2576-2580.[39]熊超,宋笔锋,袁昌盛,邵立民,张亚锋.微型扑翼飞行器机翼纵向力矩特性风洞试验研究[J].西北工业大学学报,2007,25(5):733-736.[40]侯宇,方宗德,刘岚,傅卫平.仿生微扑翼飞行器机构动态分析与工程设计方法[J].航空学报,2005,26(2):173-178.[41]昂海松,曾锐,段文博,史志伟.柔性扑翼微型飞行器升力和推力机理的风洞试验和飞行试验[J].航空动力学报,2007,22(11):1838-1845.[42]曾锐.仿鸟微型扑翼飞行器的气动特性研究[D].南京航空航天大学博士学位论文,2004.[43]肖天航,段文博,昂海松.仿鸟柔性扑翼气动特性与能耗的数值研究[J].空气动力学学报,2011,29(6):709-718.[44]肖天航,昂海松,周新春.柔性扑翼非定常流场的数值计算方法[J].航空学报,2009,30(6):990-999.[45]孟令兵,昂海松,肖天航.基于CFD/CSD的蜻蜓柔性翼气动特性分析[J].航空动力学报,2014,29(9):2063-2069.[46]肖天航,昂海松.蜻蜓前后翼相位特性和气动干扰的数值研究[J].航空学报,2009,30(7):1165-1175.[47]RaynerJMV.AVortexTheoryofAnimalFlight.Part1.TheVortexWakeofaHoveringAnimal[J].JFluidMech,1979,91:697-730.[48]WangXX,WuZN.LiftForceReductionDuetoBodyImageofVortexforaHoveringFlightModel[J].JFluidMech,2012,709:648-658.[49]WangXX,WuZN.Stroke-averagedLiftForcesDuetoVortexRingsandTheirMutualInteractionsforaFlappingFlightModel[J].JFluidMech,2010,654:453-472.152 北京理工大学博士学位论文[50]曾理江,宋德强,郝群.昆虫运动机理的研究[J].光学技术,1999,6:18-25.[51]郝群,曾理江,河内启二.昆虫翅膀的运动参数测量研究[J].北京理工大学学报,2003,23(4):444-448.[52]WangH,ZengLJ,LiuH,YinCY.MeasuringWingKinematics,FlightTrajectoryandBodyAttitudeDuringForwardFlightandTurningManeuversinDragonflies[J].JExpBiol,2003,206:745-757.[53]陈文元,张卫平.微型扑翼式仿生飞行器[M].上海交通大学出版社,2010.[54]左德参.仿生微型飞行器若干关键问题的研究[D].上海交通大学博士学位论文,2007.[55]左德参,陈文元,彭松林,李智军,张卫平.基于MEMS技术的扑翼式微飞行器的研究[J].系统仿真学报,2005,17(6):1505-1508.[56]彭松林,左德参,张卫平,陈文元.曲柄摇杆扑翼机构的联合仿真及优化设计[J].系统仿真学报,2007,19(8):1758-1761.[57]陆夕云,杨基明,尹协振,童秉纲.飞行和游动的生物运动力学和仿生技术研究[J].中国科学技术大学学,2007,27(10):1159-1163.[58]鲍麟,童秉纲.模型翼拍动中动态柔性变形效应的数值研究[J].中国科学院研究生院学报,2005,22(6):676-684.[59]余永亮,童秉纲.拍翼轨迹对昆虫前飞气动性能的影响[C].第十届全国分离流、旋涡和流动控制会议论文集,2005,155-159.[60]DemollR.DerFlugderInsektenundVögel[M].Jena,GFischer,1918.[61]HoffW.DerFlugderInsektenundderVögel[J].DieNaturWissenschaften,1919,10:162-169.[62]DemollR.ZuschriftenandieHerausgeber.DerFlugderInsektenundderVögel[J].DieNaturWissenschaften,1919,27:480-482.[63]Weis-FoghT.FlappingFlightandPowerinBirdsandInsects,ConventionalandNovelMechanisms[J].InSwimmingandFlyinginNature,1975,2:729-762.[64]Weis-FoghT.EnergeticsofHoveringFlightinHummingbirdsandinDrosophila[J].JExpBiol,1972,56:79-104.[65]Weis-FoghT.BiologyandPhysicsofLocustFlight:VIILiftandMetabolicRateofFlyingInsects[J].JExpBiol,1964,41:257-271.[66]Weis-FoghT.UnusualMechanismsfortheGenerationofLiftinFlyingAnimals[J].ScientficAmerican,1975,233:80-87.153 北京理工大学博士学位论文[67]Weis-FoghT.QuickEstimatesofFlightFitnessinHoveringAnimals,includingNovelMechanismsforLiftProduction[J].JExpBiol,1973,59:169-230.[68]MaxworthyT.ExperimentsontheWeis-FoghMechanismofLiftGenerationbyInsectsinHoveringFlight.Part1Dynamicsofthefling[J].FluidMesh,1979,93:47-63.[69]LighthillMJ.OntheWeis-FoghMechanismofLiftGeneration[J].JFluidMech.1973,60:1-17.[70]EllingtonCP.TheAerodynamicsofHoveringInsectFlight.I.TheQuasi-SteadyAnalysis[J].PhilTransRSocLondB,1984,305:1-15.[71]EllingtonCP.TheAerodynamicsofHoveringInsectFlight.II.Morphologicalparameters[J].PhilTransRSocLondB,1984,305:17-40.[72]EllingtonCP.TheAerodynamicsofHoveringInsectFlight.III.Kinematics[J].PhilTransRSocLondB,1984,305:41-78.[73]EllingtonCP.TheAerodynamicsofHoveringInsectFlight.V.AVortexTheory[J].PhilTransRSocLondB,1984,305:115-144.[74]EllingtonCP.TheAerodynamicsofHoveringInsectFlight.VI.LiftandPowerRequirements[J].PhilTransRSocLondB,1984,305:145-181.[75]WillmottAP,EllingtonCP,ThomasALR.FlowVisualizationandUnsteadyAerodynamicsintheFlightoftheHawkmoth,ManducaSexta[J].PhilTransRSocLondB,1997,352:303-316.[76]GilmourKM,EllingtonCP.PowerOutputofGlycerinatedBumblebeeFlightMuscle[J].JExpBiol,1993,183:77-100.[77]WakelingJM.,EllingtonCP.Dragonglyflight.I.GlidingFlightandSteady-stateAerodynamicForces[J].JExpBiol,1997,200:543-556.[78]WakelingJM,EllingtonCP.Dragonglyflight.II.Velocities,AccelerationsandKinematicsofFlappingFlight[J].JExpBiol,1997,200:557-582.[79]WakelingJM,EllingtonCP.Dragonglyflight.III.LiftandPowerRequirements[J].JExpBiol,1997,200:583-600.[80]VanDenBergC,EllingtonCP.TheThree-dimensionalLeading-edgeVortexofa‘Hovering’ModelHawkmoth[J].PhilTransRSocLondB,1997,352:329-340.[81]VanDenBergC,EllingtonCP.TheVortexWakeofa‘Hovering’ModelHawkmoth[J].PhilTransRSocLondB,1997,352:317-328.154 北京理工大学博士学位论文[82]EllingtonCP.TheNovelAerodynamicsofInsectFlight:ApplicationstoMicro-airVehicles[J].JExpBiol,1999,202:3439-3448.[83]TytellED,EllingtonCP.HowtoPerformMeasurementsinaHoveringAnimal"sWake:PhysicalModellingoftheVortexWakeoftheHawkmoth,ManducaSexta[J].PhilTransRSocLondB,2003,358:1559-1566.[84]DickinsonMH,LehmannFO,SaneSP.WingRotationandtheAerodynamicBasisofInsectFlight[J].Science,1999,284:1954-1960.[85]SaneSP,DickinsonMH.TheAerodynamicEffectsofWingRotationandaRevisedQuasi-steadyModelofFlappingFlight[J].JExpBiol,2002,205:1087-1096.[86]SaneSP,DickinsonMH.TheControlofFlightForcebyaFlappingWing:LiftandDragProduction[J].JExpBiol,2001,204:2607-2626.[87]BirchJM,DicksonWB,DickinsonMH.ForceProductionandFlowStructureoftheLeadingEdgeVortexonFlappingWingsatHighandLowReynoldsNumbers[J].JExpBiol,2004,207:1063-1072.[88]BirchJM,DickinsonMH.TheInfluenceofWing-wakeInteractionsontheProductionofAerodynamicForcesinFlappingFlight[J].JExpBiol,2003,206:2257-2272.[89]FrySN,SayamanR,DickinsonMH.TheAerodynamicsofFree-flightManeuversinDrosophila[J].Science,2003,300:495-498.[90]FrySN,SayamanR,DickinsonMH.TheAerodynamicsofHoveringFlightinDrosophila[J].JExpBiol,2005,208:2303-2318.[91]LehmannFO,SaneSP,DickinsonM.TheAerodynamicEffectsofWing-wingInteractioninFlappingInsectWings[J].JExpBiol,2005,208:3075-3092.[92]LentinkD,DickinsonMH.RotationalAccelerationsStabilizeLeadingEdgeVorticesonRevolvingFlyWings[J].JExpBiol,2009,212:2705-2719.[93]LiuH,EllingtonCP,KawachiK.,VanDenBergC.AComputationalFluidDynamicStudyofHawkmothHovering[J].JExpBiol,1998,201:461-477.[94]LiuH.IntegratedModelingofInsectFlight:FromMorphology,KinematicstoAerodynamics[J].JComputPhy,2009,228:439-459.[95]AonoH,LiuH.Near-fieldandfar-fieldAerodynamicsinInsectHoveringFlight:anIntegratedComputationalStudy[J].JExpBiol,2008,211:239-257.155 北京理工大学博士学位论文[96]AonoH,LiuH.VorticalStructureandAerodynamicsofHawkmothHovering[J].JBiomechSciEng,2006,1:234-245.[97]NakataT,LiuH.AerodynamicPerformanceofaHoveringHawkmothwithFlexibleWings:aComputationalApproach[J].JExpBiol,2011,279:722-731.[98]BarkerDB,FourneyME.MeasuringFluidVelocitieswithSpecklePatterns[J].OpticsLetters,1977,1(4):135-137.[99]AdrianRJ.ScatteringParticleCharacteristicsandTheirEffectonPulsedLaserMeasurementsofFluidFlow:SpeckleVelocimetryvsParticleImageVelocimetry[J].AppliedOptics,1984,23(11):1690-1691.[100]PrasadAK.StereoscopicParticleImageVelocimetry[J].ExpinFluids,2000,29:103-116.[101]SpeddingGR,PennycuickCJ.MomentumandEnergyintheWakeofaPigeoninSlowFlight[J].JExpBiol,1984,111:81-102.[102]SpeddingGR,RaynerJMV,PennycuickCJ.QuantitativeStudiesoftheWakesofFreely-flyingBirdsinaLowTurbulenceWindTunnel[J].ExpFluids,2003,34:291-303.[103]SpeddingGR,HedenströmA.PIV-basedInvestigationsofAnimalFlight[J].ExpFluids,2009,46:749-763.[104]MurphyJ,HuH.AnExperimentalStudyofaBio-InspiredCorrugatedAirfoilforMicroAirVehicleApplications[J].ExpFluids,2008,49:531-546.[105]PoelmaC,DicksonWB,DickinsonMH.Time-ResolvedReconstructionoftheFullVelocityFieldAroundaDynamically-ScaledFlappingWing[J].ExpFluids,2006,41:213-225.[106]LuY,ShenGX,LaiGJ.DualLeading-edgeVorticesonFlappingWings[J].JExpBiol,2006,209:5005-5016.[107]LuY,ShenGX.Three-dimensionalFlowStructuresandEvolutionoftheLeading-edgeVorticesonaFlappingWing[J].JExpBiol,2008,211:1221-1230.[108]CombesSA,DanielTL.IntoThinAir:ContributionsofAerodynamicandInertial-elasticForcestoWingBendingintheHawkmothManduca[J].JExpBiol,2003,206:2999-3006.[109]CombesSA,DanielTL.FlexuralStiffnessinInsectWingsI.ScalingandtheInfluenceofWingVenation[J].JExpBiol,2003,206:2979-2987.[110]CombesSA,DanielTL.FlexuralStiffnessinInsectWingsII.SpatialDistributionandDynamic156 北京理工大学博士学位论文WingBending[J].JExpBiol,2003,206:2989-2997.[111]CombesSA.,CrallJD,MukherjeeS.DynamicsofAnimalMovementinanEcologicalContext:DragonflyWingDamageReducesFlightPerformanceandPredationSuccess[J].BiolLett,2010,6:426-429.[112]WoottonRJ,EvansKE,HerbertR,SmithCW.TheHindWingoftheDesertLocustI.FunctionalMorphologyandModeofOperation[J].JExpBiol,2000,203:2921-2931.[113]HuH,KumarAG,AlbertaniR.AnExperimentalInvestigationontheAerodynamicPerformancesofFlexibleMembraneWingsinFlappingFlight[J].AerospSciTechnol,2010,14(8):575-586.[114]MountcastleAM,DanielTL.AerodynamicandFunctionalConsequencesofWingCompliance[J].ExpFluids,2009,46(5):873-882.[115]DrzewieckiS.DesHelicesAeriennesTheorieGeneraledesPropulseurs[M].LibofAero,1909.[116]OsborneMFM.AerodynamicsofFlappingFlightWithApplicationtoInsects[J].JExpBiol,1951,28(2):221-245.[117]WhitneyJP,WoodR.J.AeromechanicsofPassiveRotationinFlappingFlight[J].JFluidMech,2010,660:197-220.[118]HubelTY,TropeaC.ExperimentalInvestigationofaFlappingWingModel[J].ExpFluids,2009,46:945-961.[119]DudleyR,EllingtonCP.MechanicsofForwardFlightinBumblebees:II.Quasi-steadyLiftandPowerRequirements[J].JExpBiol,1990,148:53-88.[120]DudleyR.TheBiomechanicsofInsectFlight:Form,Function,Evolution[M].PrincetonUniversityPress,2000.[121]GrodnitskyDL.FormandFunctionofInsectWings:TheEvolutionofBiologicalStructures[M].TheJohnHopkinsUniversityPress,1999.[122]WagnerH.UberdieEntstehungdesDynamischenAuftriebesvonTragflUgeln[J].ZeitschriftfurAngewandteMathematikuntMechanik,1925,5(17):17-35.[123]陈玉璞,王惠民.流体动力学[M].清华大学出版社,2003[124]WagnerH.UberdieEntstehungdesDynamischenAuftriebesvonTragflUgeln[J].ZeitschriftfurAngewandteMathematikuntMechanik,1925,5(17):17-35.[125]UsherwoodJR,EllingtonCP.TheAerodynamicsofRevolvingWingsI.ModelHawkmothWings157 北京理工大学博士学位论文Ellington[J].JExpBiol,2002,205:1547-1564.[126]WuJH,SunM.UnsteadyAerodynamicForcesofaFlappingWing[J].JExpBiol,2004,207(7):1137-1150.[127]HongY,AltmanA.StreamwiseVorticityinSimpleMechanicalFlappingWings[J].JournalofAircraft,2007,44(5):1588-1597.[128]HongY,AltmanA.LiftFromSpanwiseFlowinSimpleFlappingWings[J].JournalofAircraft,2008,45(4):1206-1216.[129]LuttgesM.AccomplishedInsectFliers[J].FrontiersinExperimentalFluidMechanics,1989,46:429-456.[130]KlissM,SompsC,LuttgesMW.StableVortexStructures:AFlatPlateModelofDragonflyHovering[J].TheorBiol,1989,136:209-228.[131]SompsC,LuttgesM.Dragonflyflight:Novelusesofunsteadyseparatedflows[J].Science,1985,228(4705):1326-1329.[132]SrygleyRB,ThomasALR.UnconventionalLift-generatingMechanismsinFree-flyingButterflies[J].Nature,2002,420:660-664.[133]TaylorGK,NuddsRL,ThomasALR.FlyingandSwimmingAnimalsCruisataStrouhalNumberTunedforHighPowerEfficiency[J].Nature,2003,425:707-711.[134]TuncerIH,PlatzerMF.ThrustGenerationDuetoAirfoilFlapping[J].AIAAJ,1996,34:324-331.[135]EllingtonCP,VanDenBergC,WillmottAP,ThomasALR.LeadingEdgeVorticesinInsectFlight[J].Nature,1996,384:626-630.[136]LentinkD,DickinsonMH.BiofluiddynamicScalingofFlapping,WpinningandTranslatingFinsandWings[J].JExpBiol,2009,212:2691-2704.[137]SunM,WuJH.LargeAerodynamicForcesonaSweepingWingatLowReynoldsNumber[J].ActaMechanicaSinica,2004,20(1):24-31.[138]BirchJM,DickinsonMH.SpanwiseFlowandtheAttachmentoftheLeading-edgeVortexonInsectWings[J].Nature,2001,412:729-733.[139]ShyyW,TrizillaP,KangCk,AonoH.CanTipVorticesEnhanceLiftofaFlappingWing[J].AIAAJ,2009,47(2):289-293.[140]DeVoriaAC,RinguetteMJ.VortexFormationandSaturationforLow-aspect-ratioRotating158 北京理工大学博士学位论文Flat-platefins[J].ExpFluids,2012,52(2):441-462.[141]CarrZ,ChenC,RinguetteMJ.TheEffectofAspectRatioontheThree-dimensionalVortexFormationofRotatingFlat-plateWings[C].In50thAIAAAerospaceSciencesMeetingincludingtheNewHorizonsForumandAerospaceExposition,2012.[142]DeVoriaAC,RinguetteMJ.VortexFormationandSaturationforLow-aspect-ratiorotatingFlatPlatesatLowReynoldsNumber[J].ExpFluids,2012,52:441-462.[143]AnsariSA,PhillipsN,StablerG..ExperimentalInvestigationofSomeAspectsofInsect-likeFlappingFlightAerodynamicsforApplicationtoMicroAirVehicles[J].ExpFluids,2009,46:777-798.[144]WilkinsP,KnowlesK.InvestigationofAerodynamicsRelevanttoFlappingwingMicroAirVehicles[C].In:37thAIAAFluidDynamicsConferenceandExhibition,2007,[145]WilkinsP,KnowlesK.TheLeading-edgeVortexandAerodynamicsofInsect-basedFlapping-wingMicroAirVehicles[J].AeronautJ,2009,113:253-262.[146]GopalakrishnanP,TaftiDK.EffectofRotationKinematicsandAngleofAttackonFlappingFlight[J].AIAAJ,2009,47(11):2505-2519.[147]AlexanderDE.Nature’sFlyers:Birds,Insects,andtheBiomechanicsofFlight[M].TheJohnsHopkinsUniversityPress,2002.[148]BrodskyK.VortexFormationintheTetheredFlightofthePeacockButterflyInachisioL.(Lepidoptera,Nymphalidae)andsomesepectsofinsectflightevolution[J].JExpBiol,1991,161:77-95.[149]WarrickDR,TobalskeBW,PowersDR.AerodynamicsoftheHoveringHummingbird[J].Nature,2005,435(23):17-130.[150]SpeddingGR.TheWakeofaJackdaw(CorvusMonedula)inSlowFlight[J].JExpBiol,1986,125:287-307.[151]SpeddingGR.TheWakeofaKestrel(FalcoTinnunculus)inGlidingFlight[J].JExpBiol,1987,127:45-57.[152]SpeddingGR,RosenM,HedenstromA.AFamilyofVortexWakesGeneratedbyaThrushNightingaleinFreeFlightinaWindTunneloveritsEntireNaturalRangeofFlightSpeeds[J].JExpBiol,2003,206:2313-2344.[153]GrodnitskyDL,MorozovPP.FlowVisualizationExperimentsonTetheredFlyingGreenLacewingsChrysopaDasyptera[J].JExpBiol,1992,169:143-163.159 北京理工大学博士学位论文[154]AltshulerDL,PrincevacM,PanH,LozanoJ.WakePatternsoftheWingsandTailofHoveringHummingbirds[J].ExpFluids,2009,46:835-846.[155]WolfM,Ortega-JimenezVM,DudleyR.StructureoftheVortexWakeinHoveringAnna’sHummingbirds(Calypteanna)[J].ProcBiolSci,2013,280:1-7.[156]SaneSP.TheAerodynamicsofInsectFlight[J].JExpBiol,2003,206:4191-4208.[157]DickinsonMH.TheEffectsofWingRotationonUnsteadySerodynamicPerformanceatLowReynoldsNumbers[J].JExpBiol,1994,192:179-206.[158]WuJH.TheInfluenceoftheWakeofaFlappingWingontheProductionofAerodynamicForces[J].ActaMechanicaSinica,2005,21(5):411-418.[159]LentinkD,VanHeijstGF,MuijresFT,VanLeeuwenJL.VortexInteractionswithFlappingWingsandFinscanbeUnpredictable[J].BiolLett,2010,6(3):394-397.[160]LehmannFO.WhenWingsTouchWakes:UnderstandingLocomotorForceControlbyWake-wingInterferenceinInsectWings[J].JExpBiol,2008,211:224-233.[161]章梓雄,董曾南.粘性流体力学[M].清华大学出版社,2011.[162]WangZJ.VortexSheddingandFrequencySelectioninFlappingFlight[J].JFluidMech,2000,410:323-341.[163]http://wenku.baidu.com/view/29cd3d0a581b6bd97f19ea2e.html[164]童秉刚,尹协远.涡运动理论[M].中国科学技术大学出版社,2009.[165]BatchelorGK.AnIntroductiontoFluidDynamics[M].CambridgeUniversityPress,1967.[166]尹协远,孙德军.旋涡流动的稳定性[M].国防工业出版社,2003.[167]丁祖荣.流体力学(上册)[M].高等教育出版社,2003.[168]吴子牛.空气动力学(上册)[M].清华大学出版社,2007.[169]张鹏飞.扑翼飞行的理论与实验研究[D].北京理工大学硕士学位论文,2010.[170]刘岚.微型扑翼飞行器的仿生翼设计技术研究[D].西北工业大学博士学位论文,2007.[171]DudleyR,EllingtonCP.MechanicsofForwardFlightinBumblebees:I.KinematicsandMorphology[J].JExpBiol,1990,148:19-52.[172]李周复.风洞特种实验技术[M].航空工业出版社,2010[173]WagterCD,KoopmansA,CroonGD.AutonomousWindTunnelFree-FlightofaFlappingWingMAV[J].AdvancesinAerospaceGuidance,NavigationandControl,2013,45:603-621.160 北京理工大学博士学位论文[174]ZakariaMY,BayoumyAM,ElshabkaAM.ExperimentalAerodynamicCharacteristicsofFlappingMembraneWings[C].13thInternationalConferenceonaerospacescience,2009:1-18.[175]KimDK,HanJH,KwonKJ.WindTunnelTestsforaFlappingWingModelwithaChangeableCamberusingMacro-fiberCompositeActuators[J].SmartMaterStruct,2009,18:113-121.[176]LiuL,ZhangXJ,HeZX.AerodynamicAnalysisandWindTunnelTestforFlappingMavs[J].JournalofTheoreticalandAppliedInformationTechnology,2012,45(2):542-550.[177]PennycuickCJ.Wind-tunnelStudyofGlidingFlightinthePigeonColumbaLivia[J].JExpBiol,1968,49:509-526.[178]JonesMB,ValiyffA,HarveyJ.FlexibilityofanOrnithopterWingTestedinaWindTunnel[C].28thInternationalCongressoftheAeronauticalSciences,2012.[179]ChenP,JoshiS,SwartzSM,BreuerKS,BahlmanJW.BatInspiredFlappingFlight[C].22ndAIAA/ASME/AHSAdaptiveStructuresConference,2014.[180]HickleC.WindTunnelRenovation,FlowVerificationandFlappingWingAnalysis[M].AmazonDigitalServices,2012.[181]ThomasALR,TaylorGK,SrygleyRB,NuddsRL,BomphreyRJ.Dragonflyflight:Free-flightandTetheredFlowVisualizationsRevealaDiverseArrayofUnsteadyLift-generatingMechanisms,ControlledPrimarilyviaAngleofAttack[J].JExpBiol,2004,207:4299-4323.[182]BomphreyRJ,TaylorGK,LawsonNJ,ThomasALR.SmokeVisualizationofFree-flyingBumblebeesIndicatesIndependentLeading-edgeVorticesoneachWingPair[J].ExpFluids,2009,46:811-821.[183]刘政崇.风洞结构设计[M].中国宇航出版社,2005.[184]任荣林,王振羽.风洞设计原理[M].北京航空学院出版社,1985.[185]中国人民解放军总装备部军事训练教材编辑工作委员会编.低速风洞实验[M].国防工业出版社,2002.[186]PopeA,HarperJJ,WileyJ.彭锡铭等译.低速风洞实验[M].国防工业出版社,1977.[187]陈鑫.矩形叶栅风洞设计及流场品质分析[D].大连海事大学硕士论文,2011.[188]ReynoldsO.AnExperimentalInvestigationofCircumstanceswhichDetermineWhethertheMotionofWaterShallbeDirectorSinuousandoftheLawofRegistanceinParallelChannels[J].PhilosophicalTransactionsoftheRoyalSocietyofLondon,1883,174:935-982.161 北京理工大学博士学位论文[189]罗明晖,呼和敖德.国外流动显示技术进展[J].力学进展,1983,13(4):402-415.[190]束继祖,李华煜.流场显示技术在流体力学中的应用和展望[J].力学进展,1979,9(1):26-37.[191]范洁川.流动显示技术的若干现状与发展[J].气动实验与测量控制,1995,9(1):10-17.[192]EvansRA.NewMethodofFlowVisualizationforLow-DensityWindTunnels[J].JournalofAppliedPhysics,1957,28(9):130-145.[193]崔尔杰,洪金森.流动显示技术及其在流体力学研究中的应用[J].空气动力学学报,1991,9(2):190-199.[194]王振东.踪迹随风叶,程途犯斗槎——漫谈流动显示及应用[J].力学与实践,1996,18(2):75-77.[195]林建忠等.流体力学[M].清华大学出版社,2005.[196]范洁川.近代流动显示技术[M].国防工业出版社,2002.[197]KitagawaK,SakakibaraM,YasuharaM.VisualizationofFlappingWingoftheDroneBeetle[J].JournalofVisualization,2009,4:393-400.[198]ThomasE.Forloveofinsects[M].HarvardUniversityPress,2005.[199]AltshulerDL,DicksonWB,VanceJT.Short-amplitudeHigh-frequencyWingStrokesDeterminetheAerodynamicsofHoneybeeFlight[J].ProcNatlAcadSciUSA.2005,102(50):18213-18218.[200]ZhangGJ,SunJH,ChenDZ.FlappingMotionMeasurementofHoneybeeBilateralWingsusingFourVirtualStructured-light[J].SensorsandActuatorsA:Physical,2008,148(11):19-27.[201]TruongTV,LeTQ,TranHT,ParkHC.FlowVisualizationofRhinocerosBeetle(Trypoxylusdichotomus)inFreeFlight[J].JournalofBionicEngineering,2012,9:304-314.[202]VidelerJJ,StamhuisEJ.Leading-edgevortexliftsswifts.Science,2004,306:1960–1962.[203]FordCWP,BabinskyH.LiftandtheLeading-edgeVortex[J].JFluidMech,2013,720:280-313.[204]MuijresFT,JohanssonLC,BarfieldR,WolfM,SpeddingGR,HedenströmA.Leading-EdgeVortexImprovesLiftinSlow-FlyingBats[J].Science,2008,319:1250-1253.[205]WangZJ.DissectingInsectFlight[J].AnnRevFluidMech,2005,37:183-210.[206]IimaM.ATwo-dimensionalAerodynamicModelofFreelyFlyingInsects[J].JournalofTheoreticalBiology,2007,247:657-671.[207]MinottiFO.UnsteadyTwo-dimensionalTheoryofaFlappingWing[J].PhysicalreviewE,2005,66:1907-2002.162 北京理工大学博士学位论文[208]KeigoO,KosukeS,TakajiI.LiftGenerationbyaTwo-dimensionalSymmetricFlappingWing:ImmersedBoundary-latticeBoltzmannSimulations[J].FluidDynRes,2012,44:33-50.[209]LeeKB,KimJH,KimC.AerodynamicEffectsofStructuralFlexibilityinTwo-DimensionalInsectFlappingFlight[J].Jounalofaircraft,2011,48(3):894-909.[210]MillerLA,PeskinCS.AComputationalFluidDynamicsofClapandFling’intheSmallestInsects[J].JExpBiol,2005,208:195-212.[211]RamamurtiR,SandbergWC.AComputationalInvestigationoftheThree-dimensionalUnsteadyAerodynamicsofDrosophilaHoveringandManeuvering[J].TheJournalofExperimentalBiology2007,210:881-896.[212]SunM,WangJ,Xiong,Y.DynamicFlightStabilityofHoveringInsects[J].ActaMechanicaSinica,2007,23:231-246.[213]LiuH,KawachiK.AnumericalStudyofInsectflight[J].JCompPhys,1998,146:124-156.[214]YuX,SunM.AComputationalStudyoftheWing-wingandWing-bodyInteractionsofaModelInsect[J].ActaMechSin,2009,25:421-431.[215]韩占忠,王敬,兰小平.FLUENT-流体力学仿真计算实例与应用[M].北京理工大学出版社,2004.[216]张涵信,邓小刚.内点、边界、网格与物理分析相结合的非定常模拟[M].国防工业出版社,2006.[217]王福军.计算流体动力学分析:CFD软件原理与应用[M].清华大学出版社,2004.[218]管德.非定常空气动力计算[M].北京航空航天大学出版社,1991.[219]于勇.FLUENT入门与进阶教程[M].北京理工大学出版社,2008.[220]王瑞金.Fluent技术基础与实例分析[M].清华大学出版社,2007.[221]王惠民.流体力学基础[M].清华大学出版社,2005.[222]RamamurtiR,SandbergWC.AThree-dimensionalComputationalStudyoftheAerodynamicMechanismsofInsectFlight[J].JExpBiol,2002,205:1507-1518.[223]PerssonPO,WillisDJ,PeraireJ.TheNumericalSimulationofFlappingWingsatLowReynoldsNumbers[C].48thAIAAAerospaceSciencesMeetingIncludingtheNewHorizonsForumandAerospaceExposition,2010.[224]YuML.NumericalandExperimentalInvestigationsonUnsteadyAerodynamicsofFlapping163 北京理工大学博士学位论文Wings[D].IowaStateUniversity,2012.164 北京理工大学博士学位论文攻读学位期间发表论文与研究成果清单[1]第一作者.蜜蜂振翅的流场显示实验和理论分析[J].中国科学E辑:技术科学,2014,44(11):1211-1221.[2]第一作者.扑翼模型的气动力研究[J].北京理工大学学报,2014,34:172-175.(EI检索).[3]第二作者.扑翼模型运动方程的优化[J].北京理工大学学报,2011,31:178-182.(EI检索).[4]第一作者.TheDesignofaLow-speedWindTunnelforStudyingoftheFlowFieldofInsects’Flight[C].InternationalConferenceonExperimentalMechanics(ICEM2014)inconjunctionwiththe13thAsianConferenceonExperimentalMechanics,15-17November,2014,Singapore,Singapore.(已录用,EI检索).[5]第一作者.ExperimentalResearchonFlightPatternsofBees[C].2014InternationalConferenceonAppliedMechanics,Modelling,andSimulation,16-17October,2014,Xiamen,China.(已录用,EI检索).165 北京理工大学博士学位论文致谢转眼之间，六年多的求学时光悄然而逝。在这人生中最重要的求学历程中，首先要感谢导师马巍研究员，在攻读博士研究生期间，深深受益于马老师的关心、爱护和谆谆教导。马老师严谨的治学态度，兢兢业业雷厉风行的工作作风以及一丝不苟的科研精神成为指导我奋力前行的指明灯塔。能师从马老师，我为自己感到庆幸，在此谨向马老师表示我最诚挚的敬意和感谢！需要特别感谢宁建国教授，在仿生扑翼力学领域，宁老师给予了耐心的指导和极大的帮助，谨向宁老师致以我最衷心的感谢和最诚挚的祝福！攻读博士学位期间，还要感谢任会兰老师、马天宝老师、宋卫东老师、王成老师。无论是在学习上还是在生活中，各位老师的耐心指导和悉心教诲给予我无穷的动力，伴我一路前行，在此向你们表示最诚挚的谢意！感谢张鹏飞、马赟、储著鑫，在求学的道路上，大家攻坚克难，一往直前，前行的道路上洒满了我们拼搏的汗水。感谢郭婷婷，费广磊、林庚浩、方敏杰、贺月香、李健、于超、耿荻、栗建桥、孙春勇、赵海涛，感谢王婧、段妍、王星、展婷变、李艳梅、赵宇哲、刘晓俊、甘元超、褚亮、袁新鹏等人，能有幸在如此优秀的团队中，同大家一起学习工作，让我受益颇多。过往岁月的点点滴滴都让我如此难忘，在此向大家表示感谢。特别要感谢我的父母，在我求学的道路上，他们自始至终给予了无限的理解、支持和关爱，他们是我最坚实的后盾。他们的鼓励和帮助使我得以顺利完成博士的学习和研究。166