2. 哈尔滨理工大学 理学院, 黑龙江 哈尔滨 150080
高空长航时飞行器因可执行大气监测、通信中继、高空侦察、反导防御等多种军用或民用任务, 受到广泛关注。为了提高气动效率, 在追求结构轻量化的同时, 这类飞行器还普遍采用大展弦比的机翼设计[1]。展弦比的增大会导致机翼结构刚度的降低, 在气动载荷的作用下会出现较大变形和振动, 变形和振动反过来又影响机翼表面流动, 产生复杂的气动弹性现象, 如静发散和颤振[2-4]。
对于大展弦比机翼颤振问题, 现有的时域分析方法普遍采用分区思想, 将机翼的非定常气动力计算和结构动力学分析在时间层面上交替进行, 同时在机翼表面上进行数据传递。在每一个时间步内, 首先求解作用在机翼上的气动力载荷, 然后求解机翼的变形和振动响应并更新机翼的形状, 为下一时刻的气动力计算做准备。
大展弦比机翼的非定常气动力计算方法可以分为片条理论、面元法和CFD数值模拟三大类。片条理论的主要思想是将大展弦比机翼沿展向划分为若干片条单元, 采用Theodoren模型[5]、ONERA模型[6]以及有限状态理论[7]等二维翼型气动理论求解每个片条的气动力, 根据展向片条之间的气动干扰进行三维修正, 由此获得整块机翼受到的气动力载荷[8]。涡格法[9]是一种具有代表性的面元法, 该方法由势流方程导出, 将机翼离散成沿弦向和展向分布的附着涡, 并在机翼后缘布置自由涡, 利用翼面无穿透条件求解涡强, 积分得到机翼的气动特性[10]。随着计算机技术的发展, 机翼颤振时域分析逐渐采用以Euler方程或N-S方程为基础的CFD数值模拟方法来计算非定常气动力[11-12], 这类方法直接求解流动基本方程, 除了气动力以外还能获得瞬时流场信息。由此可得, 目前常用的大展弦比机翼的气动力求解方法大多基于势流方程或Euler方程, 忽略了气体黏性对气动力和流场的影响, 应用到颤振时域分析时可能无法得到合理的计算结果。
在大展弦比机翼颤振时域分析过程中, 研究人员常用梁模型来求解机翼的结构动力学响应, 该方法根据机翼的结构外形特征, 将机翼抽象为非线性细长梁, 再利用有限元方法进行数值求解, 如本征梁模型[13]和梁单元模型[14]等。Hodges-Dowell梁[13]是一种典型的本征梁模型, 该模型直接在变形后的梁坐标系下建立动力学变量特征, 极大地减少了坐标的转换。但是, 除了梁的速度以外, 该模型还需要预先求解梁的内力, 未知量个数大幅增加。共轭旋转(CR)法[14]是最常用的梁单元模型之一, 该方法的核心思想是将梁的变形分解为刚体的平动和转动, 以及单元坐标系内的线弹性变形, 动力学平衡方程在单元坐标系内建立, 通过坐标转换矩阵构造总体切线刚度矩阵[15]。CR法将梁变形本身的位移作为变量, 物理意义明确, 缺点是计算过程中需要不断更新切线刚度矩阵, 计算效率有所降低。除了梁模型以外, 另一类常用的结构动力学分析方法是模态法[16], 该方法以结构的主要固有振型为广义坐标系, 基于模态坐标建立动力学方程。相比于梁模型, 模态法的优势在于建立的方程形式更为简单, 并且广义质量和广义刚度矩阵都是对角矩阵, 计算效率也更高。
综上所述, 大展弦比机翼的颤振时域分析主要关注颤振速度和频率[17-18], 部分研究也有涉及根部攻角[19]、展弦比[20]、扭转刚度[12]等参数对机翼颤振特性的影响。机翼的变形和振动响应有助于认识和揭示颤振诱发机理, 目前仍然缺乏详细的研究。本文以展弦比为16的NACA0012弹性平直机翼为研究对象, 结合三维雷诺平均N-S方程、S-A湍流模型[21]和机翼动力学方程, 建立了机翼颤振时域分析方法。利用该方法对机翼的颤振特性、颤振过程中的变形和振动响应, 以及瞬时流场进行数值模拟和分析, 探究机翼发生颤振的原因, 研究材料刚度和密度对机翼颤振速度和频率的影响。
1 大展弦比机翼颤振时域分析方法 1.1 非定场气动力计算方法本文采用自研CFD软件求解展弦比为16的NACA0012弹性平直机翼的非定常气动力, 其控制方程为积分形式的守恒型三维非定常雷诺平均N-S方程, 直角坐标系下的形式为
|
(1) |
式中, 守恒向量项Q为[ρ ρu ρv ρw e]T, ρ为气体密度, (u v w)为直角坐标系下速度分量, e为单位质量气体总能量。F, G和H为3个方向的对流通量矢量项, Fv, Gv和Hv为3个方向的黏性通量矢量项。Ω表示控制体。对流通量项采用迎风型的Roe格式进行离散, 黏性通量项采用标准的二阶中心格式进行离散, 时间推进采用双时间步法, 为了封闭方程引入S-A湍流模型[21]。
1.2 机翼动力学方程本文以展弦比为16的NACA0012弹性平直机翼为研究对象, 通过模态法[16]求解机翼的变形和振动响应
|
(2) |
式中：Δr为机翼的变形矢量; N为所选取的固有振型阶数; Φi(x, y, z)为第i阶模态的固有振型矢量; ξi(t)为第i阶模态的广义位移。
机翼动力学方程可用矩阵形式表示为
|
(3) |
式中：
为了便于时域求解, 引入状态变量E=[ξ1, ξ2, …, ξN, 
|
(4) |
式中
|
采用基于预估-校正技术的四阶杂交线性多步法[22]对(4)式进行时域推进求解, 将(4)式展开为(5)式, Δt为计算时间步长。该方法通过气动力插值技术将气动力部分由隐式格式转变为显式格式, 实现了每个物理时间步内仅需求解一次非定常流场, 兼具隐式稳定性好和显式效率高的优点。
|
(5) |
通过有限元模态分析得到大展弦比机翼的结构参数矩阵(M, G和K)和固有振型, 通过CFD数值模拟得到气动力(F)。由于CFD翼面网格点和模态分析控制点一般是不重合的, 在每一个时间步内, 通过结构固有振型在翼面处进行数据传递, 即压力分布通过振型转化为广义气动力, 进而求解机翼动力学方程获得广义位移, 广义位移通过振型再转化为翼面物理位移, 作为CFD下一时刻求解的边界条件。模态分析得到的结构固有振型通过径向基函数(RBF)插值法[23]转换为CFD翼面网格点描述下的振型分布, 并采用基于径向基函数插值的网格变形技术[23]实现机翼变形后CFD网格的运动。
1.3 算例验证AGARD445.6机翼颤振特性风洞实验是由NASA兰利研究中心在其跨声速动态风洞中完成的[24]。目前, 该机翼模型已经成为国际上颤振计算程序考核的标准算例[25]。
该机翼模型沿流向翼型为NACA65A004, 展弦比为1.644, 根梢比为0.659 2, 1/4弦线后掠角为45°, 来流攻角为0°。机翼模型采用桃花心木制成, 颤振分析取前四阶振动模态, 其中, 第一阶和第三阶模态为弯曲模态, 第二阶和第四阶模态为扭转模态。
图 1给出了本文仿真得到的AGARD445.6机翼颤振速度和频率随马赫数Ma的变化关系, 并与实验结果[24]和文献结果[25]比较。可以看出, 颤振速度和频率均出现了明显的凹坑现象。本文结果在亚声速状态(Ma=0.499, 0.678, 0.901和0.96)时和实验结果吻合很好, 在超声速区(Ma=1.072和1.141)略高于实验结果, 但是仍然明显优于文献[25]的结果, 表明本文建立的机翼颤振时域分析方法准确可靠。
|
| 图 1 AGARD445.6机翼的颤振速度和频率 |
展弦比为16的NACA0012弹性平直机翼颤振是本文研究的重点。机翼模型的几何参数、材料参数和飞行环境如表 1所示。颤振分析取前五阶振动模态, 表 2给出了模态名称和固有频率, 机翼的第一、二和五阶模态为垂向弯曲模态, 第三阶为扭转模态, 第四阶为弦向弯曲模态。
| 类型 | 名称 | 数值 |
| 几何参数 | 弦长c/m 展弦比 |
1 16 |
| 材料参数 | 垂向弯曲刚度/(N·m2) 弦向弯曲刚度/(N·m2) 扭转刚度/(N·m2) 质量/kg |
2×104 4×106 1×104 12 |
| 飞行环境 | 飞行高度/km 空气密度ρ/(kg·m-3) 根部攻角/° |
20 0.088 9 0 |
| 阶数 | 模态名称 | 固有频率/(rad·s-1) |
| 1 | 垂向弯曲 | 2.25 |
| 2 | 垂向弯曲 | 14.101 |
| 3 | 扭转 | 30.671 |
| 4 | 弦向弯曲 | 31.344 |
| 5 | 垂向弯曲 | 39.478 |
不同来流速度U∞时大展弦比NACA0012弹性平直机翼的广义位移时间响应曲线如图 2所示。可以看出, 第一和二阶振动模态广义位移(ξ1和ξ2)在U∞为20 m/s和28.5 m/s时出现了一定幅度的静漂移, 说明由于展弦比较大, 且材质轻柔, 机翼发生了静弯曲变形。当U∞为20 m/s时, 各阶模态广义位移响应均表现为振动衰减。当U∞提高至28.5 m/s时, 广义位移均表现为近似等幅响应。进一步提高U∞至40 m/s时, 广义位移表现出明显的振动发散趋势。综合以上计算结果, 可以得出, 机翼在U∞为28.5 m/s时处于临界稳定状态, 即颤振速度Uf约为28.5 m/s。
|
| 图 2 不同来流速度U∞时大展弦比NACA0012弹性平直机翼的广义位移时间响应曲线 |
利用结构固有振型将机翼的广义位移转化为翼面物理位移。图 3a)给出了颤振速度为28.5 m/s)时大展弦比NACA0012弹性平直机翼翼尖位置的垂向弯曲位移响应, 图 3b)则给出了对应的FFT频谱。可以看出, 位移响应主频率为22 rad/s。因此, 可以得出, 机翼的颤振频率ωf为22 rad/s。
|
| 图 3 U∞=28.5 m/s时大展弦比NACA0012弹性平直机翼翼尖位置的垂向弯曲位移响应及其FFT频谱 |
图 4给出了来流速度U∞=40 m/s时机翼的整体变形特征以及z=15 m展向位置的压力云图。可以看出, 当t=0.00 s时, 机翼未发生变形, 压力关于机翼弦向平面对称, 从前缘到尾缘, 压力先减小然后增大。随着时间t的推进, 机翼的振动幅值逐渐增大, 变形愈发明显, 其局部剖面的形状和扭转角都发生变化, 机翼背风面附近的负压区也逐渐增大。
|
| 图 4 U∞=40 m/s时大展弦比NACA0012弹性平直机翼的整体变形特征以及z=15 m展向位置的压力云图 |
频率重合理论的颤振机理为: 机翼2个振动模态的频率随着来流速度的增大而相互接近, 这时这两阶模态的耦合作用增强, 最终引发颤振[16, 26]。通过比较和分析大展弦比NACA0012弹性平直机翼的颤振频率(22 rad/s)和各阶振动模态固有频率, 可以看出, 第一阶和第五阶模态频率均与颤振频率相差较大, 颤振频率基本介于第二阶和第三阶模态频率之间, 因此, 很有可能是这两阶模态耦合导致了机翼颤振, 这也符合经典颤振的诱发机理, 即弯扭耦合作用[26]。然而, 第三阶和第四阶模态频率十分接近, 第二~四阶模态对机翼颤振的影响仍不明确, 需要作进一步详细分析。
将第二~四阶振动模态分别固定, 重新仿真得到颤振速度时机翼其他模态广义位移的时间响应曲线, 如图 5所示。图 6给出了第四阶模态固定时机翼翼尖位置的垂向弯曲位移响应及其FFT频谱。可以看出, 当第二阶或第三阶模态固定时, 广义位移均振动衰减, 意味着颤振消失。当第四阶模态固定时, 广义位移和翼尖位移均呈现近似等幅响应, 机翼颤振速度和频率没有变化。以上结果说明, 第二阶和第三阶模态的耦合作用是机翼发生颤振的根本原因。
|
| 图 5 不同振动模态固定时U∞=28. 5 m/s时大展弦比NACA0012弹性平直机翼其他模态广义位移时间响应曲线 |
|
| 图 6 第四阶模态固定时U∞=28.5 m/s时大展弦比NACA0012弹性平直机翼翼尖位置的垂向弯曲位移响应及其FFT频谱 |
材料刚度和密度是机翼设计的重要参数。本节研究材料刚度对大展弦比NACA0012弹性平直机翼颤振速度和频率的影响。将机翼的垂向弯曲刚度、弦向弯曲刚度和扭转刚度等比例同时变化, 选取了6组刚度值, 如表 3所示, 其中RS为刚度比, 表征机翼材料刚度的变化。可以看出, 最小的一组为原刚度的25%, 即RS为0.25, 最大的一组刚度的RS为1.5。相应的, 图 7给出了这6组刚度下机翼前五阶振动模态的固有频率, 随着刚度的增大, 模态固有频率也在逐渐升高。
| 序号 | 垂向弯曲刚度/ (N·m2) | 弦向弯曲刚度/ (N·m2) | 扭转刚度/ (N·m2) | 刚度比RS |
| 1 | 5×103 | 1×106 | 2.5×103 | 0.25 |
| 2 | 1×104 | 2×106 | 5×103 | 0.5 |
| 3 | 1.5×104 | 3×106 | 7.5×103 | 0.75 |
| 4 | 2×104 | 4×106 | 1×104 | 1 |
| 5 | 2.5×104 | 5×106 | 1.25×104 | 1.25 |
| 6 | 3×104 | 6×106 | 1.5×104 | 1.5 |
|
| 图 7 不同刚度比下大展弦比NACA0012弹性平直机翼前五阶振动模态的固有频率 |
图 8给出了不同刚度比下大展弦比NACA0012弹性平直机翼的颤振速度和频率。可以看出, 对于最小的一组刚度, 其颤振速度比原刚度时减小了约50%, 颤振频率减小了约49%。对于最大的一组刚度, 其颤振速度和频率分别提高了约23%和20%。这说明了随着材料刚度的增大, 机翼的颤振速度和频率都会提高。
|
| 图 8 不同刚度比下大展弦比NACA0012弹性平直机翼颤振速度和频率 |
本节研究材料密度对大展弦比NACA0012弹性平直机翼颤振速度和频率的影响。定义ρs/ρ为密度比, 用于表征机翼材料密度ρs的变化。图 9给出了不同密度比下机翼前五阶振动模态的固有频率, 可以看出, 最小的ρs为原材料密度的25%, 最大的为原来的2.5倍。随着ρs/ρ的增大, 模态的固有频率逐渐降低。
|
| 图 9 不同密度比下大展弦比NACA0012弹性平直机翼前五阶振动模态的固有频率 |
图 10给出了不同密度比下大展弦比NACA0012弹性平直机翼的颤振速度和频率。可以看出, 对于最小的材料密度, 原ρs时的颤振速度和频率分别比其降低了约16%和51%。对于最大的ρs, 其颤振速度和频率分别比原ρs下降低了约2%和39%。说明了随着材料密度的增大, 颤振速度会先明显降低然后趋于不变, 颤振频率会一直降低。
|
| 图 10 不同密度比下大展弦比NACA0012弹性平直机翼颤振速度和频率 |
本文采用以三维雷诺平均N-S方程和S-A湍流模型为基础的CFD数值模拟方法来求解展弦比为16的NACA0012弹性平直机翼的非定常气动力, 结合机翼动力学方程, 建立了机翼颤振时域分析方法, 进而对机翼的颤振特性、振动响应、整体变形特征和瞬时流场特征进行了数值模拟研究, 主要得到以下结论:
1) 由于展弦比较大, 且材质轻柔, 机翼会发生一定幅度的静弯曲变形。机翼颤振由第二阶(垂向弯曲)和第三阶(扭转)振动模态的耦合作用所致, 颤振频率也介于这两阶模态的固有频率之间。
2) 通过对具有不同材料刚度和密度的机翼进行数值模拟, 对比其颤振速度和频率, 发现适当地增大材料刚度或者减小材料密度, 可以提高大展弦比机翼的颤振速度和频率。
| [1] |
刘国春. 大展弦比机翼气动外形设计方案研究[J]. 飞机设计, 2011, 31(3): 9-12.
LIU Gouchun. Research on design of high aspect wing shape[J]. Aircraft Design, 2011, 31(3): 9-12. (in Chinese) DOI:10.3969/j.issn.1673-4599.2011.03.003 |
| [2] | ARENA A, LACARBONARA W, MARZOCCA P. Nonlinear aeroelastic formulation and postflutter analysis of flexible high-aspect-ratio wings[J]. Journal of Aircraft, 2013, 50(6): 1748-1764. DOI:10.2514/1.C032145 |
| [3] |
马艳峰, 贺尔铭, 曾宪昂, 等. 基于流固耦合方法的大展弦比机翼非线性颤振特性分析[J]. 西北工业大学学报, 2014, 32(4): 536-541.
MA Yanfen, HE Erming, ZENG Xianang, et al. Studying nonlinear flutter of high-aspect-ratio wing based on fluid solid coupling[J]. Journal of Northwestern Polytechnical University, 2014, 32(4): 536-541. (in Chinese) DOI:10.3969/j.issn.1000-2758.2014.04.012 |
| [4] |
段静波, 周洲, 江涛. 基于曲梁模型的大展弦比大柔性机翼颤振分析[J]. 西北工业大学学报, 2016, 34(5): 774-782.
DUAN Jingbo, ZHOU Zhou, JIANG Tao. A method for flutter of the very flexible wing with high-aspect-ratio based on curve beam model[J]. Journal of Northwestern Polytechnical University, 2016, 34(5): 774-782. (in Chinese) DOI:10.3969/j.issn.1000-2758.2016.05.006 |
| [5] |
谢长川, 吴志刚, 杨超. 大展弦比柔性机翼的气动弹性分析[J]. 北京航空航天大学学报, 2003, 29(12): 1088-1090.
XIE Changchuan, WU Zhigang, YANG Chao. Aeroelastic analysis of flexible large aspect ratio wing[J]. Journal of Beijing University of Aeronautic and Astronautics, 2003, 29(12): 1088-1090. (in Chinese) |
| [6] | TANG D M, DOWELL E H. Comments on the ONERA stall aerodynamic model and its impact on aeroelastic stability[J]. Journal of Fluids and Structures, 1996, 10: 353-366. DOI:10.1006/jfls.1996.0023 |
| [7] | PETERS D A, KARUNAMOORTHY S, CAO W M. Finite state induced flow models; I-two-dimensional thin airfoils[J]. Journal of Aircraft, 1995, 32(2): 313-322. DOI:10.2514/3.46718 |
| [8] |
杨超, 杨澜, 谢长川. 大展弦比柔性机翼气动弹性分析中的气动力方法研究进展[J]. 空气动力学学报, 2018, 36(6): 1009-1019.
YANG Chao, YANG Lan, XIE Changchuan. Development of aerodynamic methods in aeroelastic analysis for high aspect ratio flexible wing[J]. Acta Aerodynamic Sinica, 2018, 36(6): 1009-1019. (in Chinese) |
| [9] | KATZ J, PLOTKIN A. Low-speed aerodynamics: from wing theory to panel methods[M]. New York: McGraw-Hill, 1991. |
| [10] |
孙岩, ANDREAD R, 王运涛, 等. 基于非线性涡格法的快速静气动弹性数值模拟技术[J]. 气体物理, 2020, 5(6): 26-38.
SUN Yan, ANDREA D R, WANG Yuntao, et al. Fast static aeroelasticity simulation approach based on nonlinear vortex lattice method[J]. Physics of Gases, 2020, 5(6): 26-38. (in Chinese) |
| [11] | PATIL M J, HODGES D H. Limit-Cycle oscillations in high-aspect-ratio wings[J]. Journal of Fluid and Structures, 2001, 15: 107-132. DOI:10.1006/jfls.2000.0329 |
| [12] | AFONSO F, LEAL G, VALE J, et al. The effect of stiffness and geometric parameters on the nonlinear aeroelastic performance of high aspect ratio wings[J]. Proceedings of The Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2017, 231(10): 1-27. |
| [13] | HODGES D H. A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams[J]. International Journal of Solids and Structures, 1990, 26(11): 1253-1273. DOI:10.1016/0020-7683(90)90060-9 |
| [14] | CRISFIELD M A. A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements[J]. Computer Methods in Applied Mechanics and Engineering, 1990, 81(2): 131-150. DOI:10.1016/0045-7825(90)90106-V |
| [15] |
王伟, 段卓毅, 耿建中, 等. 基于CR理论的大柔性机翼几何非线性结构建模[J]. 航空学报, 2017, 38(增刊1): 72154.
WANG Wei, DUAN Zhuoyi, GENG Jianzhong, et al. Modelling of geometrically nonlinear structure of large flexible wing based on CR theory[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(suppl 1): 72154. (in Chinese) |
| [16] |
叶正寅, 王刚, 张伟伟. 流固耦合力学基础及其应用[M]. 哈尔滨: 哈尔滨工业大学出版社, 2016.
YE Zhengyin, WANG Gang, ZHANG Weiwei. Fundamentals of fluid-structure coupling and its application[M]. Harbin: Harbin Institute of Technology Press, 2016. (in Chinese) |
| [17] | QIAO Shengjun, JIAO Jin, NI Yingge, et al. Effect of stiffness on flutter of composite wings with high aspect ratio[J]. Mathematical Problems in Engineering, 2021, 4: 1-14. |
| [18] | KEIVAN E, MORTEZA D, MOHAMMAD H P, et al. Effects of aeroelastic nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings[J]. Applied Mechanics and Materials, 2012, 110/111/112/113/114/115/116: 4297-4306. |
| [19] | ANDREA A, WALTER L. Nonlinear aeroelastic formulation and postflutter analysis of flexible high-aspect-ratio wings[J]. Journal of Aircraft, 2013, 50(6): 1748-1764. |
| [20] |
刘湘宁, 向锦武. 大展弦比复合材料机翼的非线性颤振分析[J]. 航空学报, 2006, 27(2): 213-218.
LIU Xiangning, XIANG Jinwu. Study of nonlinear flutter of high-aspect-ratio composite wing[J]. Acta Aeronautica et Astronautica Sinica, 2006, 27(2): 213-218. (in Chinese) |
| [21] | SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows[J]. Recherche Aerospatiale, 1994, 1: 5-21. |
| [22] | ZHANG W W, JIANG Y W, YE Z Y. Two better loosely coupled solution algorithms of CFD based aeroelastic simulation[J]. Engineering Applications of Computational Fluid Mechanics, 2007, 1(4): 253-262. |
| [23] | WANG G, MIAN H H, YE Z Y, et al. Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions[J]. AIAA Journal, 2014, 53(4): 1016-1025. |
| [24] | YATES E C J. AGARD Standard aeroelastic configurations for dynamic response I-wing 445.6[R]. NASA TM-1988-0017809, 1988 |
| [25] | SILVA W A, LⅡ B P, CHWALOWSKI P. Evaluation of linear, inviscid, viscous, and reduced-order modelling aeroelastic solutions of the AGARD 445.6 wing using root locus analysis[J]. International Journal of Computational Fluid Dynamics, 2014, 28(3): 122-139. |
| [26] |
郑欣, 刘宇斌, 陈璞, 等. 基于弯扭耦合理论的颤振频率计算方法[J]. 工程力学, 2018, 35(增刊1): 1-5.
ZHENG Xin, LIU Yubin, CHEN Pu, et al. Two-dimensional wing flutter frequency calculation method based on bending torsion coupling theory[J]. Engineering Mechanics, 2018, 35(suppl 1): 1-5. (in Chinese) |
2. School of Science, Harbin University of Science and Technology, Harbin 150080, China
