二维定常可压缩流场分析—NACA0006翼型气动力计算

【摘要】翼型的气动力计算是空气动力学领域中十分常见的问题，本文针对一个实际NACA0006翼型，使用CFD仿真软件Fluent软件包进行气动力计算分析。计算出了在攻角为4°，马赫数为0.8的工况下的翼型外流流场，获得流场压力分布，通过FLUENT自动对翼型表面进行压力、剪切应力的积分获得在该工况下的空气动力，从而得到在工况下的翼型升力（与来流方向垂直）、阻力（与来流方向一致）和力矩。计算结果与课本正确结果完全一致。此外，改变马赫数值，比较计算结果，进一步学习了解激波和不同马赫数下的受力情况。

1．  实例简介

Fluent软件是用于模拟具有复杂外形的流体流动以及热传导的计算机程序，可以有效地模拟NACA0006翼型空气动力问题。本实例采用CFD仿真软件Fluent对NACA0006翼型空气动力问题进行了模拟分析。将对计算域进行结构化网格划分，采用压力基耦合求解器来求解可压缩流动；对来流采用远场边界条件；湍流模型选择Spalart-Allmaras模型；并且在初始化时使用Full Multigrid初始化，以获得更好的初始流场数据；除了通过残差判别收敛，还通过力和表面上设置点监视来检查解的收敛情况。

如图1所示的是NACA0006翼型，翼宽度为1m，具体几何结构数据来自课本光盘中的NACA0006.dat文件。空气自无穷远处以马赫数0.8和攻角4度绕流此翼。攻角是翼弦与来流速度之间的夹角，抬头为正，低头为负，常用符号α表示。马赫数为翼型速度与声速的比值。雷诺数本实例研究空气绕流的流动情况，分析NACA0006在马赫数0.8、攻角4度工况下的升力、阻力、和力矩。在完成该工况下的数值模拟之后，改变马赫数或攻角，分析其所造成的影响。变形一：攻角为4度，马赫数为0.6；变形二：攻角为4度，马赫数为1.2。

 http://s7/mw690/0042rYW2gy6Hyi4i3pYc6&690
图1 NACA0006翼型截面

在一个半圆形加长方形的二维平面中充满等密度的空气等介质，半圆半径为15m，长方形宽为15m。翼型以0.8马赫数的速度往左移动，同时形成空气绕流。假设实际在常温下进行，则由雷诺数 ，可判断流动为湍流状态。计算区域示意图如图2所示，为了划分比较整齐的结构化网格，将计算域划分为两个部分。               

http://s1/mw690/0042rYW2gy6HyiEPLywc0&690

图2 计算区域示意图

2．  数值模拟方法

http://s3/mw690/0042rYW2gy6HyiI60rE02&690

    对于不同的具体问题都有其适合的流动模型。Fluent提供的模型包括：单方程（spalart-allmaras）模型、双方程（k-ε）模型系列、（k-ω）模型系列、雷诺应力模型等。其中单方程（spalart-allmaras）模型是一个相对简单的一方程模型，求解了一个有关涡粘性的输运方程，适合具有壁面限制的流动方程。例如飞行器、翼型等绕流流场分析，所以本实例选用单方程（spalart-allmaras）模型。

2.2 边界条件及网格划分

     翼型模型可以从现有的附书光盘中导入，也可以自己用GAMBIT绘图的方法自己绘制。绘制方式有从上到下的（直接建面和直接建体）和从下到上的（点—线—面—体），绘图过程中应该注意将关键的点面体进行有意义的命名，可以方便后期网格划分。也可以将PRO/E等建模软件所建模型导入到GAMBIT在进行网格划分。

网格划分全部采用四边形的结构化网格（见图3和图4），左半边网格总数为300 300，右边网格总数200 100。网格划分采用接近翼型处密，远处稀的方式。运用FLUENT软件进行数值模拟，经过迭代计算，各个方程的残差均小于10-4，计算至结果收敛。

  http://s8/mw690/0042rYW2gy6HyiRZ8ebb7&690
     图 3 计算区域网格示意图                 图4 NACA0006翼型表面附近网格示意图

翼型给定的几何形状为梭行和尺寸为长1m，宽0.06m。只需设定两种边界条件：计算域外围和翼型上下表面的边界条件。将最外圈边界设为远场边界条件（pressure_far_field），该边界条件只对可压缩流动适合。对翼型的上下表面设为避免条件（wall）。本例为定常问题，无需设置初始条件。

3. 结果与讨论

    通过FLUENT求解，观察图5-9 残差动态显示、监控点piont 4 的速度、翼型的阻力系数、翼型的升力系数、翼型的力矩系数经过不到200次均达到收敛。这表明，网格划分、边界条件的设置、算法的选择等符合要求，计算结果有效。

http://s12/mw690/0042rYW2gy6HyiVBcSL9b&690

图 5残差动态显示图

  http://s6/mw690/0042rYW2gy6Hyj0xkr375&690

  图 6监控点piont 4 的速度随迭代过程变化曲线图      图 7翼型的阻力系数随迭代过程变化曲线

http://s12/mw690/0042rYW2gy6Hyj0F8bF4b&690

  图 8 翼型的升力系数随迭代过程变化曲线           图 9 翼型的力矩系数随迭代过程变化曲线

从图10可以看出翼的压力分布大致如下：压力最大区域在翼的下表面的前部，其次是翼的下表面后半部分及上表面的后半部分，压力最小的部分在上表面的前半部分。清楚地表明了翼型的升力阻力的位置及大小。上表面中间靠后位置出现了一道明显的激波，由激波性质可知，穿过激波，压力上升、温度上升、速度下降、马赫数从大于1减小到小于1。本实例中清晰的捕捉到了激波，壁面激波附近的网格和压力云图分布（图11-12）表明，压力从蓝色区域的较低值在激波位置瞬间提高到绿色区域较高的值。

http://s8/mw690/0042rYW2gy6Hyj8wm2j07&690

 图 10翼周围压力分布云图                     图 11 壁面激波附近的网格和压力云图分布

http://s15/mw690/0042rYW2gy6Hyj8BFgOae&690

图 12翼周围压力分布等值线图

图13-14表明速度分布如下，在翼型上表面前半部分红色区域速度最大，翼型上表面附近后半部分蓝色区域速度最小。其余部分呈黄绿色速度值在这两者之间。从速度分布等值线图（图13），可以发现在激波位置最靠近翼型上表面处，速度由最大突然降为最小，这与激波性质相符，穿过激波，压力上升、温度上升、速度下降、马赫数从大于1减小到小于1。表明此处确实出现了激波。

   http://s7/mw690/0042rYW2gy6HyjfN4Ymb6&690
   图 13激波位置靠近翼型上表面处速度分布等值线图              图 14马赫数分布云图

由翼型上下表面的压力系数分布线图（图15）可知，红色的表示翼型上表面压力系数，白色的表示压力下表面压力系数。可以发现翼型下表面压力稳定，没有较大变化，没有形成激波，而上表面形成了激波，而且激波出现在0.45倍弦长处。

 http://s10/mw690/0042rYW2gy6HyjjrI1z09&690
  图 15 翼型上、下表面压力系数分布线图            图 16翼型上下表面切应力的X分量分布曲线图

    由于激波引起的逆梯度越大，边界层的分离就越严重，分离点就是壁面切应力消失的点，逆向流动可以通过壁面切应力的X分量是否为负来判断，如图16，可以看出边界层的分离点大约出现在0.45倍弦长的位置，这与图15表明的激波出现在0.45倍弦长处结论相吻合。从翼型附近的速度矢量图（图17）可以清楚的看到边界层分离情况，分离位置也是刚好出现在0.45倍弦长的位置。同时可以发现激波后方翼型上表面附近有一个明显的回流区。

 http://s10/mw690/0042rYW2gy6Hyjly7IBc9&690

图 17 翼型上壁面附近的速度矢量图

4．  改变马赫数对激波的影响

在前文工况下，已经得到了激波，也符合激波的性质。下面通过只改变马赫数，来学习不同马赫数对激波产生的影响，以及对翼型空气动力的影响。

4.1 攻角为4度，马赫数为0.6

从图18-20可以看出翼的压力分布大致如下：压力最大区域在翼的下表面的前部，但是整个下表面压力变化不大，压力最小的部分在上表面的前半部分，清楚地表明了翼型的升力阻力的位置及大小。上下表面压力都是连续变化的，翼型上、下表面压力系数分布线图没有出现类似图15中的断裂。也可以看出上下表面间的压力差，没有0.8马赫数时的大，表明升力有所下降。通过与前文图10-11的对比发现，没有产生激波，上下表面压力变化梯度都不是很大，没有像穿过激波后的压力瞬间上升的现象。进一步从反面验证了激波的性质。

 http://s2/mw690/0042rYW2gy6HyjmVumJ31&690
    图18 翼周围压力分布云图                         图19翼型上、下表面压力系数分布线图

  http://s6/mw690/0042rYW2gy6HyjpI5zn75&690   http://s10/mw690/0042rYW2gy6HyjpNliN89&690
     图20翼型上下表面切应力的X分量分布曲线图                图21 马赫数分布图

图21表明速度分布如下，在翼型上表面前半部分深橘色区域速度最大，翼型前端尖尖处淡蓝色区域速度最小。大部分呈黄绿色速度值在这两者之间。上下表面的速度差没有0.8马赫数工况下大，表明升力降低，这与图18的结论一致。攻角不变，马赫数减小，升力也随之减小。

4.2 攻角为4度，马赫数为1.2

    由图22-23可以清楚地看到，在翼型尾部有一部分压力特别低的区域，其之后压力又瞬间升高，这表明此处产生了激波，而且在马赫数增大的情况下，激波会往后移。而且在此工况下，翼型尾部的上下表面压力差最大，这种情况下可能会使飞机尾部向上翘，头部往下低，导致飞机飞行不平衡。图中翼型头部的压力最大，说明其他条件保持不变，增加马赫数，飞机受到的阻力会明显增大。

  http://s10/mw690/0042rYW2gy6Hyjt6CiR69&690
       图22 翼周围压力分布云图                                图23 马赫数分布图

    从图24-25可以看出，翼型上下表面附近的压力变化和切应力变化都不大，但是从图24可以看出，在翼型末尾点的位置，上表面的压力有一个明显的瞬间升高（红色线），从翼型上下表面切应力X分量分布曲线图中可以看出，在翼型末尾点的位置切应力明显呈下降趋势，这与图15-16中的出现激波的位置处的上表面X分量切应力变化相同，这表明此处产生了激波，与图22-23的结论相符。

 

http://s7/mw690/0042rYW2gy6HyjywFH8d6&690
  图24翼型上、下表面压力系数分布线图               图25翼型上下表面切应力的X分量分布曲线图

5．  结论

本实例是二维定常可压缩流场分析，对NACA0006翼型进行了气动力计算，所得计算结果与课本（fluent入门与进阶教程）完全一致。计算获得了流场的压力分布，算出了在4度攻角，马赫数为0.8的工况下的空气动气（升力、阻力、力矩）。计算过程中采用监测点监测求解的收敛性，并且采用多种后处理工具检查与激波相关的流动现象，发现激波出现在翼型上表面0.45倍弦长处，同时也验证的激波性质（穿过激波，压力上升、温度上升、速度下降、马赫数从大于1减小到小于1）的正确性，此外还通过改变马赫数值，学习受力的变化以及激波的变化。发现当马赫数小于一定值时，不会产生激波。当产生激波以后，马赫数越大，激波位置会向后移，阻力会逐渐增大。