气动导数作为描述飞行器机动飞行和受扰动时气动特性的关键性气动参数，在飞行器气动性能、控制系统和总体设计中扮演着非常重要的作用^[1-4]。在传统的飞行动力学相关问题的研究中，气动力的数据往往基于小扰动线性叠加原理计算出来，在这种准定常假设情况下，气动力仅仅表示为瞬时飞行状态参数的函数，并且可以以一种简单的解析函数关系式表示出来^[2-5]。

但是，现代飞行器的飞行包线普遍向大迎角区域扩展，在大迎角下飞机机动飞行产生的三维非定常分离流和涡流使得空气动力呈现高度非线性特性，气动力和力矩不仅依赖于瞬时迎角、侧滑角、姿态角等参数, 而且与它们的时间历程有关, 因此原来使用的低阶线性叠加模型将不再适用^[5-6]。同时，由于机动飞行状态涵盖了较大的迎角、侧滑角、角速率的变化范围，如果采用风洞实验或是数值计算模拟，其时间成本和经济成本都难以接受^[7-10]。因此有必要建立起较大飞行包线内普适性较好的的非定常气动力模型^{[1, 4]}。

Etkin模型是目前动导数求解时最常用的一种非定常气动力模型，Etkin模型物理意义明确，考虑了时间历史效应对气动导数的影响^[3]。但是，在非定常气动力建模时，该模型中的各项气动导数对不同运动形式的非定常气动力的影响规律和适用程度尚不清楚。为此，本文结合Etkin气动力模型，研究了气动力关于迎角的一阶和二阶导数在气动力模型的作用，希望能精确地重构出翼型单自由度或是耦合强迫运动过程中的非定常气动力，为未来发展高效的、可靠的气动力模型提供参考数据。

1 强迫运动非定常气动力模型

本文的计算采用课题组自己开发的柔性体动力学问题求解软件GMFlow^[11-13]，其中流场求解部分采用基于SA模型的有限体积法^[13]，强迫运动时的网格变形方法为弹簧网格变形方法^[14-16]。为了验证求解方法的正确性，首先计算了标准算例NACA0012翼型强迫俯仰运动的非定常气动力变化情况，将计算结果与文献中的计算结果和实验结果对比，对比结果见文献[13]。

俯仰运动的运动规律可以描述为^[15]：

$ \alpha (t) = {\alpha _0} + A\sin (\omega t) = {\alpha _0} + A\sin (2{\rm{ \mathsf{ π} }}ft) $

(1)

式中：α₀是初始位置处的迎角；A是简谐振动的振幅；ω是简谐振动的圆频率；f是简谐振动的频率。

本文定义减缩频率为：

$ k = \frac{{\omega {\rm{C}}}}{{2{V_\infty }}} $

(2)

式中C是翼型的弦长。在本文中，强迫运动时自由来流的马赫数为0.755，翼型弦长为1.0 m，强迫运动的减缩频率为0.081 4。俯仰运动的初始迎角为0.016°，俯仰振幅为2.51°。

图 1(a)是强迫俯仰运动时的升力系数和关于1/4弦点的俯仰力矩系数随时间的变化曲线，图中计算了3个周期的气动力，由图可知，在第1个计算周期的初始阶段，计算的结果收敛性较差，这主要是由于定常计算的步数不足。从第2个周期开始，力和力矩系数已经达到了较好的谐振性，可以认为计算结果已经收敛。因此，为了减小计算量，本文的所有强迫运动过程都只计算了3个运动周期。

图 1(b)是翼型强迫沉浮运动时的力和力矩系数变化情况，其运动规律为：

$ z(t) = {z_0} + {z_m}\sin (\omega t) $

(3)

式中：z₀=0是初始位置处的纵向位移；z_m=0.1 m是沉浮运动的振幅。

考虑洗流影响，在沉浮运动的任一时刻，瞬时迎角为：

$ \alpha (t) = {\alpha _0} - \omega {z_m}\cos (\omega t)/{V_\infty } $

(4)

图 1(c)是翼型强迫俯仰/沉浮耦合运动时的力和力矩系数变化情况，其运动规律为式(1)和式(3)叠加。对比图 1可知，虽然耦合运动形式是俯仰和沉浮运动的叠加，但是耦合运动的气动力和力矩并不等于俯仰运动和沉浮运动的简单叠加，这也说明了翼型强迫运动时气动力的非线性迟滞特性比较复杂，并不是简单的线性叠加关系。

1.1 一阶气动模型

根据Etkin气动力模型^[2-3]，强迫运动过程中的非定常气动力可以表示为：

$ \begin{array}{*{20}{c}} {\Delta {C_j} = {C_j} - {C_{j0}} = {C_{j\alpha }}\Delta \alpha + }\\ {{C_{j\dot \alpha }}\left( {\frac{C}{{2{V_\infty }}}} \right)\Delta \dot \alpha + {C_{jq}}\left( {\frac{C}{{2{V_\infty }}}} \right)\Delta q} \end{array} $

(5)

式中C_j0是平衡位置处的力系数或是力矩系数。由于式(5)中${{C}_{j\dot{\alpha }}}$和C_jq的量纲相同，都是空气动力系数对角度随时间一阶变化率的导数，所以在本文中称式(5)为一阶气动力模型。根据强迫俯仰运动时运动规律可知：

$ \dot \alpha (t) = q = A\omega {\kern 1pt} {\kern 1pt} \cos (\omega t) $

(6)

俯仰运动的非定常气动力可以表示为：

$ \Delta {C_j} = A\sin (\omega t) \cdot {C_{j\alpha }} + kA(\cos (\omega t) - 1)\left( {{C_{j\dot \alpha }} + {C_{jq}}} \right) $

(7)

此处需要注意，由于ΔC_j是相对于初始位置的变化量，因此右侧是(cos(ωt)-1)而不是cos(ωt)。所以：

$ {C_{j\dot \alpha }} + {C_{jq}} = \left( {\int\limits_{{T_n}}^{{T_{n + 1}}} \Delta {{\rm{C}}_j}\cos (\omega t){\rm{d}}t} \right)/\left( {kA\frac{{\rm{ \mathsf{ π} }}}{\omega }} \right) $

(8)

根据强迫沉浮运动时运动规律可知：

$ \dot \alpha (t) = {\omega ^2}{z_m}\sin (\omega t)/{V_\infty } $

(9)

沉浮运动的非定常气动力可以表示为：

$ \Delta {C_j} = {C_{j\alpha }}\frac{{ - \omega {z_m}\cos (\omega t)}}{{{V_\infty }}} + {C_{j\dot \alpha }}\frac{{{\omega ^2}{z_m}}}{{{V_\infty }}}\left( {\frac{C}{{2{V_\infty }}}} \right)\sin (\omega t) $

(10)

所以：

$ {C_{j\dot \alpha }} = \left( {\int\limits_{{T_n}}^{{T_{n + 1}}} \Delta {C_j}\sin (\omega t){\rm{d}}t} \right)/\left( {\frac{{{\omega ^2}{z_m}{\rm{ \mathsf{ π} }}C}}{{2V_\infty ^2}}} \right) $

(11)

将式(8)和式(11)相减就可以得到单独的${{C}_{j\dot{\alpha }}}$和C_jq，根据图 1(a)和图 1(b)的结果，可以求得一阶气动力模型的各个导数值，升力系数对${\dot{\alpha }}$和q的导数值分别为-38.945 7和6.674 7，力矩系数对和q的导数值分别为-2.059 5和-1.334 4。

1.2 二阶气动模型

根据Etkin气动力模型^[2-3]，非定常气动力可以表示为：

$ \begin{array}{l} \Delta {C_j} = {C_{j\alpha }}\Delta \alpha + {C_{j\alpha }}\left( {\frac{D}{{2{V_\infty }}}} \right)\Delta \dot \alpha + {C_{j\ddot \alpha }}{\left( {\frac{D}{{2{V_\infty }}}} \right)^2}\Delta \ddot \alpha + \\ {\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} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{jq}}\left( {\frac{D}{{2{V_\infty }}}} \right)\Delta q + {C_{j\dot q}}{\left( {\frac{D}{{2{V_\infty }}}} \right)^2}\Delta \dot q \end{array} $

(12)

由于式中${{C}_{j\ddot{\alpha }}}$和${{C}_{j\dot{q}}}$的量纲相同，都是空气动力系数对角度随时间二阶变化率的导数，所以在本文中称式(12)为二阶气动力模型。与1.1节类似，由俯仰运动的气动力变化规律可以得到：

$ \left\{ {\begin{array}{*{20}{l}} {{C_{j\dot q}} = {C_{j\alpha }}/{k^2} - \left( {\int\limits_{{T_n}}^{{T_{n + 1}}} \Delta {C_j}\sin (\omega t){\rm{d}}t} \right)/\left( {{k^2}A\frac{{\rm{ \mathsf{ π} }}}{\omega }} \right)}\\ {{C_{j\dot \alpha }} + {C_{jq}} = \left( {\int\limits_{{T_n}}^{{T_{n + 1}}} \Delta {C_j}\cos (\omega t){\rm{d}}t} \right)/\left( {kA\frac{{\rm{ \mathsf{ π} }}}{\omega }} \right)} \end{array}} \right. $

(13)

由沉浮运动的气动力变化规律可以得到：

$ \left\{ {\begin{array}{*{20}{l}} {{C_{j\dot \alpha }} = \left( {\int\limits_{{T_n}}^{{T_{n + 1}}} \Delta {C_j}\sin (\omega t){\rm{d}}t} \right)/\left( {\frac{{{\omega ^2}{z_m}{\rm{ \mathsf{ π} }}C}}{{2{V_\infty }^2}}} \right)}\\ {{C_{j\ddot \alpha }} = \left( {\left( {\int\limits_{{T_n}}^{{T_{n + 1}}} \Delta {C_j}\cos (\omega t){\rm{d}}t} \right)/\left( {\frac{{{z_m}{\rm{ \mathsf{ π} }}}}{{{V_\infty }}}} \right) + {C_{j\alpha }}} \right)/{k^2}} \end{array}} \right. $

(14)

式(13)减去式(14)就可以得到单独的${{C}_{j\dot{\alpha }}}$和C_jq，根据图 1(a)和图 1(b)的结果，可以求得一阶和二阶气动力模型的各个气动导数值，一阶导数与上节完全相同，升力系数对${\ddot{\alpha }}$和${\dot{q}}$的导数值分别为578.511 8和-38.275 2，力矩系数对${\ddot{\alpha }}$和${\dot{q}}$的导数值分别为18.899 6和12.904 4。式(10)、(11)与式(12)~(14)相比，一阶气动导数完全一样，这也间接说明传统上忽略高阶导数的做法对低阶气动导数的求解结果并没有影响，这是传统上普遍采用Etkin气动力模型进行小迎角、小扰动飞行包线范围内动态稳定性分析的重要原因。

2 气动力建模结果比较

为了定量考察这些气动力模型对强迫运动过程非定常迟滞效应模拟的适用程度，本节对比了这些气动力模型的计算结果与直接采用CFD进行计算得到的结果。

图 2分别是采用一阶和二阶Etkin气动力模型计算得到的强迫俯仰运动、强迫沉浮运动以及耦合运动的气动力与采用CFD方法得到的气动力迟滞曲线的对比图。由图 2(a)可知，对于强迫俯仰运动，采用二阶气动导数得到的升力系数与CFD计算值完全重合，而采用一阶气动导数得到的升力系数误差随着迎角的增加而增大，在最大迎角位置比CFD计算值大50%。对于俯仰力矩系数，一阶模型的误差较大，而二阶模型的结果与CFD计算值虽然不像升力系数那样完全重合，但是吻合程度也较好。

表 1是俯仰运动不同位置处的不同变量对该时刻非定常气动力的贡献情况，需要注意的是，强迫俯仰运动时${\dot{\alpha }}$与q的数值相等，${\ddot{\alpha }}$与${\dot{q}}$的数值相等。由表 1可知，在俯仰运动1/4周期时，到达抬头最大位置处，此时${\ddot{\alpha }}$对应的非定常升力和力矩分别占总的非定常升力和力矩的-37.98%和-45.72%，${\dot{q}}$对应的非定常升力和力矩分别占总的非定常升力和力矩的2.51%和-31.21%。由于忽略了${\ddot{\alpha }}$与${\dot{q}}$的影响，所以一阶气动力模型对应的升力和力矩系数与CFD计算值差别较大，见图 2(a)。在俯仰运动3/4周期时，到达低头最大位置处。此时${\ddot{\alpha }}$与${\dot{q}}$这2项对非定常特性的贡献也比较大，所以与1/4周期时的情况类似，一阶模型的计算结果误差较大。而在1/2周期时，翼型低头经过初始位置，由于Δα、Δ${\ddot{\alpha }}$以及${\dot{q}}$的数值刚好为0，此时非定常气动力主要由${\dot{\alpha }}$和q项产生，所以一阶气动力模型就能较为准确地重现出非定常气动力。在一个周期时，翼型抬头经过初始位置，此时的非定常气动力贡献情况与1/2周期时相反，此时由于Δ${\dot{\alpha }}$和Δq的数值刚好为0，此时非定常气动力主要由α、${\ddot{\alpha }}$以及${\dot{q}}$项产生。这些分析结果与图 2(a)的结论一致，说明在俯仰运动过程中，二阶模型才能更准确地再现出非定常气动力和力矩。

图 2(b)是采用一阶和二阶Etkin气动力模型计算得到的强迫沉浮运动的气动力与采用CFD方法得到的气动力迟滞曲线的对比图。由图 2(b)可知，对于强迫沉浮运动，采用二阶气动导数得到的升力系数与CFD计算值几乎重合，而采用一阶气动导数得到的升力系数误差较大，在最大纵向位移位置比CFD计算值大90%。对于俯仰力矩系数，一阶模型的误差较大，而二阶模型的结果与CFD计算值虽然不像升力系数那样完全重合，但是吻合程度也较好。

表 2是沉浮运动不同位置处迎角的各阶导数对非定常气动力的贡献情况，由于强迫沉浮运动时没有俯仰角速度，所以在表 2中没有出现q和${\dot{q}}$项对应的非定常气动力。由表 2可知，在沉浮运动1/4周期时，到达纵向最大位置处，此时${\ddot{\alpha }}$对应的非定常升力和力矩分别占总的非定常升力和力矩的-94.84%和148.24%。由于忽略了${\ddot{\alpha }}$的影响，所以一阶气动力模型对应的升力和力矩系数与CFD计算值差别较大，见图 2(b)。在沉浮运动3/4周期时，到达纵向最小位置处。此时${\ddot{\alpha }}$对非定常特性的贡献也比较大，其数值分别为-36.92%和-49.93%，所以与1/4周期时的情况类似，一阶模型的计算结果误差较大。而在1/2周期时，翼型下沉经过初始位置，由于Δ${\dot{\alpha }}$的数值刚好为0，此时非定常气动力主要由α和${\ddot{\alpha }}$产生。在一个周期时，翼型回到初始位置，此时的非定常气动力贡献情况与1/2周期时刚好相反，由于Δα和Δ${\ddot{\alpha }}$的数值为0，所以此时非定常气动力主要由Δ${\dot{\alpha }}$产生。这与图 2(b)的结论一致，说明在沉浮运动过程中，二阶模型才能更准确地再现出非定常气动力和力矩。

图 2(c)是采用一阶和二阶Etkin气动力模型计算得到的强迫俯仰/沉浮耦合运动的气动力与采用CFD方法得到的气动力迟滞曲线的对比图。由图 2(c)可知，对于强迫耦合运动，采用二阶气动导数得到的升力系数与CFD计算值完全重合，而采用一阶气动导数得到的升力系数误差较大，在最大纵向位移位置比CFD计算值大145%。对于俯仰力矩系数，一阶模型的误差较大，而二阶模型的结果与CFD计算值虽然不像升力系数那样完全重合，但是吻合程度也较好。

表 3是耦合运动不同位置处各导数对非定常气动力的贡献情况。由表 3可知，在耦合运动1/4周期时，到达纵向最大位置处，此时也处于抬头的最大位置，此时${\ddot{\alpha }}$对应的非定常升力和力矩分别占总的非定常升力和力矩的-45.35%和-70.82%，${\dot{q}}$对应的非定常升力和力矩分别占总的非定常升力和力矩的2.19%和-35.25%。由于忽略了${\ddot{\alpha }}$与${\dot{q}}$这2项的影响，所以一阶气动力模型对应的升力和力矩系数与CFD计算值差别较大，见图 2(c)。在耦合运动3/4周期时，到达纵向最小位置处，此时也处于低头的最大位置。此时${\ddot{\alpha }}$对非定常特性的贡献也比较大，其数值分别为-245.66%和21.26%，${\dot{q}}$对非定常特性的贡献分别为25.87%和23.00%，所以与1/4周期时的情况类似，一阶模型的计算结果误差较大。而在1/2周期时，翼型下沉经过初始位置，由于Δ${\dot{q}}$的数值刚好为0，对非定常气动力没有贡献。在一个周期时，翼型回到初始位置，此时Δq的数值刚好为0，对非定常气动力没有贡献。虽然耦合运动气动力并不是单独运动的简单叠加，但是通过对经典Etkin气动力模型的二阶延拓，能准确地再现出沉浮/俯仰耦合运动过程的非定常气动特性。

3 结论

1) 不论是强迫俯仰运动、沉浮运动，还是俯仰/沉浮耦合运动，将气动导数拓展至迎角和俯仰角的二阶导数，都可以十分精确地重现出强迫运动过程中的非定常升力变化情况。

2) 由于俯仰力矩的迟滞曲线并不是简单的椭圆形，二阶模型计算出的强迫运动过程的俯仰力矩与CFD计算值的吻合程度不像升力那么好，说明俯仰力矩的模型要比升力更加复杂。

3) 俯仰/沉浮耦合运动的非定常气动力并不是俯仰运动和沉浮运动的简单叠加，说明精确的气动力建模还需要深入考虑其他变量的影响。

本文的研究结果表明，Etkin气动力模型对于非线性较强的气动力建模仍然具有较好的适用性，但是，对于三维流动以及接近失速迎角情况下的非定常气动力的建模，需要更加深入地讨论马赫数、减缩频率、更高阶导数以及交叉导数在非定常气动力模型中的作用。