微型扑翼飞行器动力学模型参数的灵敏度分析

毕富国，何广平

(北方工业大学 机械与材料工程学院，北京 100144)

在1992年，美国国防高级研究计划局(Defense Advanced Research Projects Agency,DARPA)首次提出了微型飞行器(MAV)的概念，并认为其对未来战争格局有深远影响。在DAPRA计划的影响下，各类微型飞行器在军事中得到应用。而微型扑翼飞行器(FWMAV)相较于其他无人机[1]，因为其具有的高机动性、可悬停和隐蔽性好等优点而备受关注。随着微制造技术[2]的发展和昆虫飞行动力学[3-4]的深入研究，目前开展微型扑翼飞行器的研究已经成为机器人学领域的一个热点[5-7]。与传统固定翼和旋翼飞行器相比。微型扑翼飞行器所具有的高气动效率，抗干扰能力以及可悬停等特点依赖于非定常空气动力学现象，如延迟失速、旋转升力和尾迹捕捉[8]等。但由于尺度微小，微型扑翼飞行器的负载能力十分有限，系统设计中所需的传感器、电源、通信设备、控制设备等机载器件对飞行器的负载能力提出了较高要求。在研制微型扑翼飞行器时，其涉及的低雷诺数非定常空气动力学特点和性质是必须考虑的。

在小型扑翼飞行器研究方面，代表性地，如荷兰代尔夫特理工大学的DelFly系列[9]，美国的Nano Hummingbird[10]等。目前，有关微型扑翼飞行器的研究主要以哈佛大学的RoboBee系列[11]为主，该微型扑翼飞行器由初期的80毫克升力逐步提升到253mg升力。在以往的研究中，RoboBee研究团队主要通过实验来定量研究机翼形态和惯性参数对扑翼飞行性能的影响。通过在不同机翼形状和运动参数[12]下进行大量无源动态运动实验，比较气动力特性，逐步修正和优化飞行器的气动性能。本研究基于叶素理论建立气动力模型，引入Sobol全局灵敏度分析方法，对建立的微型扑翼飞行器气动模型中的设计参数进行定量分析，得到参数对升力的影响系数，利用数值方法研究参数与升力的关系，可以更有针对性地对主要的设计参数进行优化，为微型扑翼飞行器的设计和控制提供参考。

1 气动力建模

Sane[13]研究表明，基于准稳态空气动力学模型，根据实验结果对气动力系数进行修改，可以合理地预测由于延迟失速引起的瞬时气动力，而无需使用调节扑翼飞行空气动力学偏微分方程的复杂有限元解。因此，下面给出的空气动力学模型是基于准稳态方程的。准稳态方程不能模拟尾迹捕获这种复杂的非稳态空气动力学现象，但是已有研究表明，在类正弦运动中，尾迹捕捉所产生的气动力对整体气动力贡献较小，准稳态理论对昆虫尺度的微型扑翼飞行器的气动力分析是适用的。

在气动力计算中采用叶素法来对气动力进行计算，可以给出每个微元的气动力计算公式[12]：

(1)

其中,U为机翼相对速度，如图1所示；α为机翼与拍动平面的夹角；Caero(α)为关于攻角的气动力系数。下面确定公式中的参数，从而能计算每个微元上的气动力，如图2所示。

图2 参数化机翼示意图

图2中，R为翼根到机翼最远点的长度；r为相对于翼根的距离；c(r)为距离翼根r处的展长分布；dr为沿展向的微元。

1.1 气动力系数的确定

Dickinson[14]通过实验研究总结出了关于昆虫气动力系数的经验公式，并给出其应用范围100

哈佛大学在研究过程中所采用的气动力升力系数和阻力系数计算公式如下：

(2)

其中，CLmax=1.8，CD0=0.4，CDmax=3.4；将各参数代入式(2)可得：

CL(α)=1.8sin(2α)

CD(α)=1.9-1.5cos(2α)

并通过几何关系得：

CN(α)=cos(α)CL+sin(α)CD

(3)

再将各个系数代入微元方程(1)中，积分后获得单个机翅的升力、阻力及法向力。

伯克利的Deng等[15]基于实验数据给出了如下的法向力和切向力系数，以表达考虑了延迟失速现象的气动力系数：

(4)

并通过力平衡分析得到升力系数和阻力系数与法向力系数个切向力系数之间的关系：

CL(α)=cos(α)CN(α)-sin(α)CT(α)

CD(α)=sin(α)CN(α)+cos(α)CT(α)

(5)

进一步将各个系数代入微元方程(1)中，积分后获得单个翅膀的法向力、切向力以及对应的升力和阻力。

本文通过给定攻角函数并代入两种方法来对其准确性进行验证，获得结果如图3所示。

图3 攻角与气动力的关系

可以看出，以上通过经验公式得到的各组气动力系数，当攻角在约45°时，气动升力系数达到最大值。两者之间的差值则主要由于实验方法的不同导致的。

由于哈佛大学的Robobee系列已经展示了较好的飞行特性，因此，在下面均采用哈佛大学所给出的气动力计算公式，且攻角公式由文献[16]给出：

α(t)=αmaxsin(2πft+δ)

(6)

假设起始位置的平均攻角面偏差δ为0，αmax=45°。

由Dickinson[3]试验获得的平均升力系数表明，在发生超前和对称运动时昆虫获得的升力较大，但在飞行器设计过程中，要实现超前运动存在很大困难，因此在本文中不考虑超前滞后现象，采用名义方程：

α(t)=αmaxsin(2πft+δ)

(7)

将式(7)代入气动力系数计算公式(2)，取机翼向前扑动为正方向，假设机翼在行程变换时的旋转时瞬间完成，得到图4所示气动力系数计算结果。

从计算结果可以看出，升、阻力系数在运动前半程和后半程先增大后减小并中间时刻出现最大值。由于在前半程和后半程运动变换时，机翅的方向发生改变。升力为机翅的切线方向，升力系数为0，阻力为机翅的法线方向，阻力系数不为0。

1.2 延展向方向弦长分布的确定

在使用叶素法计算升力时，需要求解沿机翼翼展向方向的弦长分布。

图4 气动力系数计算结果

(8)

一阶半径矩是翅膀的中心区域，取值在0.4～0.6[12]，二阶半径矩描述翅膀的面积分布。

并发现翅膀形状可以用β分布来表示：

(9)

其中B(p,q)为β方程：

(10)

其中：

(11)

将无量纲方程进行推导，得到微型扑翼飞行器的翅翼轮廓方程：

(12)

1.3 机翼相对速度的确定

悬停状态下，来流速度为0。因此，相对速度取机翅压力中心速度[15]：

(13)

其中，R为机翼长度，ω为机翼拍动角速度。

将气动力系数、弦长分布以及相对速度参数方程代入微元气动力计算式(1)，得到

(14)

(15)

其中空气密度ρ=1.205 kg/m3。

沿展向积分得到气动力关于攻角的计算方程：

(16)

(17)

则单个翅膀的气动力为升力和阻力的矢量和：

F=FL+FD

(18)

2 Sobol全局灵敏度分析

上一节所述气动力模型呈高度非线性特征，且各参数之间存在耦合。传统的局部灵敏度分析方法，如直接求导法、有限差分法和摄动法等是基于线性或非线性模型，对于对非线性较强的模型，局部灵敏度分析法无法提供合理的分析结果。全局灵敏度分析方法则适用于高非线性模型、非单调模型、可分析全范围参数，并可以考察参数间相互作用，避免了局部分析方法的局限性。因此，在气动力参数灵敏度分析过程中考虑其特点采用全局灵敏度分析方法，具有重要的实际应用价值。其中，Sobol方法是研究全局灵敏度常用的一种方法[18～20]。

2.1 Sobol分析方法理论

Sobol方法是一种基于方差分解的全局灵敏度定量分析方法。气动力模型Y=f(X)，X=x1，x2，x3,…,xk,其中x1，x2，x3,…,xk为模型中待分析的输入参数。Y为模型输出值。可以将模型分解为如下形式：

f1,2,…,k(x1,x2,…,xk))

(19)

其中：

f0=E(Y)

fi(xi)=E(Y|xi)-f0

fij(xi，xj)=E(Y|xi，xj)-f0-fi-fj

从其中可以看出，fi是单独改变xi的效果(称为xi的主要效应)，并且fij表示除了它们各自变化的效果之外，同时改变xi和xj对输出的效果，这被称为二阶交互。高阶项具有类似的定义。

如果输入参数之间是相互独立的，则方程式(19)右侧各项相互正交，协方差为0。对两侧求方差得：

(20)

其中：

Vi=Varxi(Ex～i(Y|xi))

Vij=Varx～ij(Ex～ij(Y|xi,xj))-Vi-Vj

x～i表示除了xi之外的所有变量的集合。

通过计算比率：

(21)

式中，Si是xi的一阶灵敏度系数，表示参数xi对模型输出的影响。

同理定义比率：

(22)

式中，STi是xi的全局灵敏度系数，表示考虑参数耦合作用后，xi对模型输出的影响。

2.2 Sobol法计算过程

首先，进行因子筛选。根据定义设定展长为15 mm，对气动力模型进行数值代入计算得到运动参数对气动力的影响，如图5所示。

图5 运动参数对气动力的影响

由图5可以看出，在其他参数确定的情况下，气动力随最大扑动角的增大而增大；随攻角的增大先增大后减小；随扑动频率增大而增大。在上述分析中只是分别对每个参数对气动力的影响进行了控制变量分析，并未考虑在各个设计参数共同影响下气动力的变化。在确定分析参数之后，对气动力模型分析参数进行两次采样，得到AN·D和BN·D矩阵，其中N为抽样参数，D为分析参数。另构造ABi(i=1,2，…，D)即用B的第i列替换A的第i列。

(23)

V(Y)=V(A)+V(B)

(24)

其中，根据文献[21]图6及表2给出的数据表明，在参数量为5的模型全局灵敏度分析中，取样数在500左右即可达到收敛。在本次全局灵敏度分析过程中取样数为1 000，满足取样要求。参数取值范围如表1所示。

表1 各参数取值范围

3 计算结果及分析

考虑到飞行器设计过程中机翼几何参数的加工难度和加工条件要求较高，且在飞行器实际飞行过程中主要依赖于运动参数进行控制，因此在以往研究过程中关于参数灵敏度分析的研究受到更多关注。本文同时进行了几何参数的灵敏度分析，从而为机翼设计提供参考。

通过Matlab计算得到气动力模型中机翼几何参数一阶半径矩、展弦比以及运动参数最大扑动角、最大攻角、扑动频率的全参数一阶灵敏度和全局灵敏度，计算结果分别是：总方差为4.6 e-2；一阶半径距方差为2.42 e-4；展弦比方差为1.19e-2；最大扑动角方差为1.579e-3； 最大攻角方差为1.576e-3；扑动频率方差为1.7e-2。

计算可得：一阶半径矩的一阶灵敏度为5.26e-3，全局灵敏度为2.43e-2；展弦比的一阶灵敏度为2.59e-1，全局灵敏度为5.28e-1；最大扑动角的一阶灵敏度为3.42e-2，全局灵敏度为1.41e-1；最大攻角的一阶灵敏度为3.43e-2，全局灵敏度为1.08e-1；扑动频率的一阶灵敏度为3.68e-1，全局灵敏度为5.88e-1。分析结果如图6所示。

图6 全参数灵敏度分析结果

进一步，按照上述分析方法分别对机翼设计中设计的机翼几何参数和运动控制参数进行灵敏度分析。

在独立的几何参数分析中设定最大扑动角为120°；最大攻角为45°；扑动频率为120。得到分析结果为一阶半径矩的一阶灵敏度为3.65e-2，全局灵敏度为4.48e-2；展弦比的一阶灵敏度为9.95e-1，全局灵敏度为1.01。结果如图7所示。

图7 几何参数灵敏度分析结果

在独立的运动参数分析中设定一阶半径矩为0.56；展弦比为7。得到分析结果为最大扑动角的一阶灵敏度为1.26e-1，全局灵敏度为2.32e-1；最大攻角的一阶灵敏度为6.95e-2，全局灵敏度为1.09e-1；扑动频率的一阶灵敏度为5.82e-1，全局灵敏度为7.46e-1。结果如图8所示。

图8 运动参数灵敏度分析结果

通过对比可以看出，全参数分析结果与对几何参数和运动参数分别分析得出的结果保持一致，从而证明了分析结果的有效性。全参数分析显示，对于气动力影响的因素大小依次为：扑动频率>展弦比>最大扑动角>最大攻角>一阶半径矩。运动参数全局灵敏度为扑动频率>最大扑动角>最大攻角。

相比于全局灵敏度，一阶灵敏度的排序与全局灵敏度基本一致，同时也可以看出参数之间存在影响。在进行参数优化设计过程中，全局灵敏度较小的参数可以选取固定值，这样既可以减小在优化时的计算成本，也为控制律设计提供了方便。

4 结论

建立了基于叶素法的微型扑翼飞行器的气动力学模型，并采用Sobol全局灵敏度分析方法对气动力模型进行灵敏度分析，获得了机翼几何参数和运动参数对气动力的影响灵敏度。分析表明，几何参数中展弦比的全局灵敏度最大，为5.28e-1，一阶半径距全局灵敏度最小，为2.44e-2；运动参数中扑动频率的全局灵敏度最大，为5.88e-2，扑动角的全局灵敏度为1.41e-2，攻角的全局灵敏度为1.08e-1。可以看出，要提高微型扑翼飞行器的空气动力性能，在结构参数方面，翅翼的展弦比需要进行精心设计，在运动参数方面，扑翼飞行器的攻角对扑翼飞行器的升力影响最大，在控制器设计中应优先进行稳定控制。