飞行器双自由度耦合振荡的数值模拟方法

肖中云，刘刚，江雄，王建涛

中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

飞行器双自由度耦合振荡的数值模拟方法

肖中云*，刘刚，江雄，王建涛

中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

针对模拟双自由度(2DOF)运动的“PQR”强迫振荡试验装置，通过坐标变换，将机构运动转化为空间欧拉角进行表示，实现对机构俯仰/偏航、俯仰/滚转和偏航/滚转等3种耦合运动的统一描述，采用预定运动轨迹的动网格计算技术，构造了基于非定常雷诺平均Navier-Stokes (URANS)方程的飞行器双自由度振荡数值模拟方法。围绕某复杂构型飞行器低速流动，在模型迎角10°、滚转角±40°、偏航角±40°振荡预设条件下分别对单自由度偏航、滚转振荡和双自由度偏航/滚转振荡进行了模拟，得到的气动系数迟滞回线在形态和量值上与风洞试验结果吻合，证明当前方法可以作为飞行器复杂耦合运动的一种有效研究手段。针对偏航/滚转耦合振荡，从运动形式和气动特性上与单独偏航、滚转运动进行了对比分析，结果表明，耦合运动在气流角、绕体轴的旋转角速度等方面与单自由度叠加效果不同，力矩系数迟滞回线存在曲线交叉的现象，表现出与单自由度振荡不同的阻尼特性。在当前模拟状态下，偏航/滚转振荡流场以横向分离涡运动为主要特征，可以推测大迎角情况下纵横向分离涡流动结构更加复杂，耦合作用更强，需要进一步在分离流模拟方法上开展研究。

双自由度；大幅振荡；耦合运动；非定常；数值模拟

飞行器过失速机动是一种大迎角下的多自由度耦合运动，气动力呈现出严重的非对称、非线性和非定常特性，由各种实验数据线性叠加形成总的气动系数的方法，理论上缺少依据，实际上会存在很大偏差。围绕过失速机动问题，世界各航空发达国家都开发了相应的风洞试验与理论分析方法[1-5]，目的是为了获取飞行器在大迎角过失速区域的非线性、非定常气动力特性，有助于更好理解动态空气动力的复杂现象，以及帮助发展更好的飞行力学数学模型。

随着对飞行器过失速机动能力的日益重视，国内外发展了大量的风洞动态试验技术[6-8]，包括大迎角动导数、旋转天平、旋转振荡、大幅振荡试验技术等。“PQR”动态试验系统是一种典型的大振幅双自由度试验装置，可进行大幅俯仰、偏航、滚转振荡及其双自由度耦合振荡、大幅俯仰拉升及俯仰和滚转耦合拉升试验[9]。大幅振荡试验是进行飞行机动模拟，特别是过失速机动模拟的基本试验数据。目前单自由度振荡(俯仰、偏航和滚转)已经开展了较多的研究，在获取动导数、动态非定常时间迟滞效应上十分有效，同时可以用于模拟简单的机动飞行动作，如眼镜蛇机动和机翼摇滚等。对于复杂机动飞行动作来说，运动通常为多个自由度的耦合运动，如Herbst机动，飞行器在大迎角下同时包含了偏航和滚转的耦合运动。这就要求不断提高试验模拟的运动相似性，从单自由度向双自由度耦合、多自由度机动的方向发展。

双自由度振荡是一种复合运动，运动参数包括平均角度、振幅、振荡频率、相位等，各参数的组合情况非常多，运动规律及气动特性十分复杂。文献[10]较早开展了三角翼的俯仰/滚转耦合运动试验研究，后续又开展了典型飞机的偏航/滚转耦合气动力试验[11-12]，模型支撑迎角达到60°，将耦合运动时的气动特性与单自由度运动叠加得到的气动特性进行了对比。文献[13]介绍了生产型低速风洞的双自由度大幅振荡试验技术，在该技术支撑下开展了某复杂构型飞机偏航/滚转耦合运动非定常气动力特性试验研究，获得了重复性好、可靠性高的试验数据[14]。以上这些试验为研究飞行器大迎角气动特性积累了宝贵经验，同时也为非定常气动建模的发展提供了有力支撑[15-17]。

目前，用CFD方法模拟飞行器多自由耦合运动的工作还比较少见。文献[18]对二维平板机翼前冲、俯仰振荡以及两者的合成运动进行了模拟，研究大尺度涡结构对气动的影响。文献[19]研究了低雷诺数下的机翼俯仰和下沉运动，以此模拟微型飞行器的非定常空气动力效应。随着计算机运算速度的提高，非定常计算可采用的计算规模大幅提高，对分离涡流动、复杂湍流的模拟能力日益增强。发展多自由度大幅振荡的数值模拟方法，有助于摸清复杂气动特性背后的流动机理，对风洞试验起到相互验证与校核的作用。本文基于描述姿态运动的欧拉角方法和动网格数值计算，发展了模拟飞行器双自由度大幅振荡的数值模拟方法，探讨了耦合运动的特点及其与解耦运动的关系，最后以偏航/滚转振荡为例，通过对某复杂外形飞行器的模拟验证了本文方法。

1 “PQR”机构的姿态运动描述方法

图1 俯仰/滚转耦合运动Fig.1 Pitch/roll coupled motion

图2 偏航/滚转耦合运动Fig.2 Yaw/roll coupled motion

图3 俯仰/偏航耦合运动Fig.3 Pitch/yaw coupled motion

Aα=I3×3cosα+(1-cosα)e·eT-e×sinα

(1)

式中：I3×3为单位矩阵；e×定义为

将式(1)展开，可得

(2)

绕体轴旋转的方向余弦矩阵为

(3)

式中：Cα代表cosα,Sα代表sinα,其余类似，两自由度动态试验包括了俯仰/滚转、俯仰/偏航、偏航/滚转3种类型，为了统一输入输出接口，将两自由度运动定义为“绕体轴旋转+绕地面系下任意轴旋转”两次连续旋转过程。用Xd和Xt表示地轴系和体轴系下的坐标，Aφ,θ,ψ为绕体轴旋转的方向余弦矩阵，得到坐标变换关系为

等价于：

Xt=Aφ,θ,ψAαXd

(4)

式(4)说明了双自由度运动的坐标变换包含了两次矩阵左乘的运算，注意此时运算结果与乘法顺序相关，矩阵的顺序不能互换。最后得到的坐标变换矩阵为

A=Aφ,θ,ψAα

(5)

下面考虑方向余弦矩阵向欧拉角的转换，如果考虑俯仰角在-90°～+90°之间取值，滚转角和偏航角在-180°～+180°之间取值，则方向余弦矩阵到欧拉角的转换公式为

(6)

θ=arcsin(-A13)

(7)

(8)

(9)

ifu≠0 then

{α=arctan(w/u)

β=arctan(-(vcosα)/u)

elseifw≠0 then

{α=(π/2)sign(w)

β=arctan(-v/w)

elseifw=0 then

{α=π/2

β=(π/2)sign(-v)

2 数值计算方法

流动控制方程采用任意拉格朗日-欧拉(Arbitrary Lagrangian Eulerian，ALE)形式下的可压缩非定常雷诺平均Navier-Stokes(URANS)方程，允许网格的任意运动和变形，在绝对坐标系下的积分形式为

(10)

式中：Q为守恒变量;Ω表示控制体的体积；S为控制体边界面的面积;n为单位法向矢量；VΩ为边界面上的网格运动速度；HI和HV分别为对流通量和黏性通量。非定常计算采用双时间步方法，其中内层时间离散采用了点隐式LU-SGS方法，真实时间采用了三层二阶精度格式离散，在每个真实时间步内通过内层时间迭代使气动力积分值收敛。黏性离散采用全Navier-Stokes方程方法，湍流采用的是两方程k-ω剪切应力输运(SST)模型进行模拟。由于是计算低速大迎角问题，为降低低速流动的刚性，流动控制方程采用了低速预处理方法，同时为加速流场收敛还采用了3重V循环的多重网格方法，每一重网格的迭代步数为5步，每个真实时间步内多重网格循环的最大次数为200步。在当前参数设置下，基本上都能保证子迭代过程的气动力收敛，如果在迭代过程中密度残差下降到两个量级以上，则被认为是子迭代收敛，提前退出循环。

计算模型如图4所示，平均气动弦长(参考长度)为0.426 9 m，来流风速为25 m/s，雷诺数Re=1.71×106，第一层网格距离物面5×10-6m，满足湍流附面层模拟y+=Ο(1)的要求，网格采用多块对接网格(包括232块)，单元总数约为1 228万。

图4 计算模型和滚转振荡瞬时流场Fig.4 Calculation model and instantaneous flow field of roll oscillation

在开展动态计算之前，首先通过静态流场对计算方法进行校核，模拟迎角α=0°～80°。图5给出的是计算与试验的升力系数CL比较，可以看到，在所有迎角范围内计算与试验值吻合良好，最大升力系数约为1.48，对应的迎角约为35°。图6给出的是阻力系数CD曲线比较，可以看到，在迎角α=40°以前，计算与试验值吻合较好，在α>40°以后，计算值较试验值偏高，最大偏差量约为10%。图7给出了俯仰力矩系数Cm的比较，也可以看到两者之间吻合较好，并且可以看到该飞行器模型在中小迎角范围内俯仰方向是静不稳定的，俯仰力矩随迎角增大而增大，而在大迎角下(α>30°)是俯仰静稳定的。总的来说，当前计算在规律和趋势上和试验值都比较一致，反映了低速大迎角流场的主要特征。由于目前公认的URANS方法模拟大分离流动存在不足，所以认为在失速迎角以后计算与试验值存在一定偏差是合理的，在当前计算中主要表现在阻力值的预测上。

图5 升力系数比较Fig.5 Comparison of lift coefficients

图6 阻力系数比较Fig.6 Comparison of drag coefficients

图7 俯仰力矩系数比较Fig.7 Comparison of pitching moment coefficients

3 计算结果与分析

3.1 偏航、滚转单自由度大幅振荡

在计算偏航/滚转双自由度振荡之前，首先对偏航振荡和滚转振荡分别进行计算，单自由度滚转振荡的运动规律用正弦函数表示，φ=40°·sin(2πft),振荡频率f=0.8 Hz。图4给出的是计算模型和滚转振荡瞬时流场，对应瞬时滚转角为34.6°，分别给出了横向截面的马赫数Ma云图、空间流线与涡量等值面，可以看到在当前位置下，飞行器左右侧旋涡的位置和强度都发生改变，模拟得到的物面边界层与空间旋涡的形态分布合理。

图8给出了滚转和偏航力矩系数(Cl,Cn)的计算值与试验值比较，可以看到计算与试验值吻合很好，进一步验证了当前计算方法的合理性。尽管计算流体力学存在计算坐标系、飞机体轴系之分，这里为叙述方便将角度定义的方向和气动力矩的方向都统一到机体坐标系下。在转动方向与力矩方向定义一致的前提下，可以根据力矩曲线的斜率或方向来判断飞机的力矩特性[10]。具体来说，静态情况下如果力矩系数的斜率为负，说明在该方向上的力矩是静稳定的，否则为静不稳定；在动态情况下，如果迟滞回线沿逆时针方向旋转，则说明在该方向上的运动对力矩起阻尼作用，如果迟滞回线沿顺时针方向旋转，则说明在该方向上的运动对力矩起发散作用。

图8 滚转振荡的滚转与偏航力矩系数(α=10°)Fig.8 Roll and yaw moment coefficients in roll oscillation (α=10°)

如图8所示，对于滚转力矩来说，从最小滚转角对应的力矩点和最大滚转角对应的力矩点之间连线，可以得出滚转力矩呈左上右下的走势，即正的滚转角对应的滚转力矩为负，说明飞机在滚转方向上是静稳定的。但是对当前滚转大幅振荡来说，滚转力矩形成了很大的迟滞圈，即滚转力矩产生了很强的动态非定常效应。图中标注了迟滞回线的旋转方向，可以看到滚转力矩曲线呈逆时针走向，并且滚转力矩环呈O型，这说明在整个滚转角历程中滚转运动对滚转力矩的变化起阻尼作用。图中同时还给出了偏航力矩的迟滞回线，可以看到当前偏航力矩并非小量，说明飞行器滚转和偏航之间存在较强的耦合作用，同时偏航力矩迟滞回线的旋转为逆时针方向，意味着滚转引起的偏航交叉动导数起阻尼作用。在振荡过程中，当飞行器处于平衡位置时，图8中滚转力矩迟滞圈包围面积很大，对比上行和下行状态，当飞行器处于平衡位置时，两者的姿态角和气流迎角、侧滑角都一致，角速度不同、时间历程不同，造成了滚转力矩系数的较大差异。

图9 偏航振荡的滚转与偏航力矩系数 (α=10°)Fig.9 Roll and yaw moment coefficients in yaw oscillation (α=10°)

单自由度偏航振荡的运动规律用正弦函数表示，ψ=40°sin(2πft),振荡频率f=0.8 Hz。图9给出的是偏航振荡的计算与试验值比较，包含滚转和偏航力矩系数，可以看到计算与试验在趋势与量值上都比较一致，可以说计算正确反映了偏航运动下的横航向力矩特性变化，不足之处是在较大偏航角下，计算的偏航力矩与滚转力矩峰值较试验值偏小，其余地方两者都吻合较好。对于偏航力矩来说，在偏航角ψ=-20°～20°区间，偏航力矩呈左高右低的走向，说明此时偏航力矩与运动方向相反，有利于飞行器的稳定；但是在ψ<-20°和ψ>20°的区间，偏航力矩的走向与运动方向相同，不利于偏航方向的稳定。由于偏航力矩主要是飞行器垂尾产生，在较大的偏航角下，垂尾效率下降甚至失速，会导致偏航稳定的效率下降。从总的曲线变化趋势上看，偏航力矩迟滞环呈逆时针变化，动导数起阻尼作用。另外，由偏航引起的滚转力矩迟滞环呈顺时针变化，偏航引起的滚转交叉动导数起发散作用。

3.2 偏航/滚转耦合大幅振荡

偏航/滚转耦合的运动机构见图2，机构偏角定义为

(11)

图10 偏航/滚转运动中的欧拉角、迎角和侧滑角变化 Fig.10 Euler angle，angle of attack and sideslip angle variations in yaw/roll coupled motion

设定滚转和偏航的振荡幅值φA=ψA=40°,振荡频率均为f=0.8 Hz,初始相位相同，即飞行器从平衡位置开始向右偏航、向右滚转。图10给出的是偏航/滚转耦合运动下的欧拉角、迎角和侧滑角变化曲线，可以看到，在运动过程中俯仰角θ始终为0°，偏航角ψ和滚转角φ在±40°之间变化，偏航与滚转方向的欧拉角与机构运动偏角相等。图中同时给出了迎角、侧滑角的变化，其中平衡位置的来流迎角α=10°，可以看到飞行器在姿态角发生变化的同时，气流角也在发生变化，其中迎角在10°～36°之间变化，侧滑角则在±22°之间变化。

图11给出的是体轴系的3个角速度变化曲线，可以看到，除了偏航和滚转方向的角速度外，俯仰方向上也产生了角速度，并且俯仰角速度的周期是其他两个方向角速度周期的一半。由此可以看到，偏航/滚转耦合运动并不是偏航运动和滚转运动的简单叠加，其中包括了多个动态参数的变化，根据线性理论假设，将气动力/力矩系数通常表示为飞行状态参数的函数[21]，即

Cj0(α,β,V∞,h,δi)+Cjd(α,β,V∞,h,δi,p)+

Cjd(α,β,V∞,h,δi,q)+…

(12)

图12给出的是偏航/滚转振荡的升力系数CL、阻力系数CD和侧向力系数Cy比较，从图中可以看到，升力系数和阻力系数呈“∞”形变化，侧向力系数呈双“∞”形变化。计算值与试验值总体上吻合较好，一定程度说明本文算法的正确性，不足之处是最大振荡角度下计算的升力系数值略小于试验值，同时也导致了该角度下的阻力系数较试验值小。

图11 偏航/滚转运动中的体轴系角速度变化Fig.11 Angular velocity variations of body axial system in yaw/roll coupled motion

滚转力矩和偏航力矩变化是偏航/滚转耦合振荡重点关注的对象。图13给出了滚转力矩系数随滚转角变化曲线。可以看到当前力矩曲线呈两个“∞”形分布，在滚转角-18°<φ<18° 区域，曲线为逆时针旋转，说明滚转动导数起阻尼作用，而在φ<-18°和φ>18°的两个区域，曲线为顺时针旋转，说明此时滚转动导数起发散作用。

图14给出了耦合振荡下的偏航力矩系数变化曲线，可以看到当前曲线的旋转方向为逆时针，说明耦合振荡引起的偏航动导数起阻尼作用。图13和图14分别给出了滚转和偏航力矩系数的风洞试验值，从计算值与试验值的比较来看，两者在曲线形态和量值上都吻合较好，说明当前采用的计算方法是合理的。

图12 偏航/滚转振荡的升力、阻力和侧向力系数比较Fig.12 Comparison of lift, drag and side force coefficients in yaw/roll oscillation

图13 偏航/滚转振荡的滚转力矩系数Fig.13 Roll moment coefficients in yaw/roll oscillation

图14 偏航/滚转振荡的偏航力矩系数Fig.14 Yaw moment coefficients in yaw/roll oscillation

图15 耦合运动与单独运动叠加的偏航力矩系数比较Fig.15 Comparison of yaw moment coefficients between couple motion and single freedom motion addition

下面考虑耦合运动与单独运动叠加结果的比较，图15给出的是偏航力矩系数曲线，可以看到耦合结果与叠加结果具有较高的相似性，并且与单独偏航运动的结果相似，可见该飞行器的偏航力矩主要由偏航运动产生，相比之下滚转振荡产生的偏航力矩较小。图16给出的是滚转力矩系数曲线，可以看到，当前耦合运动的力矩系数在大的趋势上与图9偏航振荡比较相近，呈左下右上的分布形态。经过简单叠加得到的力矩曲线在较大滚转角范围内(-32°<φ<32°)都是沿逆时针方向旋转，在φ>32°和φ<-32°之后曲线出现了交叉，但包围的面积很小。出现交叉的原因是滚转运动和偏航运动的力矩曲线旋转方向相反，具体来说，滚转运动产生的滚转力矩系数曲线沿逆时针方向旋转，偏航运动产生的滚转力矩系数沿顺时针方向旋转。与单独运动叠加对比，耦合运动在小滚转角下表现出了类似的形态，力矩曲线沿逆时针方向旋转，不同的是当滚转角较大时(φ>18°和φ<-18°)，力矩曲线出现了交叉，旋转方向变为顺时针，并且每个顺时针旋转曲线包围的面积较大，与逆时针旋转所包围的面积相当。

另外，由于耦合运动的滚转力矩曲线形态发生改变，当飞行器从平衡位置向右侧滚转时，随滚转角增大时，滚转力矩加速上扬，在φ=34°附近达到最大值，该值较单独运动和单独运动叠加对应的滚转力矩要高出很多，力矩的作用方向是进一步加速飞行器的滚转运动。

图16 耦合运动与单独运动叠加的滚转力矩系数比较Fig.16 Comparison of roll moment coefficients between couple motion and single freedom motion addition

图17 耦合运动与单独运动叠加的侧滑角比较Fig.17 Comparison of side slip angles between couple motion and single freedom motion addition

为了弄清力矩曲线变化的原因，首先对耦合运动的姿态角与气流角变化进行了分析，从图10可以看到，耦合运动除了偏航和滚转方向的姿态角变化以外，可能对滚转力矩产生较大影响的是侧滑角的变化。为进一步分析，图17对单独运动和耦合运动的侧滑角进行了比较，可以看见，单独偏航运动的侧滑角在[-40°,40°] 区间变化，单独滚转运动的侧滑角在[-6.4°,6.4°] 区间变化，但方向相反，耦合运动的侧滑角约在[-22°,22°] 区间变化。从图9知道，偏航运动当偏航角在25°之前产生的滚转力矩单调上升，25°以后开始下降，这其中起主导作用的就是侧滑角的变化。而对于耦合运动来说，侧滑角始终在±25°内变化，该角度范围不足以使垂尾发生失速，使滚转力矩得以持续上升，这也就解释了图16耦合振荡的滚转力矩系数峰值增大的原因。从本文分析可以看到，非定常气动系数的影响因素很多，包括模型姿态角、气流角、对应角速度等，非定常运动的时间历程效应也是影响因素之一，这些因素的叠加导致了耦合运动气动特性的复杂性，同时给气动系数建模带来了难度。

4 结 论

采用URANS方法模拟预定运动轨迹,建立了专门针对飞行器双自由度大振幅振荡的数值模拟方法，完成了某复杂外形飞行器静态、偏航振荡、滚转振荡和偏航/滚转耦合振荡的模拟，结合试验数据对力矩特性的分布规律进行了分析。

1) “PQR”试验装置在进行双自由度试验时，模型方位与旋转顺序无关，在用CFD进行预定轨迹计算时，将用户输入用一次绕任意轴旋转叠加一次绕体轴旋转表示，输出用坐标变换的欧拉角及绕体轴的角速度表示，可以实现对任意俯仰/偏航、俯仰/滚转、偏航/滚转等耦合运动的描述。

2) 当前计算采用了低速预处理、多重网格等加速收敛技术，保证非定常计算每个子迭代过程气动力的充分收敛，得到的计算结果合理，在模型平衡位置迎角10°，滚转±40°偏角、偏航±40°偏角下的结果与风洞试验值吻合一致。

3) 飞行器双自由度运动不是两个单自由度运动的简单叠加，气动特性具有单自由度运动的部分特点，但迟滞回线的形态更加复杂，甚至表现出不同的运动阻尼特性。

4) 本文偏航/滚转耦合振荡的平衡位置迎角为10°，流场以横向分离涡运动为主要特征，可以推测大迎角情况下纵横向分离涡流动结构更加复杂，耦合作用更强，需要进一步在大分离流模拟方法上开展研究。

本文模拟外形为1∶16风洞缩尺模型，来流风速为25 m/s，基于平均气动弦长的雷诺数为百万量级，当前条件下本文URANS方法取得了与试验值比较吻合的结果，真实飞行条件下飞行器速度更快、雷诺数更高，对计算方法构成挑战。下一步有必要开展真实飞行条件下的计算研究，依托高雷诺数风洞试验来验证CFD方法；同时分离涡流动模拟是当前的热门研究方向，开展这类方法研究将有助于解决大迎角大分离流动的动态模拟问题。

致 谢

本研究得到了汪清研究员、刘志涛副研究员的帮助和指导，在此表示感谢。

[1] NELSON R C, PELLETIER A. The unsteady aerodynamics of slender wings and aircraft undergoing large amplitude maneuvers[J]. Progress in Aerospace Sciences, 2003, 39(2): 185-248.

[2] KRAMER B R. Experimental evaluation of superposition techniques applied to dynamic aerodynamics: AIAA-2002-0700[R]. Reston, VA: AIAA, 2002.

[3] LAN C E, LI J, YAU W, et al. Longitudinal and lateral-directional coupling effects on nonlinear unsteady aeodynamic modeling from flight data: AIAA-2002-4804[R]. Reston, VA： AIAA, 2002.

[4] WANG Z, LAN C E, BRANDON J M. Fuzzy logic modeling of lateral-directional unsteady aerodynamics：AIAA-1999-4012[R]. Reston, VA： AIAA, 1999.

[5] JENKINS J E. Nonlinear aerodynamic characteristics of a 65 degree delta wing in rolling motion: Implications for testing and flight mechanics analysis：AIAA-1997-0742[R].Reston, VA： AIAA, 1997.

[6] OWENS D B, BRANDON J M, CROOM M A, et al. Overview of dynamic test techniques for flight dynamics research at NASA LaRC: AIAA-2006-3146[R]. Reston, VA： AIAA, 2006.

[7] GILLARD W J. AFRL F-22 dynamic wind tunnel test results: AIAA-1999-04015[R]. Reston, VA： AIAA, 1999.

[8] 李伟. 飞机大振幅动态实验及其非定常气动特性研究[D]. 南京: 南京航空航天大学, 2008.

LI W. Unsteady aerodynamic experiments and characteristics for the aircraft oscillation in large amplitude[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2008 (in Chinese).

[9] 李周复. 风洞特种试验技术[M]. 北京: 航空工业出版社, 2010.

LI Z F. Special test technology in wind tunnel[M]. Beijing: Aviation Industry Press, 2010 (in Chinese).

[10] 黄达, 吴根兴. 三角翼俯仰滚转耦合运动气动特性研究[J]. 航空学报, 1999, 20(6): 485-488.

HUANG D, WU G X. Investigation of unsteady aerodynamic characteristics for a delta wing oscillating in large amplitude pitching roll motion[J]. Acta Aeronautica et Astronautica Sinica, 1999, 20(6): 485-488 (in Chinese).

[11] 黄达, 吴根兴. 飞机偏航-滚转耦合运动非定常空气动力实验[J]. 南京航空航天大学学报, 2005, 37(4): 408-411.

HUANG D, WU G X. Experiment on fighter oscillating in large amplitude yaw-roll motion[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2005, 37(4): 408-411 (in chinese) .

[12] HUANG D, WU G X. Unsteady rolling moment characteristics for fighter oscillating with yawing-rolling coupled motion[J]. Journal of Aircraft, 2006, 43(5): 1570-1573.

[13] 刘春明, 赵志军, 卜忱. 低速风洞双自由度大幅振荡试验技术[J]. 航空学报, 2016, 37(8): 2417-2425.

LIU C M, ZHAO Z J, BU C， et al. Double degree-of-freedom large amplitude oscillation test technology in low speed wind tunnel[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2417-2425 (in Chinese).

[14] 杨文, 卜忱, 眭建军. 某复杂构型飞机偏航-滚转耦合运动非定常气动力特性实验研究[J]. 实验流体力学, 2016, 30(3): 61-65.

YANG W, BU C, GUI J J. Investigation of the unsteady aerodynamic characteristics of a fighter with complex configuration undergoing yaw-roll coupling oscillation motion[J]. Journal of Experiments in Fluid Mechanics, 2016, 30(3): 61-65 (in Chinese).

[15] 孙海生, 张海酉, 刘志涛. 大迎角非定常气动力建模方法研究[J]. 空气动力学学报, 2011, 29(6): 733-737.

SUN H S, ZHANG H Y, LIU Z T. Comparative evaluation of unsteady aerodynamics modeling approaches at high angle of attack[J]. Acta Aerodynamica Sinica, 2011, 29(6): 733-737 (in Chinese).

[16] 汪清, 钱炜祺, 丁娣. 飞机大迎角非定常气动力建模研究进展[J]. 航空学报, 2016, 37(8): 2332-2347.

WANG Q, QIAN W Q, DING D. A review of unsteady aeeodynamic modeling of aircrafts at high angles of attack[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2332-2347 (in Chinese).

[17] 史志伟, 黄达, 吴根兴, 等. 耦合运动非定常气动模型对飞机飞行特性仿真的影响[J]. 航空学报, 2008, 29(6): 1424-1428.

SHI Z W, HUANG D, WU G X, et al. Effects of coupled motion unsteady aerodynamic model on flight characteristics simulation of aircraft[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(6): 1424-1428 (in Chinese).

[18] JANTZEN R, TAIRA K. Aerodynamic force modeling for unsteady wing maneuvers: AIAA-2014-0223[R]. Reston, VA: AIAA, 2014.

[19] MICHAEL V L. The high frequency, high-amplitude pitch problem: Airfoils, plates and wings: AIAA-2009-3686[R]. Reston, VA: AIAA, 2009.

[20] 张帆, 曹喜滨, 邹经湘. 一种新的全角度四元数与欧拉角的转换算法[J]. 南京理工大学学报, 2002, 26(4): 376-380.

ZHANG F, CAO X B, ZOU J X. A new large-scale transformation algorithm of quaternion to Euler angle[J]. Journal of Nanjing University of Science and Technology, 2002, 26(4): 376-380 (in Chinese).

[21] WANG F J, CHEN L. Numerical prediction of stability derivatives for complex configurations[J]. Procedia Engineering, 2015, 99(1): 1561-1575.

Numericalsimulationmethodsfor2DOFcoupledoscillationsofaircraft

XIAOZhongyun*,LIUGang,JIANGXiong,WANGJiantao

ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

Focusingonthe“PQR”forcedoscillationtestapparatusformodelingTwo-Degree-of-Freedom(2DOF)motions,coordinatetransformationisperformedtoobtainunifiedEulerangledescriptionsofthepitch/yaw,pitch/rollandyaw/rollcoupledmotions.BasedonmovinggridsforprescribedmotionandUnsteadyReynoldsAveragedNavier-Stokes(URANS)equations,anumericalalgorithmisdevelopedtosimulatecoupledoscillationofaircrafts.Surroundinglowspeedflowsofacomplexaircraft,single-degreeoscillationsofyawandrollandyaw/rollcoupledoscillationarerespectivelycomputedatangleofattackof10°andamplitudeofoscillationof40°.Thenumericalresultsagreewellwiththeexperimentaldata,showingthatcurrentalgorithmhaspromisingpotentialsinsimulationofcomplicatedcoupledmotionsofaircrafts.Asformovingprocessandaerodynamics,theyaw/rollcoupledmotioniscomparedwiththesuperpositionofsingleyawandrolloscillations,whichshowsthatdifferencesexistinbodyangularvelocitiesandflowangles,andthemomenthysteresiscurveoftheyaw/rollcoupledmotionshowscrossings,indicatingtotallydifferentdampingcharacteristicsfromsingleDOFoscillation.Incurrentyawandrolloscillations,flowsarecharacterizedbylaterallyseparatedvortexes.Itisassumedthatathighangleofattack,whenflowsaremorecomplicatedwithlateralandlongitudinalvortexesandstrongcoupledeffects,morestudiesshouldbedoneonnumericalmethodsofhighlyseparatedflows.

2DOF;largeamplitudeoscillation;coupledmotion;unsteady;numericalsimulation

2017-04-21；Revised2017-05-16；Accepted2017-06-07；Publishedonline2017-06-161118

URL：http://hkxb.buaa.edu.cn/CN/html/20171111.html

s：NationalBasicResearchProgramofChina(613221)；NationalNaturalScienceFoundationofChina(11572341);NationalKeyResearchandDevelopmentProgram(2016YFB0200704)

.E-mailscxiaozy@sina.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121345

V211.3

A

1000-6893(2017)11-121345-10

2017-04-21；退修日期2017-05-16；录用日期2017-06-07；< class="emphasis_bold">网络出版时间

时间：2017-06-161118

http://hkxb.buaa.edu.cn/CN/html/20171111.html

国家“973”计划(613221)； 国家自然科学基金(11572341)； 国家重点研发计划(2016YFB0200704)

.E-mailscxiaozy@sina.cn

肖中云，刘刚，江雄，等．飞行器双自由度耦合振荡的数值模拟方法J． 航空学报,2017,38(11)：121345.XIAOZY,LIUG,JIANGX,etal.Numericalsimulationmethodsfor2DOFcoupledoscillationsofaircraftJ.ActaAeronauticaetAstronauticaSinica,2017,38(11):121345.

(责任编辑：李明敏)