飞行力学模型是研究飞艇技术的一个重要基础，建立飞行力学模型的难点在于气动力和力矩的确定和处理。在长达一百多年的飞艇技术发展过程中，许多学者对飞艇气动力和力矩的计算、测量和处理方法进行了研究。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]中的气动力和力矩表达式并没有考虑力矩参考点可动的情况；然而在飞行力学中，运动体所受气动力矩的参考点一般是固联在运动体上的，因而是运动的。考虑到这些问题和不足，重新推导运动体的气动力和力矩表达式可得到能同时兼容涡量矩定理和无旋无环流结果的气动力F_A和力矩M_A的公式为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_{\rm{A}}} = - \frac{{{\rm{d}}{\mathit{\boldsymbol{Q}}_{f,\infty }}}}{{{\rm{d}}t}}}\\ {{\mathit{\boldsymbol{M}}_{\rm{A}}} = - \frac{{{\rm{d}}{\mathit{\boldsymbol{K}}_{f,\infty }}}}{{{\rm{d}}t}} - {\mathit{\boldsymbol{v}}_0} \times {\mathit{\boldsymbol{Q}}_{f,\infty }}} \end{array}} \right. $	（1）

式中：Q_f，∞是无穷大流场中的总动量；K_f，∞是无穷大流场中的总动量矩；${\mathit{\boldsymbol{v}}_0} = {\mathit{\boldsymbol{\dot r}}_0}$是位矢为r₀的力矩参考点O₀的速度。Q_f，∞和K_f，∞的具体表达式为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Q}}_{f,\infty }} = \int_{{R_{f,\infty }}} \rho \mathit{\boldsymbol{v}}{\rm{d}}R = \frac{\rho }{{N - 1}}\int_{{R_{f,\infty }}} \mathit{\boldsymbol{r}} \times \mathit{\boldsymbol{\omega }}{\rm{d}}R + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{\rho }{{N - 1}}\int_{{S_{\rm{i}}}} \mathit{\boldsymbol{r}} \times (\mathit{\boldsymbol{v}} \times \mathit{\boldsymbol{n}}){\rm{d}}S}\\ {{\mathit{\boldsymbol{K}}_{f,\infty }} = \int_{{R_{f,\infty }}} \rho \mathit{\boldsymbol{r}} \times \mathit{\boldsymbol{v}}{\rm{d}}R = - \frac{\rho }{2}\int_{{R_{f,\infty }}} {{\mathit{\boldsymbol{r}}^2}} \mathit{\boldsymbol{\omega }}{\rm{d}}R + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{\rho }{2}\int_{{S_{\rm{i}}}} {{\mathit{\boldsymbol{r}}^2}} (\mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{v}}){\rm{d}}S} \end{array}} \right. $

（2）

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

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

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

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Q}}_{f,\infty }} = \int_{{R_{f,\infty }}} \rho \mathit{\boldsymbol{v}}{\rm{d}}R = \frac{\rho }{{N - 1}}\int_{{S_{\rm{i}}}} \mathit{\boldsymbol{r}} \times (\mathit{\boldsymbol{v}} \times \mathit{\boldsymbol{n}}){\rm{d}}S}\\ {{\mathit{\boldsymbol{K}}_{f,\infty }} = \int_{{R_{f,\infty }}} \rho \mathit{\boldsymbol{r}} \times \mathit{\boldsymbol{v}}{\rm{d}}R = \frac{\rho }{2}\int_{{S_{\rm{i}}}} {{\mathit{\boldsymbol{r}}^2}} (\mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{v}}){\rm{d}}S} \end{array}} \right. $

（3）

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

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

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

在式(1)和式(2)中，积分区域R_{f, ∞}的内边界为S_i。在无旋无环流的情况下，由于滑移边界条件，S_i和物面边界S_b上的速度并不相同而无法互相替换。这不但影响这两种流场中气动力和气动力矩表达式的统一性，也不方便某些场合的应用。在无旋无环流的情况下，滑移边界可以用一层涡层来等效。这种涡层可以用一个涡量层在厚度逼近于零的情况下来近似。在一个很薄且边界面互相平行的涡量层内，由于涡和速度的方向均与边界面平行，因此不难建立二者之间的关系为ωdσ=－n×v_Si+n×v_Sb，其中dσ为涡量层的厚度。记S_i和S_b所包含的空间区域为R₁^－且R′_{f, ∞}= R_{f, ∞}∪R₁^－，则式(2)可改写为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Q}}_{f,\infty }} = \frac{\rho }{{N - 1}}\int_{{R^\prime }_{f,\infty }} \mathit{\boldsymbol{r}} \times \mathit{\boldsymbol{\omega }}{\rm{d}}R + \frac{\rho }{{N - 1}}\int_{{S_{\rm{b}}}} \mathit{\boldsymbol{r}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (v \times n)dS}\\ {{\mathit{\boldsymbol{K}}_{f,\infty }} = - \frac{\rho }{2}\int_{{R^\prime }_{f,\infty }} {{\mathit{\boldsymbol{r}}^2}} \mathit{\boldsymbol{\omega }}{\rm{d}}R + \frac{\rho }{2}\int_{{S_{\rm{b}}}} {{\mathit{\boldsymbol{r}}^2}} (\mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{v}}){\rm{d}}S} \end{array}} \right. $

（4）

式中：R′_{f, ∞}是运动体表面S_b外部流体所占领的无穷大空间区域。R′_{f, ∞}与R_{f, ∞}的不同之处在于其内边界为S_i而不是S_b。对于无旋无环流而言，R′_{f, ∞}包含与涡层对应的滑移边界而R_{f, ∞}则并不包含这个滑移边界。

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

2 气动系数计算与融合

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

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

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

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Q}}_{f,\infty }} = {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}} + {\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}}}\\ {{\mathit{\boldsymbol{K}}_{f,\infty }} = {\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}} + {\mathit{\boldsymbol{\lambda }}_{{M_\omega }}} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}}} \end{array}} \right. $

（5）

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

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_{{\rm{A}},p}} = - \frac{{{\rm{d}}({\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}} + {\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}})}}{{{\rm{d}}t}}}\\ {{\mathit{\boldsymbol{M}}_{{\rm{A}},p}} = - \frac{{{\rm{d}}({\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}} + {\mathit{\boldsymbol{\lambda }}_{M\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}})}}{{{\rm{d}}t}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{v}}_0} \times {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}} - {\mathit{\boldsymbol{v}}_0} \times {\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}}} \end{array}} \right. $

（6）

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

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

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{\lambda }}_{Fv}} = \left[ {\begin{array}{*{20}{c}} {{\lambda _{11}}}&0&0\\ 0&{{\lambda _{22}}}&0\\ 0&0&{{\lambda _{33}}} \end{array}} \right],{\mathit{\boldsymbol{\lambda }}_{F\omega }} = \left[ {\begin{array}{*{20}{c}} 0&0&0\\ 0&0&{{\lambda _{26}}}\\ 0&{{\lambda _{35}}}&0 \end{array}} \right]}\\ {{\mathit{\boldsymbol{\lambda }}_{Mv}} = \left[ {\begin{array}{*{20}{c}} 0&0&0\\ 0&0&{{\lambda _{26}}}\\ 0&{{\lambda _{35}}}&0 \end{array}} \right],{\mathit{\boldsymbol{\lambda }}_{{M_\omega }}} = \left[ {\begin{array}{*{20}{c}} {{\lambda _{44}}}&0&0\\ 0&{{\lambda _{55}}}&0\\ 0&0&{{\lambda _{66}}} \end{array}} \right]} \end{array}} \right. $

（7）

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

$ {\mathit{\boldsymbol{F}}_{{\rm{A}},p}} = - \left[ \begin{array}{l} {\lambda _{11}}\dot u + {\lambda _{35}}{q^2} - {\lambda _{26}}{r^2} - {\lambda _{22}}vr + {\lambda _{33}}wq\\ {\lambda _{26}}\dot r + {\lambda _{22}}\dot v + {\lambda _{11}}ur - {\lambda _{33}}w{\rm{p }} - {\lambda _{35}}pq\\ {\lambda _{35}}\dot q + {\lambda _{33}}\dot w - {\lambda _{11}}uq + {\lambda _{22}}vp + {\lambda _{26}}pr \end{array} \right] $

（8）

$ {\mathit{\boldsymbol{M}}_{{\rm{A}},p}} = - \left[ \begin{array}{l} {\lambda _{44}}\dot p - ({\lambda _{22}} + {\lambda _{33}})vw + ({\lambda _{35}} + {\lambda _{26}})vq - ({\lambda _{26}} + {\lambda _{35}})wr{\rm{ }} - ({\lambda _{55}} - {\lambda _{66}})qr\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\lambda _{55}}\dot q + {\lambda _{35}}\dot w + ({\lambda _{11}} - {\lambda _{33}})uw - {\lambda _{26}}vp - {\lambda _{35}}uq + ({\lambda _{44}} - {\lambda _{66}})pr\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\lambda _{66}}\dot r + {\lambda _{26}}\dot v + ({\lambda _{22}} - {\lambda _{11}})vu + {\lambda _{26}}ur + {\lambda _{35}}wp + ({\lambda _{55}} - {\lambda _{44}})qp \end{array} \right] $

（9）

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

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

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

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

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_v} = \bar QS({C_x}{\mathit{\boldsymbol{e}}_1} + {C_y}{\mathit{\boldsymbol{e}}_2} + {C_z}{\mathit{\boldsymbol{e}}_3})}\\ {{\mathit{\boldsymbol{M}}_v} = \bar QSL({m_x}{\mathit{\boldsymbol{e}}_1} + {m_y}{\mathit{\boldsymbol{e}}_2} + {m_z}{\mathit{\boldsymbol{e}}_3})} \end{array}} \right. $	（10）

式中：Q=ρv_D²/2为动压头；L为飞行器的特征长度；S为飞行器的特征面积；如果e₁、e₂、e₃分别是体坐标系的3个基矢量，则C_x、C_y、C_z分别称为轴向力系数、法向力系数和横向力系数；m_x、m_y、m_z分别称为滚转力矩、偏航力矩和俯仰力矩。如果将气动力在速度坐标系中投影，这时C_x、C_y、C_z分别称为阻力、升力和侧力。根据上述定义，在测量得到气动力/力矩后，即可直接计算对应的气动系数。

这种定义方法只适用于飞行器匀速直线运动且气动力和力矩处于稳态的情况。当飞行器具有加速度和角速度时，上面的许多力和力矩系数需要进一步分解。例如俯仰力矩系数m_z可进一步分解为迎角引起的静力矩系数m_z0、俯仰角速度r引起的俯仰阻尼力矩系数m_{z, r}以及迎角变化率引起的下洗延迟力矩系数m_{z,$\dot a$} 等。然而，当运动体在有黏流场中非定常运动的情况下，其气动力或力矩是难以与其当前运动参数严格关联的，也就难以严格定义气动系数的概念。举个简单的例子，把手放在静止的池塘里，手所受的水动力为零。现让手剧烈运动后突然停止，这时候手的速度仍为零但所受的水动力并不为零。这个例子有力地说明了有黏流场中运动体的非定常气动力/力矩并不能与其速度严格关联。究其根本原因，在于有黏流场受动力学方程(Navier-Stokes方程或涡量传输方程)控制而具备记忆或者受流场的历史运动信息影响，这就使得在其中运动物体所受的气动力或力矩，不仅取决于其当前运动参数，也同时与运动体的历史运动信息有关。

由于有黏流中不存在严格气动系数的概念，这给飞行器飞行力学模型的建立带来了困难。较早尝试解决这个问题的学者是Bryan^[36]，他于1911年在“缓慢运动假设”下，将气动力和力矩在平衡态附近进行泰勒级数展开，从而将运动体的扰动气动力和力矩与其扰动运动关联起来并得到了相应的气动系数。这种方法及其改进形式^[37]被许多飞行力学理论一直沿用至今。基于这种方法的改进形式，俯仰力矩可被分解为由迎角α引起的静力矩，由迎角变化率$\dot a$引起的下洗延迟力矩，由俯仰角速度r引起的俯仰阻尼力矩。后来人们也发展了一些计算俯仰阻尼力矩^[13]和下洗延迟力矩^[38]的工程估算方法，但工程算法依赖于经验参数，并且难以根据实验数据来确定。Bryan引入缓慢运动假设的目的在于用忽略历史运动信息影响的办法来实现气动力/力矩与运动体当前运动参数的关联，在某些历史运动信息已知并且其对当前气动力和力矩的影响可计算的情况下，不需要忽略历史运动信息的影响也可将气动力和力矩与当前运动参数相关。一个典型的例子是运动体以一定规律做周期性运动的情况，虽然历史运动信息会对当前气动力和力矩产生影响，但这种情况下的当前气动力和力矩仍可与当前运动参数相关联。在这个意义下仍可定义气动系数的概念，进而探讨它们在有黏流中的具体计算和测量方法。仍然以俯仰力矩为例，它可以让运动体在来流中绕平动参考点周期性俯仰旋转，先基于CFD方法计算或实验方法测量得到气动力矩，然后就可采用Wang介绍的傅里叶分析法将其分别与迎角α、俯仰角速度r相关^[14]。进一步按前面给出的公式归一化即可得到静力矩系数m_z0和与俯仰摆动过程中旋转角速度有关的气动系数m_zr。这样，俯仰或偏航摆动过程中与俯仰或偏航角速度有关的力和力矩可表示为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_{qr,v}} = \bar QS({C_{yr}}\mathit{\boldsymbol{j}} + {C_{zq}}\mathit{\boldsymbol{k}})}\\ {{\mathit{\boldsymbol{M}}_{qr,v}} = \bar QU({m_{yq}}\mathit{\boldsymbol{j}} + {m_{zr}}\mathit{\boldsymbol{k}})} \end{array}} \right. $	（11）

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

$ \left\{ {\begin{array}{*{20}{l}} {{C_{yr}} = C_y^rr,{C_{zq}} = C_z^qq}\\ {{m_{yq}} = m_y^qq,{m_{zr}} = m_z^rr} \end{array}} \right. $	（12）

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

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_{{a_y}{a_z},v}} = \bar QS({C_{y{a_y}}}\mathit{\boldsymbol{j}} + {C_{z{a_z}}}\mathit{\boldsymbol{k}})}\\ {{\mathit{\boldsymbol{M}}_{{a_y}{a_z},v}} = \bar QU({m_{y{a_z}}}\mathit{\boldsymbol{j}} + {m_{z{a_y}}}\mathit{\boldsymbol{k}})} \end{array}} \right. $	（13）

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

$ \left\{ \begin{array}{l} {C_{y{a_y}}} = C_y^{{a_y}}{a_y} = C_y^{{a_y}}(\dot v + ru - pw)\\ {C_{z{a_z}}} = C_z^{{a_z}}{a_z} = C_z^{{a_z}}(\dot w + pv - qu)\\ {m_{y{a_z}}} = m_y^{{a_z}}{a_z} = m_y^{{a_z}}(\dot w + pv - qu)\\ {m_{z{a_y}}} = m_z^{{a_y}}{a_y} = m_z^{{a_y}}(\dot v + ru - pw) \end{array} \right. $	（14）

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

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

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

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

$ \left\{ {\begin{array}{*{20}{l}} {{v_{\rm{D}}} = \sqrt {{u^2} + {v^2} + {w^2}} }\\ {\alpha = - {\rm{arctan}}\frac{v}{u}}\\ {\beta = {\rm{arcsin}}\frac{w}{{{v_{\rm{D}}}}}} \end{array}} \right. $	（15）

可进一步求得

$ \left\{ {\begin{array}{*{20}{l}} {\dot \alpha = \frac{{v\dot u - \dot vu}}{{{u^2} + {v^2}}} \approx - \dot v/{v_{\rm{D}}}}\\ {\dot \beta = \frac{{\dot w{v_{\rm{D}}} - w{{\dot v}_{\rm{D}}}}}{{{v_{\rm{D}}}\sqrt {{w^2} + v_{\rm{D}}^2} }} \approx \dot w/{v_{\rm{D}}}} \end{array}} \right. $	（16）

基于式(15)和式(16)，可以将有黏流和无旋无环流中的当前运动参数联系起来。但这样还不足以完全明确与当前运动参数相关联的气动系数，这是因为有黏流中辨识得到与俯仰角速度相关的气动系数中耦合了与迎角变化率有关的气动系数。注意到${a_y} = \dot v + ru - pw, \;{a_z} = \dot w + pv - qu$，故有

$ \left\{ {\begin{array}{*{20}{l}} {\dot \alpha = {{\dot \alpha }_1} + {{\dot \alpha }_2}}\\ {\dot \beta = {{\dot \beta }_1} + {{\dot \beta }_2}} \end{array}} \right. $	（17）

式中：${{\dot \alpha }_1} = - {a_y}/{v_{\rm{D}}}$表示运动体升沉平动所对应的迎角；${{\dot \alpha }_2} = (ru - pw)/{v_{\rm{D}}}$表示运动体俯仰转动过程中对应的迎角变化率；${{\dot \beta }_1}$、${{\dot \beta }_2}$的含义类似。这样运动体的当前运动参数可不用v_D、α、β、$\dot \alpha 、\dot \beta $、p、q、r来表示，而是采用一组新的参数v_D、α、β、p、q+ ${\dot \beta }$₂、${\dot \beta }$₁、r+ $\dot a$₂、$\dot a$₁来表示。这样处理的好处在于有黏流中由q+ ${\dot \beta }$₂、${\dot \beta }$₁、r+ ${\dot \alpha }$₂、${\dot \alpha }$₁引起的气动系数可以基于实验数据或CFD计算结果辨识得到；无旋无环流中的气动力和力矩因为有解析表达式，不难确定由q+ ${\dot \beta }$₂、${\dot \beta }$₁、r+${\dot \alpha }$₂、${\dot \alpha }$₁引起的气动系数。这样由新当前运动参数引起的气动系数在有黏流和无旋无环流中均可获得，进而也就可以解决气动系数的融合问题。

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

$ \begin{array}{l} {\mathit{\boldsymbol{F}}_{{\rm{A}},p}} = - {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot \frac{{{{\rm{d}}_r}{\mathit{\boldsymbol{v}}_{\rm{D}}}}}{{{\rm{d}}t}} - {\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot \frac{{{{\rm{d}}_r}{\mathit{\boldsymbol{\omega }}_{\rm{b}}}}}{{{\rm{d}}t}} - {\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times {\mathit{\boldsymbol{Q}}_{f,\infty }} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} - {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot \frac{{{\rm{d}}{\mathit{\boldsymbol{v}}_{\rm{D}}}}}{{{\rm{d}}t}} - {\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot \frac{{{\rm{d}}{\mathit{\boldsymbol{\omega }}_{\rm{b}}}}}{{{\rm{d}}t}} + {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot ({\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times {\mathit{\boldsymbol{v}}_{\rm{D}}}) - }\\ {{\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times ({\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}}) - {\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times ({\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}})} \end{array} \end{array} $

（18）

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

$ \begin{array}{l} {\mathit{\boldsymbol{F}}_{qr,p}} = {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot [(q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}}) \times {\mathit{\boldsymbol{v}}_{\rm{D}}}] - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {(q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}}) \times ({\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}}) - }\\ {(q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}}) \times [{\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot (q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}})]} \end{array} \end{array} $

（19）

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

$ - {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot \frac{{{{\rm{d}}_r}{\mathit{\boldsymbol{v}}_{\rm{D}}}}}{{{\rm{d}}t}} = {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot ({\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times {\mathit{\boldsymbol{v}}_{\rm{D}}}) $	（20）

式中：等号左边是表示因坐标系旋转而导致v_D在体坐标系中投影变化所对应的气动力。由于$\dot \alpha \approx - \dot v/{v_{\rm{D}}}, \dot \beta \approx \dot w/{v_{\rm{D}}}, {{\rm{d}}_r}{\mathit{\boldsymbol{v}}_{\rm{D}}}/{\rm{d}}t$实质上表征无平动加速度情况下迎角和侧滑角的变化率，因此式(20)左边表征无平动加速度情况下与迎角变化率$\dot \alpha $₂和侧滑角变化率$\dot \beta $₂有关的气动力。式(20)还表明，这种与迎角和侧滑角变化率有关的气动力，可以转换为角速度ω_b的函数。这与Cowley和Glauert在文献[12]中针对有黏流的研究结果是一致的，即运动体俯仰运动时，不仅存在俯仰角速度所引起的气动力和力矩，还存在迎角变化所引起的气动力和力矩；但二者均可视为与俯仰角速度相关的气动力和力矩。有黏流中与俯仰角速度q、r相关的动态气动力为式(11)中的第1式，这个式子并不含Ox轴向分量，因此它与无旋无环流中F_{qr, p}的法向和横向分量为

$ {\mathit{\boldsymbol{F}}_{qr1,p}} = {\mathit{\boldsymbol{F}}_{qr,p}} - {\mathit{\boldsymbol{F}}_{qr,p}} \cdot \mathit{\boldsymbol{i}} $	（21）

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

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

$ \begin{array}{l} {\mathit{\boldsymbol{F}}_{{a_y}{a_z},p}} = - {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot \frac{{{\rm{d}}(v\mathit{\boldsymbol{j}} + w\mathit{\boldsymbol{k}})}}{{{\rm{d}}t}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot (\dot v\mathit{\boldsymbol{j}} + \dot v\mathit{\boldsymbol{k}}) - {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot [{\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times (\dot v\mathit{\boldsymbol{j}} + \dot w\mathit{\boldsymbol{k}})] \end{array} $

（22）

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

$ {\mathit{\boldsymbol{F}}_{{a_y}{a_z}1,p}} = {\mathit{\boldsymbol{F}}_{{a_y}{a_z},p}} - {\mathit{\boldsymbol{F}}_{{a_y}{a_z},p}} \cdot \mathit{\boldsymbol{i}} $	（23）

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

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

$ \begin{array}{l} {\mathit{\boldsymbol{M}}_{{\rm{A}},p}} = - {\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot \frac{{{\rm{d}}{\mathit{\boldsymbol{v}}_{\rm{D}}}}}{{{\rm{d}}t}} - {\mathit{\boldsymbol{\lambda }}_{M\omega }} \cdot \frac{{{\rm{d}}{\mathit{\boldsymbol{\omega }}_{\rm{b}}}}}{{{\rm{d}}t}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot ({\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times {\mathit{\boldsymbol{v}}_{\rm{D}}}) - {\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times ({\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}}) - }\\ {{\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times ({\mathit{\boldsymbol{\lambda }}_{M\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}}) - }\\ {{\mathit{\boldsymbol{v}}_0} \times {\mathit{\boldsymbol{\lambda }}_{Fv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}} - {\mathit{\boldsymbol{v}}_0} \times {\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot {\mathit{\boldsymbol{\omega }}_{\rm{b}}}} \end{array} \end{array} $

（24）

根据式(24)可知，无旋无环流中与俯仰角速度q、r相关的气动力矩包含q、r引起的气动力矩和$\dot a$₂、$\dot \beta $₂引起的力矩，即

$ \begin{array}{l} {\mathit{\boldsymbol{M}}_{qr,p}} = {\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot [(q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}}) \times {\mathit{\boldsymbol{v}}_{\rm{D}}}] - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {(q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}}) \times ({\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot {\mathit{\boldsymbol{v}}_{\rm{D}}}) - }\\ {(q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}}) \times [{\mathit{\boldsymbol{\lambda }}_{M\omega }} \cdot (q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}})] - }\\ {{\mathit{\boldsymbol{v}}_0} \times [{\mathit{\boldsymbol{\lambda }}_{F\omega }} \cdot (q\mathit{\boldsymbol{j}} + r\mathit{\boldsymbol{k}})]} \end{array} \end{array} $

（25）

其法向和横向分量为

$ {\mathit{\boldsymbol{M}}_{qr1,p}} = {\mathit{\boldsymbol{M}}_{qr,p}} - {\mathit{\boldsymbol{M}}_{qr,p}} \cdot \mathit{\boldsymbol{i}} $	（26）

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

$ \begin{array}{l} {\mathit{\boldsymbol{M}}_{{a_y}{a_z},p}} = - {\mathit{\boldsymbol{\lambda }}_{Mv}}\frac{{{\rm{d}}(v\mathit{\boldsymbol{j}} + w\mathit{\boldsymbol{k}})}}{{{\rm{d}}t}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot (\dot v\mathit{\boldsymbol{j}} + \dot w\mathit{\boldsymbol{k}}) - {\mathit{\boldsymbol{\lambda }}_{Mv}} \cdot [{\mathit{\boldsymbol{\omega }}_{\rm{b}}} \times (v\mathit{\boldsymbol{j}} + w\mathit{\boldsymbol{k}})] \end{array} $

（27）

其法向和横向分量为

$ {\mathit{\boldsymbol{M}}_{{y^a}{z^1},p}} = {\mathit{\boldsymbol{M}}_{{a_y}{a_z},p}} - {\mathit{\boldsymbol{M}}_{{a_y}{a_z},p}} \cdot \mathit{\boldsymbol{i}} $	（28）

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

经过上面处理后，无黏流中的气动力/力矩被分解为与俯仰或偏航旋转运动有关的部分和升沉或侧滑平动运动有关的部分。前者不仅包含与角速度有关的气动力/力矩，还包含因俯仰或偏航运动而引起迎角或侧滑角变化率所对应的部分气动力/力矩；后者仅包含与升沉或侧滑所引起迎角或侧滑角变化率有关的气动力/力矩。这两个结果刚好和前面所介绍有黏流气动系数测量或计算方法得到的结果对应，双方在融合过程中直接取有黏流的结果即可。根据这个分析结果，文献[13]中同时保留有黏流中的俯仰阻尼力/力矩和无旋无环流对应结果的融合方法相当于重复计算俯仰阻尼力/力矩。文献[19]中摒弃有黏流中动导数和无旋无环流中对应结果的融合方法，将使得飞艇俯仰运动模拟过程中不存在阻尼而无法收敛。文献[20-23]仅采用无旋无环流中的对应结果，没有考虑黏性对动导数的影响。文献[18]在采用实验测量得到的与俯仰角速度相关气动系数的同时，又保留了无旋无环流中与迎角变化率有关的气动系数，这相当于重复计算了与$\dot \alpha $₂或$\dot \beta $₂有关的气动系数，将对扰动系统的动态特性分析产生影响。

3 算例分析

为了评估这种新融合方法的必要性，分析两种不同气动系数融合方法对飞艇纵向扰动运动的影响。方法Ⅰ是本文建议的方法，即把动导数分为两部分，一部分与俯仰角速度有关，另外一部分与升沉平动过程中的迎角变化率$\dot \alpha $₁有关。方法Ⅱ采用文献[18]的方法，把动导数分解为与俯仰角速度有关的部分和与迎角变化率$\dot \alpha $有关的部分。与角速度有关的部分采用了基于俯仰摆动实验数据辨识得到的结果，即同时包含俯仰角速度r和俯仰过程中迎角变化率$\dot \alpha $₂引起的气动系数。与迎角变化率有关的部分，采用无旋无环流理论中的结果；方法Ⅰ采用与dv/dt相关的部分，即与$\dot \alpha $₁相关的部分；方法Ⅱ采用与d_rv/dt相关的部分，即迎角变化同时包含由俯仰运动所引起的部分$\dot \alpha $₂和升沉运动所引起的部分$\dot \alpha $₁。以升力系数为例，二者之间的区别如图 1所示，上半部分表示无旋无环流中的升力系数分解，下半部分表示有黏流中的升力系数分解；双箭头表示两种成分是对应的。从图 1可以明显看出，方法Ⅱ比方法Ⅰ重复计算了与$\dot \alpha $₂相关的气动系数。

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

$ {\mathit{\boldsymbol{M}}_{{\rm{ln}}}}{\mathit{\boldsymbol{\dot x}}_{{\rm{ln}}}} = {\mathit{\boldsymbol{a}}_{{\rm{ln}}}}{\mathit{\boldsymbol{\dot x}}_{{\rm{ln}}}} + {\mathit{\boldsymbol{b}}_{{\rm{ln}}}}{\mathit{\boldsymbol{u}}_{{\rm{ln}}}} $	（29）

式中：

$ \begin{array}{l} \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{x}}_{{\rm{ln}}}} = [\Delta {v_{\rm{D}}},\Delta \alpha ,r,\Delta \vartheta ]}\\ {{\mathit{\boldsymbol{M}}_{{\rm{ln}}}} = } \end{array}\\ \left[ {\begin{array}{*{20}{c}} {m + {\lambda _{11}}}&0&{m{c_y}}&0\\ 0&{ - m{v_{\rm{D}}} - \bar QSC_y^{\dot \alpha }}&{ - m{c_x} - \bar QSC_y^{\dot r}}&0\\ {m{c_y}}&{m{c_x}{v_{\rm{D}}} - \bar QUm_z^{\dot \alpha }}&{ - \bar QUm_z^{\dot r} + {J_z}}&0\\ 0&0&0&{m{v_D}} \end{array}} \right] \end{array} $

（30）

$ \begin{array}{l} {\mathit{\boldsymbol{a}}_{{\rm{ln}}}} = \\ \left[ {\begin{array}{*{20}{c}} { - (\rho {v_{\rm{D}}}{C_{x0}} + 3\bar QC_x^{{v_{\rm{D}}}})S + 2{T^v}{\rm{cos}}{\mu _{\rm{e}}}}&0&{m{v_{\rm{D}}}{\alpha _{\rm{e}}}}&{ - (mg - \rho U){\rm{cos}}{\vartheta _{\rm{e}}}}\\ {2{T^v}{\rm{sin}}{\mu _{\rm{e}}}}&{\bar QSC_y^\alpha }&{\bar QSC_y^r + m{v_{\rm{D}}}}&{(mg - \rho U){\rm{sin}}{\vartheta _{\rm{e}}}}\\ {2{T^v}{d_x}{\rm{sin}}{\mu _{\rm{e}}} - 2{T^v}{d_y}{\rm{cos}}{\mu _{\rm{e}}}}&{\bar QUm_z^\alpha }&{\bar QUm_z^r + m{v_{\rm{D}}}{c_x} - m{v_{\rm{D}}}{c_y}{\alpha _{\rm{e}}}}&{\begin{array}{*{20}{l}} {(mg{c_y} - \rho U{b_y}){\rm{cos}}{\vartheta _{\rm{e}}} - }\\ {( - mg{c_x} + \rho U{b_x}){\rm{sin}}{\vartheta _{\rm{e}}}} \end{array}}\\ 0&0&{m{v_{\rm{D}}}}&0 \end{array}} \right] \end{array} $

（31）

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{u}}_{{\rm{ln}}}} = {[\Delta {\delta _z},\Delta p]^{\rm{T}}}\\ {\mathit{\boldsymbol{b}}_{{\rm{ln}}}} = \left[ {\begin{array}{*{20}{c}} 0&{2{T^p}{\rm{cos}}{\mu _{\rm{e}}}}\\ {\bar QSC_y^{{\delta _z}}}&{2{T^p}{\rm{sin}}{\mu _{\rm{e}}}}\\ {\bar QUm_z^{{\delta _z}}}&{2{T^p}{d_x}{\rm{sin}}{\mu _{\rm{e}}} - 2{T^p}{d_y}{\rm{cos}}{\mu _{\rm{e}}}}\\ 0&0 \end{array}} \right] \end{array} \right. $

（32）

其中：Q=ρV_D²表示动压头；U、S别表示飞艇的体积和特征面积；Δv_D、Δα、r、Δϑ分别表示平动速度、迎角、俯仰角速度和俯仰角的扰动值；T和μ_e分别表示螺旋桨合推力的大小和方向；d_x、d_y别表示推力的作用点，T^v、T^p分别表示推力大小随飞行速度和操纵杆位置的变化率；c_x, c_y, c_z=0 m表示重心在体坐标系中的位置；b_x, b_y, b_z=0 m表示体心在体坐标系中的位置；mg表示重力。这个线性化结果和文献[40]的纵向扰动模型是一致的，只不过文献[40]采用的扰动变量是[Δu, Δv, r, Δϑ]^T，而这里采用的扰动变量是[Δv_D, Δα, r, Δϑ]^T。

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

$ \left\{ {\begin{array}{*{20}{l}} {C_y^r = C_y^{\bar r}l/{v_{\rm{D}}} + C_y^{{a_y}}{v_{\rm{D}}}}\\ {C_y^{\dot \alpha } = - C_y^{{a_y}}{v_{\rm{D}}},C_y^{\dot r} = C_y^{\dot {\bar r}}l/{v_{\rm{D}}}}\\ {m_z^r = m_z^{{\bar r}}l/{v_{\rm{D}}} + m_z^a{v_{\rm{D}}}}\\ {m_z^{\dot \alpha } = - m_z^{{a_y}}{v_{\rm{D}}},m_z^{\dot r} = m_z^{\dot {\bar r}}l/{v_{\rm{D}}}} \end{array}} \right. $	（33）

式中：l=U^1/3表示飞艇的特征长度；r=rl/v_D；q=ql/v_D。对于气动系数融合方法Ⅱ，推导的结果表明式(33)应改为

$ \left\{ {\begin{array}{*{20}{l}} {C_y^r = C_y^{\bar r}l/{v_{\rm{D}}},C_y^{\dot \alpha } = - C_y^{{a_y}}{v_{\rm{D}}},C_y^{\dot r} = C_y^{\dot {\bar r}}l/{v_{\rm{D}}}}\\ {m_z^r = m_z^{\bar r}l/{v_{\rm{D}}},m_z^{\dot \alpha } = - m_z^{{a_y}}{v_{\rm{D}}},m_z^{\dot r} = m_z^{\dot {\bar r}}l/{v_{\rm{D}}}} \end{array}} \right. $	（34）

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

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

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{\lambda _{11}} = {K_{11}}\rho U = 35.34{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{kg/}}{{\rm{m}}^{\rm{3}}}}\\ {C_y^{\dot \alpha } = - {v_{\rm{D}}}C_y^{{a_y}} = - {v_{\rm{D}}}\frac{{ - {K_{22}}\rho U}}{{\rho {U^{2/3}}v_{\rm{D}}^2/2}} = 1.338{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}} \end{array}\\ \begin{array}{*{20}{l}} {C_y^{\dot r} = \frac{{ - {K_{26}}\rho {U^{4/3}}}}{{\rho {U^{2/3}}v_{\rm{D}}^2/2}} = 0.098{\kern 1pt} {\kern 1pt} {\kern 1pt} 95{\kern 1pt} {\kern 1pt} {\kern 1pt} {{\rm{s}}^{\rm{2}}}}\\ {m_z^{\dot \alpha } = - {v_{\rm{D}}}m_z^{{a_y}} = - {v_{\rm{D}}}\frac{{ - {K_{62}}\rho {U^{4/3}}}}{{\rho Uv_{\rm{D}}^2/2}} = - 0.137{\kern 1pt} {\kern 1pt} {\kern 1pt} 8{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}} \end{array}\\ \begin{array}{*{20}{l}} {m_z^{\dot r} = \frac{{ - {K_{66}}\rho {U^{5/3}}}}{{\rho Uv_{\rm{D}}^2/2}} = - 0.485{\kern 1pt} {\kern 1pt} {\kern 1pt} 5{\kern 1pt} {\kern 1pt} {\kern 1pt} {{\rm{s}}^2}}\\ {{C_{x0}} = 0.055{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{m/s}},C_x^{{v_{\rm{D}}}} = - 0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 4{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{m/s}}} \end{array}\\ \begin{array}{*{20}{l}} {C_y^\alpha = 0.82,m_z^a = 0.52}\\ {C_y^r = C_y^{\bar r}\frac{l}{{{v_{\rm{D}}}}} = 0.581{\kern 1pt} {\kern 1pt} {\kern 1pt} 5{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}}\\ {m_z^r = m_z^{\bar r}\frac{l}{{{v_{\rm{D}}}}} = - 0.7538{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}} \end{array} \end{array} \right. $

（35）

对于方法Ⅰ，式(35)中的C_y^r和m_z^r应分别改写为

$ \left\{ {\begin{array}{*{20}{l}} {C_y^r = C_y^{\bar r}\frac{l}{{{v_{\rm{D}}}}} - C_y^{\dot \alpha } = - 0.756{\kern 1pt} {\kern 1pt} {\kern 1pt} 5{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}}\\ {m_z^r = m_z^{\bar r}\frac{l}{{{v_{\rm{D}}}}} - m_z^{\dot \alpha } = - 0.616{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}} \end{array}} \right. $	（36）

纵向系统(29)的动态特性主要取决于其特征矩阵A=M_ln^－1a_ln的特征根。根据上面所罗列的数据并使v_D在5~65 m/s之间变化，采用这两种气动系数融合方法得A的特征根随飞行速度v_D变化的根轨迹如图 2所示。

从图 2可以看出，纵向扰动系统有4个特征根；其中一对复数根和两个实数根。这个计算结果和文献[40]中所介绍典型飞艇特征根的结果是一致的。从图 2还可以看出，不同的动导数融合

方法对共轭复数根有较大的影响；采用本文介绍的气动系数融合方法共轭复数根所对应的实部绝对值减小，这相当于减小了系统的阻尼，这个结果是容易理解的；在文献[18]处理方法中，同时在m_z^r和m_z^{$\dot \alpha $}中考虑了俯仰摆动过程中与$\dot \alpha $₂有关的气动力矩，这相当于增加了一部分俯仰阻尼力矩，因而所对应的纵向扰动系统阻尼增加。

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

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

4 结论

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

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

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