当前位置：首页 > 涨幅分析 >飞机机翼三维湍流流场数值模拟

飞机机翼三维湍流流场数值模拟

2023-05-10 14:56:27

飞机机翼三维湍流流场数值模拟

黄振(山东豪迈化工技术)

引言

长期以来地面的风洞试验都是测量飞机气动力数据和试验新的气动力现象以及概念的主要手段，通过反复的设计和实验，并借助经验的判断确定最终的设计方案，所以对经验的依赖性较强，设计的周期较长，费用较高。随着计算机技术和数值模拟软件技术的发展，数值模拟开始广泛应用于机翼的设计和流场分析中，近十几年来的研究已显示其特有的能力和潜力，因而受到 NASA (美)， AGARD(欧)，NAL (日) 等研究机构的普遍重视。一直以来，国内外对翼型仍在进行相当多的研究^[1]，本文对某型飞机机翼的三维外部流场进行模拟，分析了机翼的气动特性和三维机翼外部流场的一般规律。

控制方程

完整的雷诺应力模型包括平均运动的1个连续性方程，3个动量方程。雷诺应力的6个方程，k方程和ε方程，总共由包含12个未知量的12个微分方程组成封闭的方程组^[1^,2]，如果要涉及传热，则需要计算温度和其他有关标量的分布，还要加上1个平均温度方程和其它3个与扰动量相关的关联项的模型方程，总共16个方程，对于工程中的湍流问题，这个方程组实在是太庞大了，对于复杂流动，若要求解这个方程组，需要占用极大的计算机资源，所以在工程应用中一般对雷诺应力模型作一些简化。

1.1 三维流动的 Navier-Stokes 方程

动量平衡定律也是任何流动系统都必须满足的基本定律。该定律可表述为：微元体中流体的动量对时间的变化率等于外界作用在该微元体上的各种力之和。该定律实际上是牛顿第二定律。按照这个定律，可导出 x, y和z 3 个方向的动量平衡方程( Momentum conservation equation) ，也称为Navier-Stokes 方程。

若微元体受到的力只有重力，且 Z 轴竖直向上，则Fz=0, Fy=0,Fz=-ρg。

1.2 湍流模型及方程

目前两方程模型中使用最为广泛、最基本的两方程模型是标准κ-ε模型，即分别引入关于湍动能k和耗散率ε的方程。但是对于强旋流，弯曲壁面流动或弯曲流线流动时，会产生一定的失真，原因是在标准模型中，对于Reynolds应力的各个分量，假定粘性系数是相同的，即假定μ_t是各向同性的标量，而在弯曲流流的情况下，湍流是各向异性的，μ_t应该是各向异性的张量^[3]。为了弥补标准κ-ε模型的缺陷，许多研究学者提出了对标准κ-ε模型的修正方案，本文使用了Realizable κ-ε 模型。

在Realizable κ-ε 模型中，湍动粘度计算式中的系数Cu不是常数，而与应变率联系起来。Realizable κ-ε 模型中关于κ和ε的输运方程如下：

其中，

公式 (6) 中， u_t与 Cu按下式计算：

Realizable 模型已被有效地用于各种不同类型的流动模拟，包括旋转均匀剪切流、包含有射流和混合流的自由流动、管道内流动、边界层流动，以及带有分离的流动等。但这种标准的模型并不依赖于分子粘性，在模拟近壁面的流动时会失效，但是近壁面的流动特性对机翼的流动特性，特别是对阻力的影响很大，因此在求解近壁面流动时使用了加强壁面函数法。

1.3加强壁面函数法

对于有固体壁面的充分发展的湍流流动，沿壁面法线的不同距离上，可将流动划分为壁面区和核心区，对于核心区的流动是完全的湍流区，在壁面区流体的流动受壁面的影响比较明显，在壁面区上到壁面的距离决定了粘性力的作用大小^[3,4]。加强的壁面函数法就是综合考虑近壁面流动的特性，形成一种求解整个壁面区域流动特性的方法，这种方法是综合求解层流和湍流的壁面流动的方法。

其中：

同样，可以得到：

在公式中分别求解层流和湍流的情况：

其中 y⁺_s为对数律层的分界线，一般的 y⁺_s=60，系数α表示压力梯度的影响，β， γ表示热传导的影响，在本文研究过程中不考虑传热，所以β，γ的大小为0。而对于层流部分：

对于层流部分只考虑压力梯度的影响，热传导和压缩性的影响忽略不计^[5]。

机翼三维湍流数值模拟

对某型无人机机翼三维湍流流场进行数值模拟，将三维无人机机翼模型放入一个外部流场中，在机翼的前端设置速度入口，后端设置压力出口，机翼的根部设为对称面，翼梢外部延伸部分流场。

2.1问题描述

三维机翼在外部气体流场中的模型如图1所示，机翼的根部设为对称面，机翼的前端15倍弦长的位置为外部流场的入口，机翼的后端20倍弦长处为气流出口，翼梢外部有2倍弦长的外部流场。

图 1 三维机翼模型

2.2网格生成

对于复杂三维的飞行器外形 (或其他复杂流动外形)，生成单块 (贴体) 的计算网格很困难，即使勉强生成，网格质量也不能保证，还影响流场数值求解的效果。因此，目前常用分区网格及分区计算技术，即根据外形特点将总体流场分成若干个子区域，对每个子域分别建立网格。区域分解的基本原则是：尽量使每个子域的边界简单以便于网格的建立，同时尽量实现在关键流动区域重点网格，三维机翼部分网格如图2。对于复杂外型物体，表面网格和拓扑结构很重要。表面网格是塑造模型形状的重要体现；拓扑结构则决定如何正确分区，合理规划网格。分区网格生成步骤通常可归结为：

(1) 根据外形和流动特点分区，并确定每区中的拓扑结构。

(2) 生成表面网格。表面网格的疏密有时不仅仅要考虑流场变化情况，还要考虑模型的外型复杂度，及曲线 (面) 的曲率问题。外形复杂或曲率大的地方就要适当的密一些，这样才可以较准确地反映出模型的形状。

图 2 机翼外部流场的部分网格图

2.3结果分析

计算状态取马赫数为0.1，雷诺数1.0×10⁶。计算参考面积为机翼法向投影面积，参考长度为翼根弦长，湍流模型的选择和离散方法，先用一方程Spalart-Allmaras模型，离散方程采用一阶迎风格式进行计算，再用二方程Realizable 模型，离散方程采用二阶迎风格式，在求解近壁面流动时使用加强壁面函数法进行计算，直到计算收敛。图3为三维计算模拟的残差曲线图。

图 3 三维模型在计算过程中残差的变化曲线

（攻角为0°。速度为33 m/s）

图3中前1100多步是一方程Spalart-Allmaras模型，离散方程采用一阶迎风格式进行计算，图中后半部分为二方程Realizable模型，离散方程采用二阶迎风格式进行计算。从图中可以看出一方程计算时，残差能够很快收敛，但是收敛到一定程度后就很难继续收敛，一直在小幅振荡，换成二方程计算后，残差收敛的速度慢，但是精度要高，经过长时间的迭代后，残差能大幅下降。

图4图5给出了马赫数为0.1时的三维机翼气动特性曲线。

图 4 翼型升力系数曲线

图 5 翼型阻力系数曲线

对比升力曲线可以发现，三维模拟的结果基本上略小于实验数据，但是比较接近，特别是在小攻角情况下计算结果的准确度比较高，当攻角较大时，计算的误差开始有增大的趋势。这可能是因为计算网格较复杂，在大攻角下翼型流场接近失速区，对大攻角复杂流动的模拟有偏差的缘故。N-S 方程的阻力计算问题一直是个难题，计算精度往往不高^[6]，计算的结果基本上在实验值的 2 倍左右，这与其它的模拟软件的计算结果相似。对比阻力曲线可以发现，模拟计算结果与实验结果的变化趋势基本一致，能够较好地反映三维机翼的阻力变化特性，但是在大攻角接近失速区时，阻力误差开始增大。

机翼表面的压力分布如图 6 所示，可以看出上表面的压力等值线并不是平行的，靠近翼梢的部分压力受绕流的影响。

图 6 机翼表面的压强分布等值图

沿机翼的展向分别在 z/c=0.01，z/c=0.48，z/c=0.9，z/c=0.98 取机翼剖面，比较不同区域机翼上下表面的静态压力分布。图 7、 8、 9、 10 分别是机翼展向的表面压强分布曲线图。比较 4 幅图可以发现，机翼剖面越接近于翼梢，剖面压强系数差越小，上下表面的压力差也越来越小，造成了升力的损失，这对机翼的升力系数有直接影响，特别是对于小展弦比机翼的升力系数影响会比较大。这与理论相符，由于机翼三维效应的作用，在翼尖效应影响下，翼尖附近剖面的压强差减小。

结束语

(1)利用Fluent软件对无人机机翼的三维数值模拟可以直观地研究三维机翼外部流场的气动特性，从图中可以直观看出机翼外部气流流动特性。求解N-S 方程，通过对流场特性的分析，结果合乎理论分析与实际流场特征。通过与实验对比发现，本文所得结果准确可信，证明本文方法可以用来进行无人机机翼流场数值模拟计算。