二维机翼副翼受PID控制偏转的数值模拟

邱福生，江 辉

(沈阳航空航天大学 航空宇航学院，沈阳 110136)

通过风洞试验获取机翼在非定常状态下的气动特性比较困难，而试飞试验所需要的成本比较高且风险大，但随着计算机的发展，可利用计算流体力学(computational fluid dynamics,CFD)实时获取机翼在非定常状态下受控运动的流场时序变化。对二维机翼襟翼、增升装置运动的非定常数值模拟在早期研究中已经展开[1-2]。近几年，对气动/运动/控制一体化耦合的研究开始增多[3-5]，但是都主要集中在飞行器受控机动行为的研究，对机翼副翼受控运动与流场数值模拟相结合的研究相对较少，因此对其进行深入探究有重要意义，其中所用模拟算法是基于开源软件包OpenFOAM[7]开发实现的。本文首先介绍运用动态重叠网格技术耦合PID控制的副翼偏转过程中的NS方程的数值方法，然后模拟二维NACA0012翼型定常状态的表面压力系数的分布和俯仰振荡运动过程的升力系数变化，并利用Fluent软件包模拟结果验证该数值算法的可靠性，最后对受PID控制的二维机翼副翼偏转过程进行数值模拟，初步建立机翼的流场计算与机翼副翼受PID控制运动相结合的直接数值模拟方法。

1.1 基于动态重叠网格技术的数值求解方法

为实现副翼的偏转动作，采用动态重叠网格技术[8]，其原理是把流场计算域划分成各个子计算域，然后各个子计算域通过动态装配与插值来实现流场信息的交换，最终计算出流场信息，主要涉及网格点的“寻点[9-10]”、“挖洞[11-13]”、“插值[14-15]”等关键技术：“挖洞”即去除不参与计算的网格单元；
“寻点”即寻找重叠网格之间实现信息交换的网格单元；
“插值”即建立贡献单元与插值单元的插值关系。

本文假设流场为湍流流动，控制方程采用NS方程，空间离散采用有限体积离散，选用两方程k-ω剪切应力运输(shear stress transfer,SST)模型并行计算加速数值求解过程。

1.2 PID控制算法

飞机在作机动动作时需要控制副翼达到某一偏转角度，从而需要一套控制算法来控制副翼的偏转，因此副翼偏转的控制模拟采用了增量式PID控制算法[16]。增量式PID控制算法具体如下

Δu(k)=Kp·Δe(k)+Ki·Δe(k)+Kd[Δe(k)-Δe(k-1)]

(1)

Δe(k)=e(k)-e(k-1)

(2)

其中:k表示采样序号；
Δu(k)表示第k次采样时刻的计算输出增量；
e(k)表示第k次采样时刻的输入偏差值；
e(k-1)表示第k-1次采样时刻的输入偏差值；
Kp表示比例系数；
Ki表示积分系数；
Kd表示微分系数。

1.3 气动/PID控制运动耦合计算

气动/PID控制运动耦合计算计算流程如图1所示。首先使用网格划分工具生成背景网格和二维带副翼机翼的重叠网格，然后根据初值条件计算初始定常流场。当流场场量收敛之后，气动与PID控制运动开始结合计算，达到设定时间时结束计算。当非定常NS方程数值计算时间点小于PID控制器的采样时间点时，非定常NS方程的数值计算继续向前推进，当时间推进满足采样时间时才开始采样，然后更新控制副翼的控制输出。副翼根据控制输出开始运动，利用重叠网格技术更新重叠网格与背景网格的关系，再进行非定常NS方程的数值计算，如此循环推进副翼受控运动的非定常流场时序变化过程。

图1 气动/PID控制运动耦合计算流程图

在非定常NS方程的数值计算中的每个时间步都需要先利用重叠网格技术进行重新装配、确定洞单元、插值单元、参与计算的单元，然后才开始求解计算。某时刻背景网格与网格的重叠关系如图2所示，重叠网格技术处理之后确立的单元类型如图3所示。

图2 背景网格与重叠网格的重叠关系

图3 单元类型

2.1 NACA0012翼型的定常模拟

为了验证重叠网格技术和数值计算算法的定常模拟可靠性，首先模拟了NACA0012翼型在初始迎角为10°、来流v=30 m/s时的定常工况。其中湍流模型采用k-ωSST模型，机翼与流场网格的划分和它们之间的重叠关系如图4所示，随后采用Fluent软件包在相同的条件下，对其进行对比模拟，NACA0012机翼表面压力系数分布如图5所示，Cp表示压力系数，X/L表示相对弦长，结果表明该数值算法计算结果与Fluent软件包模拟的结果基本吻合，因此该数值算法可用来进行机翼定常数值模拟。

图4 NAC0012机翼重叠网格(定常)

图5 NACA0012机翼表面压强分布

2.2 NACA0012翼型的非定常模拟

为了验证重叠网格技术和数值计算算法的非定常模拟可靠性，模拟了以翼型前缘点为旋转中心开始做俯仰振荡运动的非定常工况。模拟设置NACA0012翼型的初始俯仰角为0°，俯仰振荡运动的攻角变化规律为

(3)

NACA0012翼型的非定常模拟的初始条件为来流v=30 m/s，采用k-ωSST湍流模型，机翼与流场网格的划分和它们之间的重叠关系如图6所示。机翼在俯仰振荡过程中的升力系数随时间的变化如图7所示，Cl为升力系数，并利用Fluent软件包进行了对比模拟，该数值计算结果与Fluent计算结果比较吻合，因此该数值算法可用来进行机翼非定常数值模拟。

图6 NACA0012机翼重叠网格(非定常)

图7 升力系数随时间的变化

机翼在俯仰振荡过程中的压强云图如图8所示，机翼在顺时针变化过程中，可观测高压点逐渐从前缘点移至机翼的上表面；
机翼在逆时针变化过程中，可观测高压点逐渐从机翼的上表面移至前缘点，再移至机翼下表面。

图8 机翼在俯仰振荡过程中的压强云图

2.3 副翼受PID控制偏转模拟

选取了一个带副翼的二维机翼，其中机翼、副翼、流场网格的划分及它们之前的重叠关系如图9所示。由于副翼尺寸比较小，且需要满足副翼网格与背景网格插值匹配，因此对副翼网格所在背景网格进行加密。副翼受控偏转过程中，简化PID控制模型，只对偏转角度进行控制，设置了副翼的偏转角速度每秒不超过30°，PID控制器在采样时间内增量大于副翼在采样时间内能够偏转的最大角度时就不对其进行控制。副翼保持最大角速度偏转，当增量小于副翼在采样时间内能够偏转的最大角度，则受PID控制器控制，最终结果反映的是副翼总的偏转角度，控制机翼的副翼偏转的PID控制器的采样时间周期0.02s,控制器输入的信号是期望副翼偏转角度为11°，反馈信号是副翼总的偏转角度，在实际运用中需要角位移传感器来获取其偏转角度。机翼在受控运动过程中的压强云图如图10所示。

图9 机翼带副翼重叠网格

图10 机翼在受控运动过程中的压强云图

机翼在PID控制下的偏转角度随时间的变化如图11所示。从图11中可看出，PID参数为Kp=0.38、Ki=0.72、Kd=0.007时，副翼在0.35 s后开始接近期望偏转角度，再逐渐衰减振荡，最后接近期望偏转角度。PID参数为Kp=0.38、Ki=0.006、Ki=0.007时，副翼在0.35 s后开始接近期望偏转角度，副翼偏转角度与期望偏转角度存在稳态误差，与前一组PID参数比较可得，增大值可加速收敛，减少调节时间。PID参数为Kp=0.79、Ki=0.72、Kd=0.007时，副翼在0.35 s后开始接近期望偏转角度，在接近期望偏转角度之后，副翼偏转存在振荡现象，与第一组PID参数比较可得，减小Kp值可减小副翼转角在接近期望偏转角度之后的振荡。

图11 机翼受控偏转角度随时间变化图

针对机翼副翼在PID控制下的偏转问题，基于动态重叠网格技术的非定常计算方法，初步建立了机翼的流场计算与其副翼受PID控制运动相结合的直接数值模拟方法，并模拟了带副翼的二维机翼的流场特性，得出如下结论：

(1)合适的PID控制参数可提高副翼偏转的稳定性，并使偏转到期望角度的收敛过程加速。

(2)依据现有计算设备和数值计算方法能够实现气动与PID控制运动相结合的直接数值模拟。

(3)该数值方法能够为控制算法的设计提供参考，并能够反映副翼受控运动的非定常流场时序变化。