NACA0012翼型跨声速强迫运动非定常气动力模型

2021-01-22张庆叶正寅

哈尔滨工程大学学报 2020年11期

张庆，叶正寅

(1.西安航空学院飞行器学院，陕西西安 710077； 2.南洋理工大学机械与航空工程学院，新加坡 639798； 3.西北工业大学航空学院，陕西西安 710072)

气动导数作为描述飞行器机动飞行和受扰动时气动特性的关键性气动参数，在飞行器气动性能、控制系统和总体设计中扮演着非常重要的作用[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]：

α(t)=α0+Asin(ωt)=α0+Asin(2πft)

(1)

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

本文定义减缩频率为：

(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)=z0+zmsin(ωt)

(3)

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

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

α(t)=α0-ωzmcos(ωt)/V∞

(4)

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

图1 不同运动过程升力和力矩系数随时间变化Fig.1 History of lift/moment coefficients in different motions

1.1 一阶气动模型

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

ΔCj=Cj-Cj0=CjαΔα+

(5)

(6)

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

(7)

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

(8)

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

(9)

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

(10)

所以：

(11)

1.2 二阶气动模型

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

(12)

(13)

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

(14)

2 气动力建模结果比较

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

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

图2 不同运动过程升力和力矩系数迟滞曲线Fig.2 Comparison of lift/moment coefficients indifferent motions

表1 俯仰运动不同位置非定常气动力分布情况

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

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