基于非定常曲面涡格法的气动灵敏度和线性化分析

陈金池, 黄研昕

(中国航空工业空气动力研究院高速高雷诺数实验室, 辽宁沈阳 110034)

引 言

颤振分析属于气动弹性动稳定性问题, 是高空长航时(high altitude long endurance, HALE)无人机设计的关键。这类无人机普遍采用大展弦比机翼设计, 具有结构质量小、 柔性大的特点, 在气动载荷作用下容易发生大幅弹性变形。大变形会引起结构刚度的非线性变化, 还会使气动面发生曲面变形, 影响气动载荷的大小和分布, 进而导致颤振特性与小变形时显著不同[1-3]。

传统的概念设计往往在小变形范围内进行线性颤振分析, 采用偶极子格网法(double-lattice method, DLM)等平面气动模型即可获得良好的精度[4-5]。然而, 平面气动模型不能准确描述曲面变形的影响, 在大柔性机翼的分析中存在气动载荷失真, 颤振预测不准的问题[6]。

与平面气动模型相比, 非定常曲面涡格法(unsteady vortex-lattice method, UVLM)通过四边形涡格模拟气动面及尾流区域, 并通过对涡格的更新实现对曲面气动面、 曲面尾流区域的模拟。在此基础上, UVLM还可以结合自由尾涡模型[7], 考虑包括翻转在内的复杂尾流外形(见图1)。由于建模简单、 计算效率高、 对流场外形捕捉较准确, UVLM在大变形颤振分析中具有突出优势[8-9], 适用于概念设计阶段。近年来, 许多学者还对UVLM在阵风响应[10-11]、 颤振抑制[12]、 可变形机翼飞行器[13-14]等领域的应用进行了研究。此外, 片条理论和CFD技术也具备曲面气动建模的能力[15-16], 但分别受限于二维假设和计算成本, 在概念设计中存在限制。

图1 UVLM自由尾涡Fig. 1 Free wake of UVLM

标准的UVLM气动载荷求解是时域推进形式, 在处理大型颤振计算时会涉及大量的迭代及气动载荷更新, 这无疑会削弱UVLM在计算效率方面的优势。为了加快求解速度, 可以在大变形平衡位置对UVLM的状态空间方程进行线性化处理, 利用气动载荷灵敏度实现气动载荷的快速更新[17]。同时, 线性化的气动载荷方程也可以与线性结构动力学方程耦合, 利用特征值进行颤振分析[18-20], 提高计算效率。尽管在平衡位置附近进行线性分析, 但UVLM模拟的尾流区域仍可表现为曲面状态, 气动载荷计算相比平面气动模型更加准确。此外, UVLM气动载荷灵敏度和线性化方程还可与其他模块耦合, 为气动弹性控制及优化设计提供参考[21-23]。

作为线性化处理的一部分, 早期的气动载荷灵敏度常常通过有限差分法求解。这种方法对自变量扰动的要求较高, 计算效率较低。近年来, 基于伴随矩阵的灵敏度求解思路备受CFD研究人员的青睐[24-26]。但是, 这种思路涉及方程的变分, 对存储空间要求较高。考虑到UVLM方程的气动载荷与输入量、 状态量存在明确的因式分解关系。通过链式法则即可推导UVLM气动载荷的解析灵敏度。国内外学者针对面元法的气动载荷灵敏度开展过一些研究[27-30]。其中, 针对UVLM的研究主要考虑了气动载荷有关节点形状和位置的灵敏度[21,31], 缺少对节点速度和环量强度的考虑, 对UVLM线性化方程中的Jacobi矩阵设计不够完整。

针对这一情况, 本文将首先构造线性化的UVLM状态空间方程, 明确Jacobi矩阵包含的元素, 并着重推导气动载荷有关节点位置、 节点速度和环量强度的解析灵敏度。在此基础上, 本文通过数值仿真探讨了这些灵敏度矩阵的分布规律, 利用有限差分法验证了解析灵敏度的计算效率和精度。

1 计算方法

1.1 UVLM气动载荷模型

如图2所示, UVLM将机翼的平均气动面离散为一系列结构网格(黑色实线), 在结构网格下游1/4弦长处布置附着涡涡格(红色虚线), 通过对附着涡的更新模拟气动面曲面效应。随着计算过程的推进, 附着涡会从后缘脱落形成尾涡(黑色虚线), 集合附着涡与尾涡的位置即可较准确地反映气动面及流场的几何特征。UVLM在每个涡格内设置涡环, 通过环量强度完成气动载荷计算。

图2 UVLM涡格划分示意图Fig. 2 Discrete vortex panels of UVLM

附着涡涡格的几何中心可视为气动载荷的控制点。在计算过程中, 需要在附着涡涡格的控制点处施加不可穿透边界条件, 用以获得附着涡的环量强度。与平面气动模型不同, UVLM附着涡沿曲面分布, 边界条件中的诱导速度并非竖直向下, 而是沿着涡格的法向量向下。以t+1时刻为例, 不可穿透边界条件可表示为

(1)

其中,Γbv和Γwk为两组列向量, 分别表示附着涡和尾涡的环量强度; 矩阵Wbb和Wbw分别表示附着涡和尾涡关于气动载荷作用点的影响系数, 用以计算诱导速度,Wbb和Wbw中的元素由涡格节点位置决定, 可通过Biot-Savart定律获得;Unc表示附着涡控制点处的非环量速度, 包括刚体运动速度、 弹性振动速度、 空气扰动速度等。上标n表示涡格的法向量方向。

尾流区域的演化包括两部分, 即附着涡的脱落和已有尾涡的扩散。演化方程如下

(2)

其中,Abv和Awk是尾涡涡强演化矩阵;Bbv和Bwk是尾涡节点演化矩阵, 均由1和0组成;ζbv和ζwk分别表示附着涡和尾涡节点的位置,vwk表示尾涡节点的局部速度, 由流场的诱导速度及外部扰动组成。为避免尾涡数量过多带来的额外计算成本, 可在尾流区域达到一定长度时设置截断, 忽略远处尾涡的影响。

在每个控制点处, UVLM的气动载荷升力分量L垂直于非环量速度向上, 诱导阻力D沿非环量速度方向。L和D的矩阵表达式如下

(3)

其中,E1(2)是由0, 1, -1组成的稀疏矩阵, 用以计算有效涡强。G1(2)是对角矩阵, 用以计算投影面积, 由附着涡的几何信息和局部攻角决定。U1(2)(3)是诱导速度对角矩阵。根据本文需要,G1(2)及U1(2)(3)的表达形式如下

(4)

1.2 线性化状态空间方程

将方程(1)～(4)整合, 可得到UVLM气动载荷在离散时刻的状态空间形式

fA(xt+1,ut+1)=gA(xt,ut)
yt+1=hA(xt+1,ut+1)

(5)

其中,u,y,x分别表示方程的输入量、 输出量、 状态量, 分别为

(6)

方程(5)可以与结构动力学方程耦合, 形成气动弹性方程。

图3 UVLM气动力线性化过程的参数关系图Fig. 3 Parameter relationship of aerodynamic linearization by the UVLM

将方程(5)中的变量表示为在平衡位置x0处的基准量与扰动量之和

(7)

再将方程(7)代入方程(1)～(3), 基于上述参数关系可得到线性化的状态空间方程

(8)

其中, Δt为时间步长。

将方程(8)中的Jacobi矩阵整合, 线性化的状态空间方程可表示为

EΔxt+1=AΔxt+BΔut+1, Δyt+1=CΔxt+1+DΔut+1

(9)

1.3 气动载荷解析灵敏度

(10)

(11)

其中

(12)

(13)

(14)

其中, 下标i表示附着涡涡格节点索引, 每个元素表示节点i发生扰动对附着涡涡格k的影响。

按类似的方法可推导方程(12)～(13)中其余偏导项, 并获得方程(11)中气动载荷解析灵敏度的表达式, 在此不一一赘述。相比其他计算方法, 采用链式法则推导的解析灵敏度可直接使用矩阵运算符, 避免了数值计算中的循环过程。

2 计算结果

将长直机翼置于均匀分布的流场中, 忽略弹性变形和重力的影响, 通过气动载荷系数验证UVLM建模的准确性, 再对平衡状态下的机翼施加小扰动, 验证解析灵敏度的计算精度及效率。

2.1 气动载荷验证

构建5个攻角为5°, 弦长为1 m的长直机翼, 这些机翼的展弦比分别为4, 8, 12, 20, 100。将这些机翼依次置于速度50 m/s的均匀流场中, 通过UVLM模型获得升力系数与诱导阻力系数, 并与Katz等[32]提供的参考值进行对比。计算中, 对UVLM尾流区域在9 m处进行截断, 忽略远处尾涡的影响。

在第1步迭代计算中, 流场从静止开始起风, 机翼会受到一个较大的脉动载荷, 导致升力系数出现较大的初始值。随着尾涡的生成, 这一脉动载荷带来的影响逐渐减小。为便于比较, 本文将初始升力系数设置为0.5。对比结果如图4所示, UVLM升力在迭代过程中逐渐达到与参考值接近的稳定状态, 并且迭代步数随着展弦比的增大而增加, 稳定时的升力系数也随着展弦比的增大而接近5°攻角下的二维翼型计算值0.55。

图4 UVLM升力系数与展弦比的关系Fig. 4 Variation of lift coefficient with aspect ratio

在诱导阻力系数的对比中, 选取展弦比为8的机翼, 分析了诱导阻力系数及分量的变化情况。图5显示计算结果与Katz和Poltkin的结果吻合, 进一步验证了UVLM建模的准确性。

图5 UVLM诱导阻力系数及分量瞬态变化(展弦比8)Fig. 5 Variation of induced drag coefficient and components at an aspect ratio of 8

2.2 气动载荷灵敏度验证

将机翼的平均气动面划分为20个展向单元和2个弦向单元。在UVLM时间推进计算中, 选取0.025 s作为UVLM时间推进步长。计算结果表明, 算例机翼的气动载荷在起风1.5 s后达到稳定。因此, 选取2 s时刻稳定状态作为灵敏度计算的平衡状态。在机翼向下游方向30 m处设置截断, 忽略30 m外的尾涡影响。尾流区域如图6所示。计算中, 附着涡及尾涡的编号顺序如图7所示。

图6 UVLM尾涡在2 s时刻的涡格情况Fig. 6 UVLM wake at 2 s

图7 UVLM涡格及节点编号Fig. 7 Numbering scheme of UVLM vortex panels and nodes

(a) Analytical sensitivity

表2 升力关于涡格节点位置灵敏度结果及对比(y)

表3 升力关于涡格节点位置灵敏度结果及对比(z)

表4 升力关于涡格节点速度灵敏度结果及对比(y)

表5 升力关于涡格节点速度灵敏度结果及对比(z)

解析灵敏度与有限差分法的计算时间对比见表6, 采用解析灵敏度的计算效率显著高于有限差分法。本文最大仅涉及40×63的灵敏度矩阵计算, 在这一情况下, 解析解灵敏度的计算效率大约是有限差分法的5～6倍。这种优势在进行更大型的计算时会更加明显, 这是两种方法的理论基础决定的。解析灵敏度在计算时直接使用矩阵运算符, 不需进行任何循环或迭代。而有限差分法作为一种数值计算方法, 必须对一个变量进行扰动, 再通过多次循环计算才能获得灵敏度结果, 过程中的计算消耗远大于解析解。

表6 解析解与有限差分法计算时间对比

3 结论

本文详细描述了UVLM的状态空间建模过程, 并推导了UVLM升力关于附着涡强度以及涡格节点运动信息的解析灵敏度。经验证, 解析灵敏度的数值大小、 分布情况均与有限差分法一致, 各灵敏度项的计算误差不超过0.4%。相比有限差分法, 解析灵敏度直接进行矩阵运算, 避免了扰动量选取不当带来的收敛性问题, 也避免了大量数值计算带来时间消耗, 体现出高精度、 高效、 高鲁棒性的优势。

气动力解析灵敏度可进一步与结构响应方程耦合, 构建气动弹性分析系统。除了可大幅提高大型系统的稳定性分析效率, 解析灵敏度在气动弹性优化设计及流动控制领域都具有良好的应用价值。