超高速飞行器平尾大迎角气动弹性特性研究

超高速飞行器由于飞行速度快、突防能力强而具有巨大的军事价值，成为实现“全球快速打击(PGS）”的利器。近年来各航空航天大国将其列为发展重点，进行了广泛研究［1-2］，其中较具代表性的是美国SR-72无人机计划［3］。

全动平尾具有操纵效率高，结构紧凑等优点，常作为超高速飞行器的气动控制面，靠翼轴末端作动器驱动起到稳定、配平和控制作用。同时，为减小阻力，构型多采用小展弦比切尖三角翼；翼型多选用超薄翼型，如菱形或双楔形。

另外，超高速飞行器的任务和设计需求要求全动平尾采用轻质材料，如C/SiC材料、镍合金材料和C/C复合材料等。平尾结构刚度较小、固有频率较低，造成弹性效应更为明显，使得全动平尾的气动弹性问题相对于机身而言更突出。超高速飞行器的操纵舵面特性对于平台设计和飞行安全非常重要，因此对于全动尾翼的气动弹性研究十分必要。

超高速飞行器气动弹性问题的研究从20世纪50年代开始兴起，由于计算条件等诸多限制，此时的研究主要集中于风洞试验以及超声速非定常气动力的工程算法。Lighthill［4］采用活塞理论求解二维翼型在超声速下的气动力分布，并与当时已有的其他方法对比，对结果进行了验证。Ashley和Zartarian［5］将活塞理论应用于翼型的气动弹性计算，分析了超声速工况下的弯曲/扭转耦合及控制面颤振特性。Lauten等［6］对X-15的全动平尾缩比模型进行了风洞试验研究，在马赫数Ma=6.86下未观测到任何颤振情况。Goetz［7］采用三阶活塞理论对双楔形翼型进行了颤振分析，结果表明，增加翼型厚度和迎角会导致颤振边界下降。

直到20世纪90年代美国实施国家空天飞机计划(NASP），掀起了吸气式超高速飞行器研究的热潮。Heeg等［8］对6种超高速飞行器全动尾翼模型进行了风洞试验。结果表明，试验颤振马赫数均在15以上，试验颤振速度均低于二阶活塞理论的计算值；薄翼型更易发生颤振，且前缘钝化有利于提高颤振速度；流场中的非线性现象容易引发极限环。

NASP前期，研究者多依赖于非定常活塞理论、激波膨胀波法、牛顿理论等进行初步分析。Liu等［9］提出了统一升力面理论，结合了活塞理论和升力面理论的优势，考虑了厚度和上洗流的影响，对平直机翼的颤振特性进行了研究。Chen和Cao［10］提出了当地流活塞理论，以物面当地速度代替来流速度、以附加下洗代替下洗速度，扩大了活塞理论应用的迎角和厚度范围。Yang和Song［11］将当地流活塞理论用于超音速工况下大迎角翼面颤振分析，并将计算结果与风洞试验结果进行对比。结果表明，当地流活塞理论适用于超声速工况下的大迎角有限厚度翼面颤振问题。

随着计算能力的不断提高，NASP后期，非定常气动力的计算方法逐渐由依靠理论分析转向了与计算流体力学(CFD）相结合的方法。Mcnamara等［12-14］基于活塞理论、牛顿理论、激波膨胀波等非定常气动力方法计算了双楔形翼型的气弹稳定性，并与CFD方法进行了对比。结果表明，一阶、二阶活塞理论的误差较大，且黏性效应对于二维薄翼型的影响可忽略。Zhang和Ye［15］发展了一种基于CFD的当地流活塞理论，并应用于高超声速的气动弹性分析，该方法权衡了计算精度和计算效率，可用于大迎角和高马赫数工况的气动力分析。

此后，随着计算能力的提高和研究的深入，气动学科采用CFD方法、结构学科采用有限元分析(FEA）方法的气动弹性分析方法成为了新的研究趋势，主要研究对象为X-30、X-33、X-34和X-43验证机。Gupta等［16-18］对X-43A进行了气动弹性分析，非定常气动力分别采用了基于活塞理论、CFD以及参数识别等方法，结构采用有限元建模，结果表明，飞行器在工作区间内不会发生颤振现象。Zeng等［19］研究了鸭式布局超高速无人机的气动伺服弹性问题。结果表明，翼轴弹性对耦合系统稳定性有重要作用，增加翼轴刚度有利于抑制颤振。Shi等［20］采用当地流活塞理论，研究了转轴位置对舵面气动弹性的影响。结果表明，随着转轴位置后移，舵面失稳形态由颤振变为静发散。

由于超高速舵面多使用超薄对称楔形翼型，在气动载荷作用下，翼面变形会导致翼型形状由直线变为曲线，常规的理论算法或是基于CFD的数值算法无法考虑翼型的变形。对于结构模型的求解，理论方法假设结构在变形平衡位置附近做小幅振动，采用线性结构动力学方程。近年来，国内外很多学者提出了CFD/CSD耦合方法，开始将其应用在亚/超飞行器气弹分析中［21-23］。CFD采用Euler或N-S流动控制方程，CSD采用非线性结构动力学方程，基于CFD/CSD耦合分析精度高，可用于复杂问题的研究。

到目前为止，关于全动尾翼的气弹研究大多处在较小迎角范围内，尚未见到针对大迎角的气弹特性研究。事实上，超高速飞行过程中，在大气紊流等外部干扰的作用下，飞行器会出现大迎角工作状态，此时尾翼需要大角度偏转，气动力也会急剧增大，带来的气动弹性问题也更为严重。综上，本文针对超高速飞行器气动弹性特点最明显的全动尾翼，采用CFD/CSD/CTD耦合的方法建立其气动/结构/热耦合模型，计算并分析了不同迎角下的气动弹性响应，重点研究大迎角下平尾的气动响应和结构变形特点。

1 物理模型

1.1 气动模型

本文全动尾翼采用双楔形翼型，翼面形状参考NASP中X-30的全动尾翼［8］。翼面形状、翼型及坐标系定义如图1所示：c r为翼根弦长；c为翼型长度；t为翼型厚度；V∞为来流速度。z轴垂直于x轴和y轴，按右手定则给定正方向，详细参数如表1所示。

图1 全动尾翼模型
Fig.1 Model of all-movable tail

表1 全动尾翼几何参数
Table 1 Geometric parameters of all-movable tail

1.2 结构模型

本文超高速飞行器方案采用转轴式全动平尾，即平尾与转轴组成整体相对于机身内的轴承中心轴线旋转。按转轴的形式分为直轴式和斜轴式。直轴式的转轴垂直于飞行器对称轴线，如图2所示。

图2 转轴形式
Fig.2 Form of rotating shaft

直轴式容易布置在机身内，操纵机构重量较轻，X-43A的全动平尾选用的是直轴式［16］。超声速翼型焦点靠后，为减小铰链力矩，通常将轴线布置在平均气动弦长的40%左右处。本文将转轴位置定于翼根65%处，换算位置在平均气动弦长的43%处。尾翼/翼轴系统结构如图3所示。

图3 尾翼/翼轴系统
Fig.3 System of tail/rotating shaft

全动平尾结构形式采用比刚度较高的蒙皮全高度蜂窝结构。为了抵抗超高速飞行带来的结构温升，蜂窝多采用耐高温合金。转轴式全动平尾将分散的受力形式转成由翼轴集中承受，因此对于转轴的结构刚度要求较高。本文翼轴为钛合金材料，翼面蒙皮为镍合金材料，尾翼内部为全高度镍合金蜂窝结构。钛合金和镍合金的材料参数如表2和表3所示。内部蜂窝结构的等效弹性常数如表4所示。

表2 钛合金材料参数
Table 2 Parameters of titanium alloy material

表3 镍合金材料参数
Table 3 Parameters of nickel alloy material

表4 等效弹性常数
Table 4 Equivalent elastic constants

2 计算方法

CFD/CSD/CTD耦合进行气动弹性研究，其计算过程可分为3部分：流体场气动力/气动热计算、固体场热响应/热应力/热应变计算、流固边界数据传递和网格变形计算。流体场气动力/气动热计算通过CFD模块实现，固体场热响应/热应力/热应变计算通过CSD/CTD耦合模块实现。

CFD求解采用基于标准k-ε湍流模型的雷诺平均Navier-Stokes方程的有限体积法，在近壁区域采用壁面函数进行修正；对于空间离散方法，扩散项无黏通量项采用AUSM通量差分分裂格式，黏性通量项采用中心差分格式，对流项使用二阶迎风格式；时间积分使用高斯-赛德尔方法。CSD方面采用牛顿-拉普森迭代法求解非线性结构动力学方程。CTD方面采用加权残差伽辽金法求解基于傅里叶定律的热传导方程。通过热-结构刚度矩阵实现热-结构控制方程的耦合。

2.1 耦合界面相容条件

CFD/CSD/CTD耦合求解时，耦合界面上网格通常不匹配，为保证流体-结构-热耦合的一致性，耦合界面的数据传递需要满足力学平衡条件、几何相容条件、热流平衡条件和温度连续条件，即

式(1）为力平衡条件，σS和σF为结构应力张量和流体黏性应力张量，n为耦合界面单元的单位外法向量；p为流体压力。式(2）为几何相容条件，u S和u F为耦合界面上的结构位移矢量和流体位移矢量。式(3）为热流平衡条件，Q S和Q F为耦合界面上的结构热流量和流体热流量。式(4）为温度连续条件，T S和T F为耦合界面上的结构温度和流体温度。

2.2 数据交换方法

采用守恒插值方法［24］实现耦合界面的数据交换。该方法在数据交互过程中满足能量守恒，包括力学能量守恒和热力学能量守恒。

力学能量守恒即耦合界面上流体力、固体力在界面位移上虚功相等，用式(5）表示：

式中：δu和δu为耦合界面上固体f s和f f流体虚位移；f s和f f为耦合界面上固体和流体表面力。

热力学能量守恒即耦合界面上流体边界和结构边界热流量相等，用式(6）表示：

式中：S为整个耦合界面；q s和q f为耦合界面上固体和流体热流密度。

耦合界面的节点-单元搜索算法采用桶式搜索［25-26］和强力搜索相结合的方法。

2.3 动网格技术

本文网格运动采用基于壁面距离的扩散光顺方法，扩散方程为）

式中：为哈密顿算子；γ为扩散系数；u为网格位移速度。γ与正则壁面距离Δd及扩散参数D有关：

D越大，远离壁面处可吸收更多网格变形，壁面附近网格变形越小，从而保证运动边界附近的网格质量。

2.4 CFD/CSD/CTD耦合计算流程

本文采用CFD/CSD/CTD耦合的方法进行气动-结构-热耦合计算，应用插值接口模块进行耦合界面信息传递。在每个物理时间步内进行预估-校正迭代计算(见图4），使得流体场与固体场在时间积分上达到同步，提高了计算精度，取得了较好的结果。

图4 预估-校正迭代方法
Fig.4 Prediction-correction iterative method

3 算例验证

选取圆管气动加热与结构传热试验［27］作为算例，验证本文计算方法的准确性，该算例已被多次用于验证高超声速流体-固-热多场耦合方法。试验圆管外半径38.1 mm，厚度12.7 mm，材料为AISI321不锈钢，ρ=8 030 kg/m3，E=206 GPa，λ=0.3，α=16.8(10-6/K，k=16.27 W/(m·K），C=502.48 J/(kg·K）。试验来流马赫数6.47，温度241.5 K，压强648.1 Pa；初始壁面温度294.4 K。

流体与固体网格交界面为耦合界面，通过耦合界面进行数据传递实现多场耦合。图5为计算网格，图6和图7分别为圆管表面压力分布归一化(P/P0）和热流分布归一化(q/q0）的计算结果与试验结果的对比，P和q分别为当地压强和当地热流；P0和q0分别为驻点压强和驻点热流；θ为界面到圆心连线与x轴夹角。可以看出，二者的一致性很好，验证了本文计算方法的准确性。

图5 计算网格
Fig.5 Computational grid

图6 表面压力分布对比
Fig.6 Comparison of surface pressure distribution

图7 表面热流分布对比
Fig.7 Comparison of surface heat flux distribution

4 结果与分析

4.1 计算条件

计算来流马赫数Ma=5，高度为30 km，大气温度为226.5 K。尾翼上下表面、前后缘及翼端面设为耦合面，计算时间步长设为0.001 s。翼轴末端采用固定边界条件。分别计算尾翼迎角为10°～40°，间隔10°的情况，依次分析尾翼的气动响应与结构变形特性。

4.2 气动力响应分析

图8为不同迎角下升力系数CL和阻力系数CD变化趋势。可以看出，10°～40°迎角时尾翼的升阻力系数均出现波动，且无相位差，同时达到波峰和波谷，整体呈收敛趋势，收敛值相对初始值均减小。

表5为不同迎角下气动力响应的统计值。从图8和表5可看出，迎角越大，初始振幅越大，气动力减小的比例越大，随时间衰减得越快。

图8 不同迎角下气动力响应曲线
Fig.8 Aerodynamic response curves at different angles of attack

表5 不同迎角下气动力响应
Table 5 Aerodynamic response at different angles of attack

4.3 结构变形后流场分析

图9为不同迎角下初始流场0 s与平衡位置(2.55 s，10°）、(1.91 s，20°）、(1.5 s，30°）、(1.49 s，40°）所对应的尾翼下表面压力系数Cp分布云图。

从图9可以看出结构变形导致下表面的压力分布发生变化。相比初始时刻，10°时平衡位置下表面后缘低压膨胀区压力较低，且前缘高压区向翼根方向移动，使得翼根前缘压力增大而翼尖前缘压力减小。与10°迎角相似，20°迎角时平衡位置尾翼下表面后缘膨胀区低压区域扩大，尾翼前缘高压区更加靠近前缘，使得翼尖处的高压区域明显缩小。与10°、20°迎角不同，30°迎角时初始位置的高压区不再集中在尾翼前缘，而是分布在中部和翼根前缘两片区域。平衡位置时后缘膨胀区压力减小，中部和翼尖附近高压区向翼根前缘移动，因此下表面仅翼根前缘小部分区域压力增大，其他区域压力均减小。与30°迎角相似，40°迎角时平衡位置后缘膨胀区压力减小，且靠近中部的区域下降尤为明显。翼根中部高压区向前缘移动，下表面仅翼根前缘压力增大，其他区域压力均减小。

图9 表面压力系数分布云图
Fig.9 Distribution contour of surface pressure coefficients

图10为10°迎角初始位置与平衡位置尾翼不同展向站位处的下表面压力系数分布曲线对比，分别为沿y轴正向20%w、50%w和80%w处，w为尾翼展长，l为当地弦长。可以看出，在20%w处，平衡位置相比初始位置，仅前缘处压力小幅增大，归一化后的当地弦长x/l＞0.5时压力均减小，而在50%w和80%w处，平衡位置的压力均低于初始位置，且靠近翼尖处降低更多。可见翼根前缘处压力增大的效果并不明显，无法弥补翼尖处的压力损失，从而造成整体升力系数降低。

图10 下表面压力系数分布曲线
Fig.10 Distribution curves lower surface pressure coefficients

迎角20°及以上，尾翼变形较大，各截面位置改变较大，以流场压力分布的正视图和后视图观察更为直接。图11所示为尾翼变形后流场压力分布的正视图(front）和后视图(back）。

从图11可以看出，相比初始位置，平衡时尾翼发生弯曲/扭转耦合现象。从正视图可以看出，展向位置小于50%w的部分前缘向上偏转，相对迎角增大，从而使得激波偏折角增大，波后压力增大，因此造成高压区前移。展向位置大于50%w的部分前缘向下偏转，相对迎角减小，相应的激波后压力减小。从后视图可以看出，尾翼后缘均向上偏转，相对迎角增加，从而增大偏折角，但此时前方激波经过最大厚度处产生膨胀波，偏折角增大反而造成压力减小。

图11 流场压力分布正视图和后视图
Fig.11 Front and back view of pressure distribution in flow field

此外，从正视图和后视图均可看出翼尖部分向上弯曲，从而减小了翼尖部分的相对迎角，耦合翼尖部分的向下偏转，使得此处压力大幅降低。迎角越大，变化幅度越大。综上，结构变形导致整体压力减小，从而减小升力系数。

图12给出了40°迎角时一个典型周期内尾翼结构变形以及尾翼表面压力分布云图的变化过程。取t为0.03～0.08 s，间隔0.01 s的6个时刻。可以看出尾翼前缘扭转和翼尖弯曲变形。变形最大时对应的升力系数最小，变形最小时升力系数最大。

图12 尾翼结构变形和压力分布云图
Fig.12 Structure deformation and pressure distribution contour of tail

4.4 结构变形分析

图13为不同迎角下变形最大时对应的等效应力(Von Mises Stress）分布云图。可以看出，靠近翼轴处应力较集中，翼轴上应力较小。迎角越大，最大应力越大。尾翼最大应力在30°迎角时达到1.2 GPa，已达到所用镍合金材料的屈服强度极限，理论上材料已破坏。应在结构方案设计时在翼轴与内部蜂窝接触部位附近加强，或在控制方案设计时限制全动平尾的工作角度。

图14为10°迎角时尾翼3个点的位移d随时间的变化曲线。可以看出，B点变形最大，A点次之，C点最小，且3条曲线间不存在相位差。为了更清楚地表示尾翼变形，定义位移差如下：

式中：Lz为沿z轴位移；r A和r C为尾翼的扭转和弯曲。

图15为不同迎角下r A和r C随时间的变化曲线。可以看出，尾翼存在弯曲/扭转耦合现象，弯曲/扭转变形随时间逐渐减小，最终达到平衡。迎角越大，初始变形越大，收敛得越快。

图13 结构应力分布云图
Fig.13 Structural stress distribution contour

图14 10°迎角位移随时间变化曲线
Fig.14 Displacement curves at 10°angle of attack over time

图15 不同迎角下r A和r C随时间变化曲线
Fig.15 Curves of r A and r C at different angles of attack over time

图16为40°迎角时A点的温度与热流随时间的变化曲线。可以看出，随着时间的推进，气动加热导致结构温度逐渐升高，耦合界面两侧的温差减小，热流随之降低，温升逐渐趋缓。由于存在热辐射，耦合界面两侧的温差不可能为0，故热流趋于一个稳定值，气动-结构-热耦合达到平衡状态。

图16 A点温度和热流变化曲线
Fig.16 Temperature and heat flux curves of point A

图17为40°迎角平衡位置时上下表面的静温(Static temperature）分布云图。可以看出，此时下表面(迎风面）最高温度达1 290 K，边缘处温度较低，上表面(背风面）最低温度达157 K，边缘处温度较高。为研究气动力、气动热对结构变形的影响，分别对结构施加3种载荷：耦合载荷(气动力气动热同时加载）、气动热载荷和气动力载荷。选取B点沿y轴位移LyB表征结构轴向变形，3种载荷对结构轴向变形的影响如图18所示。可以看出，结构轴向变形主要由气动热引起，气动力影响很小，可忽略不计。选取B点沿z轴位移LzB表征结构法向变形，3种载荷对结构法向变形的影响如图19所示。可以看出，结构法向变形由气动力和气动热共同引起。

图17 40°迎角温度分布云图
Fig.17 Temperature distribution contour at 40°angle of attack

图18 结构轴向变形
Fig.18 Structure deformation along shaft direction

图19 结构法向变形
Fig.19 Structure deformation along normal direction

5 结 论

针对超高速飞行器全动平尾大迎角飞行时需大角度偏转全动尾翼配平而带来的尾翼气动弹性问题，采用计算流体力学/计算固体力学/计算热力学(CFD/CSD/CTD）耦合方法分析了一种超高速飞行器全动平尾的气动弹性响应特性，重点针对大迎角下的气动响应及结构变形特点进行研究。结果表明：

1）CFD/CSD/CTD耦合方法可以有效地分析超高速飞行器全动平尾的气动弹性特性，能够用于研究全动平尾大迎角下的气动响应及结构变形特点。

2）各迎角时的气动力系数曲线均出现波动，随时间变化逐渐衰减至平衡位置。迎角越大，气动力系数曲线初始振幅越大，气动力系数减小的比例越大，随时间衰减得越快。

3）全动平尾存在弯曲/扭转耦合现象，导致下表面压力分布发生变化，靠近翼根部分前缘上偏，压力增大；靠近翼尖部分前缘下偏，压力减小；后缘均上偏，压力减小。结构变形使得平尾表面整体气动压力减小，升力系数降低，迎角越大现象越明显。

4）靠近翼轴处应力较为集中，而翼轴上的应力较小。全动平尾最大应力在迎角30°时达到1.2 GPa，已达到所用镍合金材料的屈服强度极限。应在结构方案设计时在翼轴与全动平尾接触部位附近加强，或在控制方案设计时限制全动平尾的工作角度。

5）全动平尾结构发生轴向与法向变形，其中轴向变形主要由气动热引起，法向变形由气动力和气动热共同引起。

6）热力载荷耦合作用对全动平尾部件气动结构特性的影响在超高速飞行器总体设计时需重点关注，同时飞行器总体设计需对弹性气动修正和结构热防护设计等方面提出设计约束。