飞艇动导数与附加质量相互融合的方法

林献武，王仕超，李智斌，兰维瑶

1. 厦门大学 航空航天学院，厦门 361005

2. 山东科技大学 电气与自动化工程学院，青岛 266510

飞行力学模型是研究飞艇技术的一个重要基础，建立飞行力学模型的难点在于气动力和力矩的确定和处理。在长达一百多年的飞艇技术发展过程中，许多学者对飞艇气动力和力矩的计算、测量和处理方法进行了研究。19世纪及以前基于无旋无环流假设估算运动体的气动力和力矩，其理论由Lamb[1]总结；这个研究主要关注非定常气动力和力矩以及稳态的Munk力矩。飞艇长时间匀速直线平动所对应的气动力和力矩估算主要是基于20世纪20年代Munk提出的划片法[2]，于20世纪中叶由Allen和Perkins[3-4]推广到有黏流的情况并由Hopkins[5]进行改进；到了20世纪80年代，随着计算技术的发展面元法也被应用于飞艇定常气动力和力矩的数值计算中[6]。有黏流情况下，飞艇长时间匀速直线平动气动力和力矩实验测量的代表性工作有文献[7-10]等。有黏流中非定常气动力和力矩实验测量的代表性工作有文献[9-12]等。有黏流中的非定常气动力计算，早期的代表性工作有文献[12-13]，近年来随着CFD技术的发展，越来越多的学者开展基于CFD计算数据来辨识动导数的研究[14-17]。关于有黏流和无黏流对应气动系数的融合，代表性的工作有文献[13，18-20]等。

在刚体假设下，运动体上任一点的速度可以用运动体的角速度和其上一点(不妨称为平动参考点)的平动速度来表征。因此飞行力学模型中一般取这些当前运动参数为未知状态变量，并建立与未知数个数相对应的微分方程。为了在飞行力学模型中计入气动力或力矩的影响，一般需要将其与这些当前运动参数关联，这样才能使方程组封闭；这种表征方法是通过气动系数的概念来实现的。对于无旋无环流场而言，Lamb的理论表明运动刚体的气动力和力矩均取决于其当前运动参数，因此非定常气动系数的概念是存在的，一般称为附加质量。对于有黏或有旋流而言，运动刚体的非定常运动气动力或力矩不仅与其当前运动参数有关，同时也受其历史运动信息影响；因此传统意义上的气动系数并不是严格存在的。基于这个原因，现在飞艇飞行力学模型中的相当部分非定常气动系数仍采用无黏流中的对应结果。与俯仰、偏航、升沉或侧滑运动有关的气动力或力矩，用于确定飞艇的扰动特性；这时其所对应的运动接近周期性振荡过程，因此这些气动力和力矩仍可近似与其当前运动参数相关，如文献[14]中介绍的傅里叶分析法所示。正是基于这个原因，与飞艇俯仰、偏航、升沉或侧滑运动有关的气动系数存在有黏流中的结果，它们与对应当前运动参数的关联系数一般称为动导数。这时就存在一个问题：在飞艇飞行力学建模中，动导数与无旋无环流中对应的附加质量均表示非定常气动力或力矩与当前运动参数之间的关系，它们应该如何融合，现有代表性工作[13，18-20]中所采用的方法并不一致。例如文献[13]在气动系数融合时，同时采用了有黏流中的俯仰阻尼力/力矩和无旋无环流中对应的非定常气动力/力矩；在文献[19]中，既不采用有黏流中的动导数，也不采用无旋无环流中的对应结果；文献[18]在采用部分有黏流中的动导数后，摒弃了作者认为与之对应的无旋无环流中的结果；在文献[20]中，非定常气动系数仅采用无旋无环流中的结果，这种融合方法也被近期的许多专著所采用[21-23]。这说明人们在这个问题的认识上不但存在争议，而且多数文献回避了这个问题，因此需要进一步探讨。由于潜航器或传统重于飞行介质的飞行器建模中也存在类似问题，因此研究这个问题具有重要意义。

导致这种不一致的根本原因有两点：其一是缺乏统一的理论框架来计算和对比无旋无环流以及有黏流中的气动力和力矩。在无旋无环流中，气动力或力矩根据Lamb的理论取决于流场中的速度分布；在有黏流中，气动力或力矩一般是根据物体表面的压力和摩擦力分布积分或加权积分得到；这样人们就难以比较两种流场中同成分气动力/力矩的从属关系。其二是两种流场中气动力/力矩的成分划分不明确，一般将气动力/力矩按其与运动体当前运动参数的关联性划分为不同成分。对于轴对称体而言，与纵向非定常气动力/力矩有关的当前运动参数有俯仰角速度和迎角变化率等，而与俯仰角速度相关的气动力/力矩既可以是因俯仰角速度引起，也可以是因部分迎角变化率引起[12]。这种特点使得人们无法从实验测量数据或数值计算(CFD)结果中区分俯仰角速度引起的气动力/力矩成分和迎角变化率引起的气动力/力矩成分[10,14]，给气动力/力矩的成分划分带来困难。为了解决这两个问题，本文首先介绍能兼容Lamb无旋无环流理论的涡量矩定理，使得飞艇的气动力和力矩可以在统一的框架下计算和分析，进而明确同成分气动系数在两种流场中的从属关系；然后通过重构运动体的当前运动参数，将俯仰角速度和迎角变化率转换为俯仰角速度和平动加速度，使得按这种新运动参数划分的气动力和力矩成分是明确并可计算或测量的。进而根据这两个研究结果，提出了一种动导数与附加质量的新融合方法。文章最后通过一个实际算例来评估这种新融合方法和文献[18]中融合方法的差异对飞艇纵向扰动特性的影响，说明采用这种新方法的必要性。

1 不可压缩流中通用的气动力和力矩表达式

关于非定常气动力和力矩的计算和分析理论，最早被Lamb用于流体无黏假设系统阐述中。这种无黏流场中既不存在涡也不存在环量，因而和空气动力学中所说的无黏流和势流概念又有所不同，这里简称其为无旋无环流。与无旋无环流对应的是有旋流或涡流，涡是由黏性引起的，在强调黏性影响的情况下也将涡流称为有黏流。

Lamb[1]所介绍无旋无环流场中对应的气动力和力矩方法是基于能量守恒的原理推导得到的；因为有黏流场中存在黏性损耗，因此这种方法难以被推广到有黏流中去。运动体在有黏流场中运动时所对应的气动力和力矩可以用Wu提出的涡量矩定理[24-26]来计算和分析，这种方法是基于动量守恒方法得到的。然而在流场黏性减少为零的情况下，涡量矩定理并不能收敛于Lamb的研究结果[27]。这种不兼容性的原因之一在于涡量矩定理的推导过程中将运动体和流场视为一体。在无旋无环流的情况下，滑移边界条件将导致速度等变量分布在流体和固体所占领的空间区域上不连续，从而导致涡量矩定理的推导过程和结果不成立。涡量矩定理不能兼容无旋无环流结果的另一个原因是目前许多涡动力学理论[28-31]中的气动力和力矩表达式并没有考虑力矩参考点可动的情况；然而在飞行力学中，运动体所受气动力矩的参考点一般是固联在运动体上的，因而是运动的。考虑到这些问题和不足，重新推导运动体的气动力和力矩表达式可得到能同时兼容涡量矩定理和无旋无环流结果的气动力FA和力矩MA的公式为

(1)

(2)

式中：N=2,3为流场的维度；Si为与物面Sb相接触的流体表面，Si与Sb之间的间距为无穷小但不重合；Rf,∞是流体所占领的无穷大空间区域；n是Si的法向方向并指向Rf,∞的外部，以指向运动体为正；r=x-x0是流场中某一点相对参考点O0的位矢；ρ是流体密度；ω是流场中的涡量分布。

注1Rf,∞是流体所占领的无穷大空间区域，它不同于一个趋于无穷大的空间区域。无穷大的空间区域Rf,∞是没有外边界的，因此式(2)在使用导数矩转换(Derivative Moment Transformation， DMT)公式[28]展开体积分后，只有内边界Si上的积分。若是趋于无穷大的空间区域，因其有外边界，在使用DMT公式展开体积分后还会多出一个与内边界Si类似的外边界Se上的面积分。

注2对于无旋无环流场，式(2)可简化为

(3)

在无旋无环流中，Si上的滑移速度可以根据Hess[32]的面元法来计算。对于旋成椭球体的计算结果表明，用式(1)计算得到的气动力和力矩，与Lamb[1]和Munk[2]所得到的结果一致。

注3基于Lamb的理论，在求得物面速度势的情况下，可以计算任意外形运动体的附加质量或无旋无环流气动力/力矩。式(3)实际上提供了一种更简便的方法，在获得物面速度分布的情况下，并不需要积分求速度势就可以直接求气动力，这是式(3)的第1个优点。其次，Lamb 公式中并没有明确力矩参考点的位置，从上下文来看力矩参考点应该是体坐标系的原点或平动参考点。事实上，力矩的参考点并不一定会和平动参考点或体坐标系的原点重合，因此式(3)给出的公式更具一般性。最后，Lamb公式中的v0是指平动参考点或体坐标系原点的速度。而式(1)中的v0是指力矩参考点的速度。在一般情况下，力矩参考点的速度是可以任意指定的，并不一定会与旋转中心或体坐标系原点的速度一致。这些都表明式(3)是对Lamb理论的推广，同时也说明了Lamb公式的不足。

注4对于有黏流场而言，根据无滑移边界条件Sb和Si上的流动参数相同，不难看出式(1)中的气动力表达式与涡量矩定理[26]的气动力表达式一致。然而对比式(1)中的气动力矩表达式与涡量矩定理中的力矩表达式可发现，前者多了一项v0×Qf，∞。涡动力学理论中缺少这一项的根本原因在于推导动量矩定理的时候，并没有考虑参考点速度的影响。在飞行力学的研究中，力矩的参考点一般是固联在运动体上的，即速度不为零。因此考虑力矩参考点速度影响的结果更实用。

(4)

式(4)清楚地表明，无旋无环流和有黏流中流场的动量或动量矩享有相同的表达式。这样根据式(1)可知，运动在有黏流场和无旋无环流场中以相同速度运动时，其气动力或力矩的差异是因为二者的涡量分布不同所致。对于有黏流而言，物面上所产生的涡将在黏性的作用下逐渐扩散到流场中去形成边界层，并在对流的作用下最终形成尾迹。随着黏性的减小，扩散作用减弱，物面所产生的涡更加集中在物面附近并形成边界层。当流场黏性进一步减少到零，这个边界层就会最终缩减为滑移边界。这表明无旋无环流场只是有黏流场在流体黏性为零情况下的一个特例。因此，在无旋无环流气动力或力矩与有黏流中对应成分融合时，应当仅取有黏流中的结果。需要注意的问题是，只有相同成分的气动力或力矩才可以按上面的研究结论融合。然而有黏流和无旋无环流中的气动力或力矩的成分划分方法并不相同。无旋无环流中的气动力或力矩一般是根据运动体平动参考点速度和角速度划分，而有黏流中的气动力或力矩一般是根据运动体的迎角、迎角变化率以及角速度等参数来划分的。这样在气动力或力矩融合之前，首先需要明确这两种气动力或力矩划分结果的成分对应关系，或者说有黏流场中气动系数与无旋无环流场中附加质量所表征气动力或力矩的对应关系。

2 气动系数计算与融合

只有基于可获得的气动系数讨论融合才有意义，也只有先明确气动系数的特点，才能讨论解决气动系数融合的方法。因此本节先分别介绍气动系数在无旋无环流和有黏流中的算法和特点，然后讨论气动系数在有黏流和无旋无环流中的融合方法。

2.1 无旋无环流中的气动系数计算

根据前面的分析，无旋无环流中气动力可根据式(1)和式(3)来计算，但这个结果并不能将气动力或力矩与运动体当前的运动参数关联起来。在初始时刻流场是静止的情况下，根据面元法的理论[32]、广义的Biot-Savart定理[26]或者文献[1]中所介绍的Kirchhoff公式，均不难证明在无旋无环流的情况下Si上的滑移速度与运动体的当前运动参数成正比。根据式(4)可知，流场中的总动量和动量矩可分别表示为

(5)

式中：vD为平动参考点的速度；ωb为运动体的旋转角速度；λFv、λFω、λMv、λMω为相关系数，它们均为二阶张量。再根据式(1)可知，运动体在无旋无环流场中的气动力和力矩可表示为

(6)

式(6)表明，无旋无环流场中运动体的气动力和力矩仅取决于运动体的当前运动参数，因此可以建立气动系数的概念。这种气动系数可以取λFv、λFω、λMv、λMω在体坐标系中的投影分量，一般被称为附加质量。

对于旋成体而言，若体坐标系原点位于主轴上且Ox轴与主轴重合，则张量λFv、λFω、λMv、λMω在体坐标系中的投影具有如下简化形式：

(7)

在这种情况下，气动力和力矩在体坐标系中的投影矩阵可分别展开为

(8)

(9)

式中：u、v、w为平动参考点速度vD在体坐标系中的3个投影分量；p、q、r为当前角速度ωb在体坐标系中的3个投影分量。

根据式(6)可知，附加质量可以根据不同运动参数下的气动力/力矩差分得到。以计算λ66为例，先给运动体一个仅拥有r分量的转动角速度，然后基于面元法计算运动体气动力矩在Oz轴方向上的分量Mz；由于在速度和角速度均为零时运动体的气动力矩为零，因此根据式(6)可知λ66=Mz/r。其他附加质量系数可按类似计算[33-35]。已有的算例表明，针对椭球体，这种新计算方法与经典结果一致；针对一般外形刚体，这种新计算方法比经典方法更直观、通用。

2.2 有黏流中的气动系数计算

在有黏流中，飞行力学著作一般采用如下方式来定义气动系数：

(10)

(11)

式中：U=SL为飞艇的特征体积；j、k分别是体坐标系Oy轴和Oz轴的基矢量。对应的气动系数可分别表示为

(12)

还有一种周期性扰动是升沉运动或侧滑平动，与其加速度有关的力和力矩可表示为

(13)

且相关的气动系数可分别表示为

(14)

式中：ay、az分别表示运动体平动参考点沿着体坐标系Oy、Oz轴的平动加速度；式(14)中各子式的最后一个等号使用了泊松公式。升沉运动中运动体无俯仰角速度而只有迎角变化，因此ay实际上表示某种迎角变化率。这里又出现了新问题，那就是在测量或计算与俯仰角速度有关的气动系数mzr时，运动体的迎角也发生了变化，这是否意味着辨识结果中也包含一部分迎角变化率引起的力矩呢？这个问题的答案是肯定的，人们很早就发现，实验测量得到与俯仰角速度相关的力矩中同时包含俯仰角速度引起的阻尼力矩和迎角变化率引起的下洗延迟力矩[12]；Wang在2012年基于CFD 和傅里叶分析法计算动导数时也无法区分俯仰阻尼力矩和下洗延迟力矩这两类气动系数[14]。因此，基于数值计算结果或实验数据辨识得到与俯仰角速度r相关的俯仰力矩系数mzr并非俯仰角速度所引起的俯仰阻尼力矩系数mz,r，因为前者还多包含一部分因迎角变化率所对应的气动力矩。这种现象的存在，给气动系数的成分划分带来了困难，进而阻碍了气动系数融合方法的研究，其解决方法将在2.3节具体讨论。

2.3 动导数与附加质量的融合

动导数和附加质量分别在有黏流和无旋无环流中将飞艇的气动力和力矩与其当前运动参数关联起来，这就产生了动导数和附加质量的融合问题。根据前面所介绍的不可压缩流中气动力和力矩的通用分析理论可知，同成分气动系数融合时，应当取有黏流中的结果并同时摒弃对应的附加质量。如果能在有黏流中计算所有的动导数，那么气动系数的融合就不存在困难。从前面的气动系数计算分析中可以看出，无旋无环流中的附加质量成分要比有黏流中的动导数丰富。对于那些有黏流中无对应成分的无旋流附加质量，在融合过程中应当保留。因此研究气动系数融合的目的在于明确无旋无环流中附加质量需要保留的部分。为了实现这个目的，需要明确这两种流场中气动系数成分的划分及其对应关系。

有黏流中的气动系数一般表示为平动参考点速度的幅值vD、迎角α、侧滑角β和运动体旋转角速度ωb的函数，而无旋无环流中的气动力或力矩一般表征为平动参考点的平动速度在体坐标系中分量u、v、w和ωb的函数。为了融合有黏流和无旋无环流的气动系数，首先需要建立这两类参数之间的联系，在小迎角假设下，这种关系为

(15)

可进一步求得

(16)

(17)

接下来讨论如何将无旋无环流气动力与新当前运动参数相关联。利用泊松公式可将式(6)中的气动力表达式改写为

ωb×(λFv·vD)-ωb×(λFω·ωb)

(18)

其中：dr表示相对导数，即对矢量求导数时，只考虑坐标的变化而不考虑参考坐标系基矢量随时间的变化[39]。从式(18)最后一个等号右边可以明显看出，与运动体角速度相关的气动力为最后3项而不是最后2项。这是因为drvD/dt采用相对导数的符号，本身受运动体角速度影响；而dvD/dt采用了绝对导数的符号，不受运动体角速度影响的缘故。这样无旋无环流中与俯仰角速度q、r相关的气动力为

Fqr,p=λFv·[(qj+rk)×vD]-

(qj+rk)×(λFv·vD)-

(qj+rk)×[λFω·(qj+rk)]

(19)

另一方面，在运动体没有平动加速度，即dvD/dt=0 的情况下，由式(18)可知

(20)

Fqr1,p=Fqr,p-Fqr,p·i

(21)

式中：i为体坐标系的Ox轴基矢量。根据有黏流和无旋无环流中的气动系数关系，式(21)与式(11) 融合后应当仅取式(11)中的第1式。

由式(18)还可知，无旋无环流中与横向和法向平动加速度有关的气动力为

(22)

有黏流场中与平动加速度相关的气动系数为式(13) 中的第1式。类似地，这个式子并不含Ox轴向分量，因此它与无旋无环流中Fa,p的法向和横向分量为

Fayaz1,p=Fayaz,p-Fayaz,p·i

(23)

根据有黏流和无旋无环流中的气动系数关系，式(23) 与式(13)融合后应当仅取式(13)中的第1式。

类似地，式(6)中的气动力矩表达式可改写为

λMv·(ωb×vD)-ωb×(λMv·vD)-

ωb×(λMω·ωb)-

v0×λFv·vD-v0×λFω·ωb

(24)

Mqr,p=λMv·[(qj+rk)×vD]-

(qj+rk)×(λMv·vD)-

(qj+rk)×[λMω·(qj+rk)]-

v0×[λFω·(qj+rk)]

(25)

其法向和横向分量为

Mqr1,p=Mqr,p-Mqr,p·i

(26)

无旋无环流中与法向和横向平动加速度相关的气动力矩为

(27)

其法向和横向分量为

Mayaz1,p=Mayaz,p-Mayaz,p·i

(28)

有黏流中与俯仰角速度q、r相关且与Mqr1,p对应的气动力矩为式(11)中的第2式，与运动体横向和法向平动加速度ax、ay相关且与Mayaz1,p对应的气动力矩为式(13)中的第2式。类似地，与同样当前运动参数相关联的，即同成分气动力和力矩的融合方法是直接取有黏流的结果。因此与周期性摆动有关的气动力和力矩融合结果分别直接取式(11)和式(13)中的第2式。

3 算例分析

图1 两种气动系数融合方法对比

取平衡态为匀速直线运动，将飞艇的飞行力学模型线性化可得其纵向扰动运动模型为

(29)

式中：

xln=[ΔvD,Δα,r,Δϑ]

Mln=

(30)

aln=

(31)

(32)

对于气动系数融合方法Ⅰ，推导的结果表明

(33)

(34)

取文献[18]中飞艇的结构参数，如表1所示。

表1 飞艇的参数

对于方法Ⅱ，气动参数计算如下：

(35)

(36)

图2 纵向扰动系统的根轨迹

为了进一步分析各特征根的差异，将复数根的实部和虚部分别绘制于图3，将两个实数根分别绘制于图4。从图3可以看出，共轭复数根的实部差异较大，但这种差异随vD变化并不显著。在vD∈(5 m/s,65 m/s)的范围内，按文献[18]融合方法所引起的复数根实部相对误差约在15%～20%之间。从图3还可以看出，共轭复数根的虚部在vD较小时差异显著，随着vD增加这种差异逐渐减小。在飞行速度为5 m/s时，复数根虚部相对误差可达33%。

图3 复特征根随飞行速度的变化

图4 实特征根随飞行速度的变化

从图4可以看出，纵向扰动系统拥有一大一小两个实数根。其中绝对值较大的实数根在飞艇飞行速度较小的时候差异较大，随着飞行速度增加，这种差异逐渐减小；在飞行速度为5 m/s时，这个实数根的相对误差约为39%。绝对值较小的实数根在飞艇飞行速度较小时几乎没有差异；随着速度增加，这种差异逐渐增大，在飞行速度约为25 m/s时，这个实数根的相对误差达到最大值12%；但随着速度的进一步增加，这个实数根的相对误差又逐渐减小。

4 结 论

气动系数是保证飞行器飞行力学模型封闭的前提条件。在无旋无环流中，飞艇的气动力/力矩可以直接与其当前运动参数相关，非定常气动系数即附加质量的概念是存在的，气动力/力矩或气动系数均可根据当前运动参数明确划分。对于有黏流而言，只有一些与周期性运动相关的非定常气动力/力矩才可以与当前运动参数相关，进而定义非定常气动系数即动导数的概念。

为了研究这两种气动系数的融合问题，首先建立了能同时适用于无旋无环流和有黏流的气动力和力矩表达式；然后在同一框架下对这两种流场中的气动系数进行分析和成分划分；最后明确了这两种气动系数的融合方法。就纵向非定常气动力/力矩而言，它一般被表示为俯仰角速度和迎角变化率的函数，这时应特别注意基于实验数据或数值计算结果辨识得到与俯仰角速度相关的气动力/力矩中也含有部分与迎角变化率引起的成份；因此与迎角变化率相关的气动力/力矩应当只考虑升沉平动加速度所对应迎角变化率引起的部分。

为了将动导数与附加质量融合，无旋无环流中的气动力/力矩也应该按类似的方法划分，这样在动导数和附加质量融合时，若动导数存在，则融合结果取动导数，同时摒弃对应附加质量所表征的气动力/力矩；若动导数不存在，则取附加质量对应的结果。忽视与俯仰角速度相关气动力/力矩中含有迎角变化率引起的成份这个事实并使用融合后的气动系数分析飞艇的动态特性，将使特征根产生较大的误差。这个分析结果对于横向动导数和附加质量的融合同样适用。