翼面静气动弹性发散数值分析与试验技术

刘 凯, 苑凯华, 田海涛

(北京机电工程研究所, 北京 100074)

引 言

气动弹性问题一直是飞行器设计中备受关注的重要问题, 在未来的飞行器设计中仍将占有重要地位。现代飞行器设计日益追求高速度、 高机动性, 使得飞行器越来越呈现出轻结构特点, 气动弹性问题也日渐突出[1]。随着高速飞行器研究工作的逐渐深入, 飞行器的性能不断提升, 高速飞行器的结构质量占比需要大幅降低, 使得结构设计人员将面临巨大的压力和挑战[2-8]。显然, 结构设计的余量越来越小, 结构越来越“纤细”, 导致结构的刚性降低, 飞行器的气动弹性影响日益突出。

发散是气动弹性力学中的一个重要问题, 属于静气动弹性稳定性问题。它的现象是: 在一个完全确定的临界风速下, 弹性升力系统受定常升力的作用, 使其扭转变形不断扩大直至破坏[9]。在进行模型临界发散风洞试验时, 可以通过观察模型在受到干扰后是否能够保持稳定来判断发散现象, 即在不断增加风速的过程中, 分别在各个固定风速下给模型一个扰动, 如果在某个风速下受扰动后的模型变形不再回复到初始平衡位置, 并出现变形扩大, 那么这个动压即为发散临界动压[10]。

气动弹性发散边界预测方面, 主要集中为前掠翼风洞试验验证[11-12], 细长体火箭弹采用工程气动力结构近似方法获取弹体的发散边界[13], 工程应用方面往往取足够的安全裕度来保证飞行器不发生静气弹发散[14]。采用CFD/CSD耦合的方法获取结构静发散边界研究很少, 本文提出的气动力/结构耦合的方法可准确预测结构发散边界, 能有效地应用于飞行器气动外形和设计中。

1 数值分析流程

静气动弹性计算中, 在气动力作用下, 结构的平衡方程可写作[15]

Ku=F

其中,K为结构刚度矩阵,u为结构自由度位移矢量,F为气动载荷矢量。

依据初始构型获得的气动力, 根据公式可以计算出相应的结构变形, 结构形状改变影响原有的气动外形, 从而又使气动力发生改变, 即气动力F随着结构外形u改变而变化。由结构变形引起的气动力变化须通过重新进行流场计算才能获得, 所以要求解结构的真实变形, 须采用循环迭代方法求解,如图1所示。

图1 静气动弹性计算流程图Fig. 1 Static aeroelastic calculation flowchart

对于静气弹发散边界预测方面, 通过计算不同动压下的静气弹收敛历程。假设动压Q当前步计算得到最大位移为un, 下一步最大变形位移为un+1, 下两步最大变形位移为un+2。通过式(1)求得耦合迭代过程位移增长比qs为

(1)

绘制各个动压Q～qs的曲线, 并利用1阶最小二乘法对数据进行处理, 求得qs=1对应的动压为发散边界。

2 静气动弹性发散计算模型

2.1 气动外形

本文计算的模型为矩形安定面, 长度为300 mm, 高度125 mm, 厚度为2.5 mm, 外形如图2所示。

图2 安定面外形Fig. 2 Configuration of the stabilizer

2.2 结构形式

安定面结构由两部分组成, 安定面和安装底座, 材料均为铝合金, 底座与安定面采用焊接形式连接, 底座前缘与安定面前缘距离为163 mm, 底座打螺纹孔, 通过螺纹安装固定。

2.3 气动力计算模型

本次计算网格采用ICEM划分, 共186个模块, 网格规模为3.5×106, 15层附面层, 附面层首层高度为2×10-5m, 附面层总厚度为1.5×10-3m。图3给出了物面网格和空间拓扑结构。

图3 物面网格Fig. 3 Surface mesh

数值模拟采用Fluent软件来求解流场, 控制方程为定常三维N-S方程, 空气状态方程按理想气体考虑, 空气黏性系数采用Sutherland公式修正, 湍流模型为标准k-ε模型, 壁面速度分布采用壁面函数修正, 计算时选用基于密度的显式格式求解。来流边界给定远场压力条件, 固定壁面给定无滑移、 无渗透边界条件。

2.4 结构计算模型

结构采用六面体纯实体模拟, 网格规模为8 000, 底座采用SPC约束六自由度。结构有限元求解采用Nastran静力分析模块求解, 图4是结构有限元模型。

图4 结构有限元模型Fig. 4 Finite element model of the structure

3 静气动弹性发散数值分析

3.1 无攻角状态计算结果

分别计算了Ma=1.5状态, 攻角0°, 60, 65, 70, 80, 85, 90 kPa这6个动压下静气弹特性, 这里挑选4种列出。在初始数值扰动下, 静气弹收敛情况如图5所示。

(a) 60 kPa

通过式(1)求得不同动压下的位移增长比, 拟合动压和位移增长比曲线见图6, 获得当前结构的静气弹发散动压为83.3 kPa。

图6 不同动压与位移增长比拟合曲线Fig. 6 Fitted curve of displacement growth ratio at different dynamic pressures

3.2 带攻角状态计算结果

分别计算了Ma=1.5状态, 攻角-2°, 65 kPa和70 kPa两个动压下的静气弹特性。在-2°攻角气动力作用下, 静气弹耦合迭代步收敛情况如图7所示。

图7 -2°攻角位移收敛曲线Fig. 7 Convergence curve of displacement at -2° angle of attack

65 kPa和70 kPa状态安定面位移都呈收敛状态, 其中65 kPa状态位移收敛值为13.2 mm, 70 kPa状态收敛值为25.7 mm, 如表1所示。

表1 两个状态位移对比

表2 发散动压对比

对比65 kPa状态和70 kPa状态, 动压增加了7.69%, 位移增长了94.7%。位移增长比例比动压增长比例大了10倍多。在实际计算状态, 70 kPa状态在耦合迭代收敛后的位移, 安定面大面区域已经超出铝合金许用应力, 见图8。

图8 70 kPa状态位移和应力云图Fig. 8 Displacement and stress contour map at 70 kPa state

4 静气动弹性发散风洞试验

4.1 风洞试验设计

本文的试验在0.6 m风洞开展, 侧壁支撑, 为了避免冲击, 采用侧壁插入机构, 如图9所示, 待流场稳定后插入均匀流场。

图9 气弹风洞试验插入机构Fig. 9 Insertion mechanism for aeroelastic wind tunnel test

试验在安定面部分粘贴应变片, 测量风洞试验过程中的应变特性, 如图10所示。

图10 风洞试验应变测量Fig. 10 Strain measurement in wind tunnel test

4.2 无攻角状态风洞试验

试验采用阶梯变动压, 0°攻角, 选取了4个动压(64.5, 77.3, 80.1, 84.2 kPa), 其中在80.1～84.2 kPa动压变化过程中安定面出现静发散破坏。试验现场照片如图11所示。

图11 0°攻角风洞试验高速摄像Fig. 11 High-speed camera recording of wind tunnel test at 0° angle of attack

试验过程中安定面有小幅度振动, 动压由80.1 kPa增压至84.2 kPa安定面突然出现破坏, 如图12所示, 破坏时刻动压为81.7 kPa。应变信号表现为随动压增加应变幅值增加, 在安定面破坏时刻, 应变突增。

图12 安定面破坏后试验照片Fig. 12 Experimental photo after destruction of the stabilizer

4.3 带攻角状态风洞试验

安定面预置-2°攻角, 如图13所示, 采用2个阶梯变动压进行试验, 实际动压由64.5 kPa到71.5 kPa安定面出现塑性变形, 如图14所示。

图13 -2°攻角风洞试验高速摄像Fig. 13 High-speed camera recording of wind tunnel test at -2° angle of attack

图14 安定面塑性变形Fig. 14 Plastic deformation of the stabilizer

5 计算与试验对比

通过计算值和风洞试验结果对比, 发现采用CFD/CSD耦合进行静气弹发散分析具有良好的数值精度。在0°攻角静发散对比方面, 偏差小于2%。

同样在-2°攻角状态, 风洞试验和计算均出现结构材料失效的现象, 表明计算与试验吻合良好。

6 结论

本文采用气动力/结构耦合的方法对一种安定面结构静气动弹性发散进行分析, 并开展了相关风洞试验研究, 得出以下结论:

1) CFD/CSD耦合可准确预测静气弹发散边界, 本文计算的安定面发散边界与试验值偏差小于2%;

2) 本文选用的安定面模型为静气弹不稳定结构, 计算和试验过程中表现为应变与动压关系呈非线性;

3) 静气动弹性不稳定结构, 弹性变形造成气动力进一步增加, 有攻角状态首先出现弹性载荷超出结构强度极限导致的破坏, 破坏时动压离发散边界仍有一定余量。