跨声速非定常空气动力计算与分析

跨声速非定常空气动力计算
Computation on Transonic Unsteady Aerodynamics
北京大学力学与工程科学系
理论与应用力学专业 00级陈雪梅
摘要
颤振问题一直是高速飞行器设计中的一大难题，特别在跨声速区段。

本文利用FLUENT6.1对一模型机翼的颤振行为进行了数值模拟，仿真机翼在高速气流中受激后扭曲变形最后发展成颤振的全过程，并对这一计算结果进行了初步分析，所得的算法具有普遍意义。

关键词：颤振，空气动力学，动网格
[引言]
早期的飞行器设计中的空气动力学分析都是将机翼﹑机身和其他气动部件当作刚体来处理。

但自第一架飞机诞生以来，空气动力学与飞机结构弹性的相互作用问题已经对航空技术的发展产生了重大影响，特别在‘彗星号’失事以后，人们对此倍加关心。

飞机在空气载荷作用下会出现可观的变形，这种变形将改变空气动载荷的分布，而它反过来又使变形发生变化。

在这种相互作用过程中，会引起振动，学术界称之为颤振。

这是一种自激振荡，它不断从气流中吸收能量。

当飞机发生颤振时，轻则出现不稳定和振动现象，重则因它引起材料‘疲劳’从而导致飞机在空中解体，以至机毁人亡。

在莱特兄弟首次试飞前，兰利的“空中旅行者”作了两次不成功的飞行试验。

第二次试飞时机翼和尾翼毁坏了，失败原因众说纷纭，气动弹性可能是第二次失败的罪魁祸首。

第一次世界大战中，英国的DH-9飞机尾翼颤振导致了飞行员死亡。

对此，英国空气动力学家贝尔斯托（L. Bairstow）首先对颤振进行了理论研究。

随着飞机速度的提高，空气动力增大，而重量小的结构形式使机翼抵抗变形的能力下降，所以气动弹性问题便得严重起来。

20世纪30年代初英国一家飞机连续发生有气动弹性引起的颤振事故，促使航空工程界对气动弹性问题普遍重视起来[摘自参考文献3，P118]。

其间的理论研究颇有成效。

美国力学家西奥多森（T. Theodorson）提交的研究报告对美国航空工业界建立颤振分析方法起了巨大作用。

50年代中后期，特别是60年代，一方面空气动力学理论的突破为非定常空气动力学研究提供了新方法；另一方面风洞技术高度发展，使振荡机翼非定常气动理论有了新的突破。

但由于理论方法的局限性以及风洞试验的高耗能及周期
长等问题，计算空气动力学应运而生。

由于涉及到非定常空气动力学，颤振及气动弹性问题的研究十分困难。

目前国内关于颤振的研究主要还是基于试验，理论仅限于线性划分析。

近年来由于计算技术的飞速发展以及CFD的实际解题能力大大扩大，用数值方法解决这样复杂的问题已是可能。

采用计算流体力学方法可缩短周期、降低费用，特别在初选阶段，优化选型需要不断改变参数、重复计算。

对那些目前不能在特定的飞行状态下进行试验的未来飞行器来说，数值模拟方法可以减少其设计风险，并在风洞实验前预先筛选设计方案。

本文的研究是与北京航空航天大学合作进行的，机翼的振型和刚度由北航提供，我们负责气动力以及气动力与机翼振动的耦合计算。

气动力计算本身已是高度非线性问题，现在还需与机翼的弹性振动相耦合，其难度相当高。

这样的研究在国内尚属首次。

在计算种我们利用了Fluent6.1软件UDF功能，这是我们的工作大大简化。

如今已完成计算的主过程，结果表明我们设计的颤振问题的直接算法是可行的，有通用价值的。

一．计算模型
1. 机翼几何外形
如图（Fig1、Fig2）所示的后掠机翼，展长为1.5m，根弦长为1.0m，梢弦长为0.6m，前缘后掠角为40︒，机翼平面位于xoy平面内，根部固支。

翼剖面为对称薄翼型。

3．网格结构
几何实体及网格用Fluent自带前处理软件Gambit2.0生成。

计算区域采用前段、翼梢段、上下为6倍展长，尾部为10倍展长的长方体区域。

由于只计算一侧机翼的情形，因此取翼根部为对称面。

流场中包含机翼的一小部分采用非结构化网格（Fig3），以便于使用Fluent中提供的动网格手段进
行计算。

其余部分采用结构化网格，以提高计算精度及加快收敛速度。

Fig3 (机翼网格图)
4. 机翼振动的数学模型
研究中假定机翼为弹性体，由于只研究小振幅下的振动，可假设振型为抛物线形式。

末梢位置呈周期性变化：
sin()
h A t ω=
其中 A 为末梢的振幅，ω为机翼周期性振动的圆频率，具体取为机翼振动的第一阶主频率
2, 4.634f f Hz ωπ==
假定机翼最大偏角为
1α=
，则振幅振幅A 具体取值为A=Lsin α，L 为
机翼展长。

按照假定，整个机翼的振型为抛物线型，于是有
2222
(,,)/sin()/z x y t y h L y A t L
ω==
其中x ，y ，z 为机翼坐标，t 为运动时间。

5. 气流主控方程
对机翼进行数值模拟的主控方程可以采用欧拉方程，也可以采用N-S 方程。

对升力和力矩的计算采用欧拉方程即可提供足够的精度，欧拉方程计算收敛速度相对较快，结果与NS 方程差别不大。

但欧拉方程无法准确计算阻力，而NS 方程对升力阻力以及力矩均能计算得很好，但收敛速度相对较慢。

本文研究机翼在跨声速段的颤振问题，8Re 10 。

因此，应该将此问题看成
充分发展的湍流问题。

采用标准的k ε-模型进行计算。

数值模拟采用的雷诺应力平均NS 方程为
为了使上述方程组封闭，采用Boussinesq 假设：
'
'2
()()3j i i i
j
t t ij
j i i
u u u u u k u u x ρμρμδ∂∂∂-=+-+∂∂∂
其中，k 为湍流动能,t μ是湍流粘滞系数，可以通过湍流动能k ,湍流耗散率ε计算：
2
t k
C μ
μρε
=
k 湍流动能,ε湍流耗散率满足的传输方程：
t k b M i k i Dk k G G Y Dt x x μρμρεσ⎡⎤⎛
⎫∂∂=+++--⎢⎥ ⎪∂∂⎝⎭⎣⎦
()2
132t k b i i D C G C C C Dt x x k k εεεεμεεεερμρσ⎡⎤⎛⎫∂∂=+++-⎢⎥ ⎪∂∂⎝⎭⎣⎦
这样就组成了封闭的方程组，即是标准k ε-湍流模式，其中的经验常数为：
121.44, 1.92,0.09, 1.0, 1.3
k C C C εεμεσσ=====
6. 网格重构算法
计算中采用了Fluent 提供的Dynamic-Mesh Modal ，机翼部分为非结构化网格，因此使用动网格中的弹性系数算法和局部重构算法对网格进行重构。

在弹性系数算法中，将任意两个网格节点之间的边等效为一根弹簧。

先由用户给定的UDF （User-Defined-Function ）函数，计算出边界上的结点位移。

此位移将在与此结点相连的任意一条边上产生一个弹性力，弹性力的大小与位移大小成正比。

由此，边界结点的位移就被传递到整个体网格。

在平衡状态下，每个结点所受的所有与它相连边上的弹性力之和为零。

该平衡条件产生了一个循环的计算结点位移的方程：
1
i
i
n m
ij j
j
m i n ij
j
k x x k
+∆=
∑∑ 其中，
i x 是结点i 的位移，i
n 是与结点i 相连的节点数目，ij
k
结点i 与相邻结
点j 之间的弹性常数，定义为：
1
ij k =
当边界结点位移已知时，就可以用Jacobi 扫描算法求解上述方程。

得到收敛解后，内部结点的位置被更新。

1,n n k converged i i i
x x x +=+ 其中，n+1和n 分别代表下一个时间步和当前时间步。

当边界结点的位移相对局部网格的尺寸很大时，
网格的质量将变得很差。

为避免这一问题，Fluent 提供局部重构算法对坏质量网格进行合并或拆分。

此时坏质量网格定义为超过某一给定体，低于某一指定体积或者网格倾斜率大于某一数值的那些网格。

二．计算过程
采用Segregated 求解器，SIMPLE 算法，利用有限体积法离散控制方程，采用一阶迎风差分格式,时间步长取为物理周期的1/20便可取得相当好的计算结果。

计算过程中监控机翼的升阻力以及力矩的变化趋势，阻力趋于稳定以及升力、力矩周期与物理振动周期吻合时结束计算。

实际计算中取来流马赫数0.8M =，边界均设为远场压力条件。

分别取攻角为
0,2,4
进行计算。

在数值模拟中的每一个计算步用UDF 函数接口给定机翼
上网格结点的位移，并利用弹性常数算法及局部重构算法对网格进行重构。

三．计算结果和分析
1.
攻角时的计算结果：
阻力系数图及升力系数图（Fig4和Fig5）
Fig4 (阻力系数曲线) Fig5 (升力系数曲线)
力矩系数图（Fig6）
Fig6（零攻角力矩系数曲线）
可以看到，阻力系数在计算稳定后保持不变，这是因为机翼只有沿垂直方向位移，对气流只产生竖直方向扰动，对水平方向的力并无太大影响。

升力与力矩的变化趋势与物理振动的变化趋势完全吻合（ 2/0.21T s πω=≈），说明了计算的合理性。

计算结果表明初始时刻升力及力矩具有最大值。

此时机翼处于平衡位置，具有最大的运动速度，即对流场有最大的扰动，理论分析上此时的升力及力矩也应取得最大值，再一次证实数值模拟方法的正确性。

2. 2 攻角计算结果：
2 攻角阻力及升力系数图（Fig7，Fig8）
Fig7 (2度攻角阻力系数曲线) Fig8 (2度攻角升力系数曲线) 2 攻角力矩系数图(Fig9)
Fig9 (2度攻角力矩系数曲线)
3.4 攻角计算结果
4 攻角阻力及升力图（Fig10,Fig11）
Fig10(4 攻角阻力系数曲线) Fig11（4 攻角升力系数曲线）
4 攻角力矩系数图（Fig12）
Fig12(4 攻角力矩系数曲线)
可以明显看到，攻角变大时，阻力系数并没有太大的变化。

而升力及力矩系数则有一个数量级以上的突变，这是由于攻角变大时对流场的扰动也相应增大。

这在前人的实践分析中也早已证明。

[参考书目3，P117，“迎角的增大正是机翼破坏的原因。

显然，……由于压力作用而顺坏]。

这就说明了计算的可靠性以及计算流体力学的可行性。

四. 研究回顾及意义展望
在早期的研究中，由于没有考虑到Fluent中数值算法与理论分析的差异性，试图用转换坐标的方法来解决问题。

计算中没有给定机翼的往复运动，而是试图设定周期性的边界条件来解决问题。

不胜遗憾，得到的计算结果与实际情况大相径庭，计算周期远小于物理运动周期，且具有相当大的随机性，即当时间步长改变时，计算所得周期也作无规律变化。

细究Fluent关于湍流模式的算法发现，其关于固壁的边界条件为一给定的插值函数，因此预定义的动边界条件在计算中被抹掉，并未起实际作用。

于是得到不合理的计算结果。

后又试图把湍流模式改为层流模式，但在如此高的Re（约为8
10的量级）下，对网格尺度的要求为1/Re（即8
10 量级），网格的数目及其奇变性又形成了另一无法调和的矛盾，仍然不能得到物理上合理的结果。

最后采用了Fluent中提供的动网格功能，直接给定机翼及其周围网格的变化规律，使模型与实际物理情况吻合。

此时问题才得到了彻底的解决，得出了合理的计算结果。

本文只是对一机翼按一具体振型作已知振动时的数值模拟的结果，它将是模拟真实机翼在流固耦合下的气动响应计算的基础。

其中对于动网格的研究，改变了以往的依靠转换坐标来完成运动边界的数值模拟的做法，得到较好的结果。

现在我们正以此为基础，考虑多振型的弹性机翼的耦合振动，即真实机翼颤振问题的数值模拟，预期不久即可有全面的结果。

致谢
本文是在陈耀松教授和李植副教授的指导下完成的。

我特别感谢陈教授，正是在他的指导、鼓励和帮助下，这篇论文才能得以完成。

陈耀松教授渊博的学术知识，敏锐的抓住问题核心的能力，简单清晰的物理建立模型的能力，“活到老学到老”的学习精神和对力学发展方向的洞察力都是我学习的榜样。

李植教授也一直关心我的研究进展，并且提供了许多好的建议。

同时，实验室里浓厚的学术气氛，轻松开放的学术氛围，都是我得以完成这一论文的必要条件。

在这一年期间，程暮林师兄给我以十分巨大的帮助，很多问题是在与他进行反复讨论的基础上完成的。

此外，实验室里的其他师兄，安亦然、李向群等，都给我以很大的地帮助。

在此我向他们以及其他所有关心和帮助过我的人表示衷心的感谢。

参考文献：
1. S.V.帕坦卡著，张政译：传热与流体流动的数值计算，科学出版社，北京，1984
2. John C. Tannehill, Dale A. Anderson, Richard H. Pletcher ：Computational Fluid Mechanics and Heat Transfer （Second Edition）, Taylor & Francis publishers, U.S.,1996
3. 李成智，李小宁：征服天空之翼—跨世纪的航空技术，湖北教育出版社，
武汉，1998
4. 管德：非定常空气动力计算，北京航空航天大学出版社，北京，1991
5. Fluent Inc.：Fluent
6.1 Tutotial Guide，U.S.，2002
6. Fluent Inc.：Fluent6.0 Dyanamic Mesh Manual，U.S.， 2001
7. E. H. 道尔，H. C. 小柯蒂斯等著，陈文俊，尹传家译：气动弹性力学现代教程，宇航出版社，北京，1991
附录：
本次数值模拟中所用到的UDF函数：
/******************************************************************** Node motion based on simple parabolic equation
********************************************************************/ #include "udf. h"
#include "dynamesh_tools. h"
#define ALPHA (1.0)
#define L (4.0)
#define A （L*sin (ALPHA/180*M_PI)
#define FRE (4.634)
#define OMEGA (2.0*M_PI*FRE)
DEFINE_GRID_MOTION (wing, domain, dt, time, dtime) {
face_t f;
Thread *tf = DT_THREAD (( Dynamic_Thread * ) dt ); int n;
Node *v;
real disp;
real y;
/*Set deforming flags*/
SET_DEFORMING_THREAD_FLAG(tf->t0);
Disp = OMEGA*A/L/L*cos(OMEGA*time)* dtime;
begin_f_loop (f, tf )
{
f_node_loop (f, tf, n)
{
v=F_NODE (f, tf, n)
if(NODE_MARK(v)==1)
{
NODE_MARK (v) = 0;
y=NODE_Y (v);
NODE_Z(v) += disp*y*y;
}
}
Update_Face_Metrics(f,tf);
}
end_f_loop(f,tf);
}
作者简介：
陈雪梅，女，1982年12 月出生于美丽的山城重庆。

2000年9月从重庆市巴蜀中学考入北京大学力学与工程科学系。

在校期间始终保持积极的学习态度，多次获三好学生和三好标兵等荣誉，并每年或本班最高奖学金，GPA总排名本班第二。

感悟与寄语：
“路漫漫其修远兮，吾将上下而求索”, 是我在研究中最大的感触。

虽然本科期间成绩一直名列前茅，但课程的学习与实际的研究是有区别的。

初进实验室时，觉得自己什么都不懂，要学的太多，经过不懈的努力与探索，现在已基本能担当起独立完成研究任务的工作。

总的说来，这是一段很有意义的经历，在其中我不仅学到了更多的理论及实践的知识，我也学会在面对困难时应该采取的态度，这些都将对我今后的研究、学习乃至生活产生不可磨灭的作用。

指导教师简介：
陈耀松，男，力学与工程科学系教授。

毕业于清华大学，曾在MIT、Cornell、Caltech等校做访问学者，目前主要从事计算流体力学的研究工作。

李植，男，力学与工程科学系副教授。

1999在莫斯科大学获得博士学位。

归国后主要从事非线性波理论的研究工作。

21。