基于GKS理论与笛卡尔网格的旋翼流场数值研究

Numerical Simulation of Rotor Flowfield Based on GKS Theory and Cartesian Grid

摘要: 直升机的旋翼气动力模型之所以比固定翼的复杂，体现在其许多参数都是非定常的，如直升机的旋翼升力和相对气流的速度大小、方向以及桨叶的切面迎角密切相关，即使直升机在定常飞行状态下，这些因素也是随时间发生变化的。现有的欧拉方程法在模拟运动的旋翼运动时，受到网格及其变形等的诸多限制，并不能精确地模拟旋翼的实时运动，因此，本文依据Boltzmann方程，运用GKS理论，采用笛卡尔自适应网格技术来实现翼型位置的实时变化，对NASA12翼型进行运动分析，最终获得旋翼的模拟结果。通过分析模拟结果，即翼型的宏观特性，翼型的升力以及旋涡等的分布，并将数值模拟结果和经典的试验结果进行对比。对比结果表明，整体模拟结果和Caradonna-Tung试验结果误差很小，因此，本文的方法能够较为精确地模拟直升机翼型的运动。

Abstract: The reason why the rotor aerodynamic model of helicopter is more complex than that of fixed wing is that many of its parameters are unsteady. For example, the rotor lift of helicopter is closely related to the velocity and direction of the relative air flow and the section angle of attack of blade. Even if the helicopter is steady, these factors change with time. When the existing Euler equations simulate the moving rotor motion, it receives many restrictions such as grid and its deformation, and cannot accurately simulate the motion of the rotor in time. Therefore, based on Boltzmann equation, GKS theory and Cartesian adaptive grid technology, this paper realizes the real-time change of airfoil position and analyzes the motion of the wing nasa12. Finally, the simulation results of rotor are obtained. The numerical simulation results are compared with the classical test results, such as the characteristics of the airfoil, the lift of the airfoil and the distribution of vortices. The comparison results show that the errors between the overall simulation results and caradonna Tung test results are few. Therefore, the method in this paper can accurately simulate the motion of helicopter airfoil.

1. 引言

直升机飞行动态特性不仅决定于飞行状态，更主要的是旋翼升力面。与固定翼飞机不同，直升机旋翼结构复杂，在飞行时其不单有绕桨毂的转动，还有桨叶的挥舞运动。直升机前飞时，由于旋翼的周向来流速度不对称而引起整个旋翼锥体后倾，从而使整个旋翼升力的方向发生了变化；在考虑旋翼与机身的耦合作用时，旋翼除了给机身一个反扭矩之外，当旋翼这个陀螺的自转轴发生进动时，其对机身还会产生陀螺效应。因此，旋翼和整机的动态数学模型必须考虑桨叶的气动弹性、动态特性、旋翼和机身之间气动干扰、接近地面时地面效应、桨毂结构形式、桨叶挥舞和摆振等 [1] [2] [3]。

基于计算流体动力学的方法分析直升机旋翼的空气动力特性，按照其发展共经历了跨音速小扰动速势方程法、全速势方程法、欧拉方程法和Navier-Stokes (N-S)方程法，CFD方法不但考虑到气流的粘性，还考虑了直升机高速前飞时，由于旋翼的周向来流速度不对称而引起前行桨尖的激波现象和后行桨尖的动态失速现象 [4]。

描述气体的运动可以有两种途径。第一种是从宏观角度出发，利用守恒定律，结合流体的应力–应变、热流–温度梯度等本构关系，建立宏观物理量所满足的控制方程，如Euler方程和N-S方程；第二种是从微观气体运动出发，结合气体分子或粒子间的碰撞模型，建立速度分布函数所遵循的控制方程，如Boltzmann方程 [5] [6]。后者就是统计力学，或者气体动理论的研究范畴。两种描述途径之间可以通过对分布函数在相空间的积分，并借助Chapman-Enskog展开联系起来。

高精度的、NS无法计算的瞬态气动力(如大攻角失速状态下大分离流动的气动力计算问题)、稀薄气体气动力等问题和地面交通行业的全细节发动机舱散热问题、气动和水动力噪声问题，从微观气体运动出发的Boltzmann方程从物理层面进行求解，精度更高。

本文依据Boltzmann方程，运用GKS理论，获得翼型的宏观特性，最后将数值模拟结果和试验结果进行对比。本文的方法能够较为精确地模拟直升机翼型的运动。

2. 数值方法

2.1. GKS理论

描述气体的运动可以有两种途径。第一种是从宏观角度出发，利用守恒定律，结合流体的应力–应变、热流–温度梯度等本构关系，建立宏观物理量所满足的控制方程，如Euler方程和N-S方程；第二种是从微观气体运动出发，结合气体分子或粒子间的碰撞模型，建立速度分布函数所遵循的控制方程，如Boltzmann方程。后者就是统计力学，或者气体动理论的研究范畴。两种描述途径之间可以通过对分布函数在相空间的积分，并借助Chapman-Enskog展开联系起来。

由于Boltzmann方程的复杂性，实际中用得最多的是由Bhatnagar，Gross和Krook于1954年提出的BGK模型方程。其核心思想是速度分布函数在平衡态附近回到平衡态的速率与偏离平衡态的大小成正比。半个多世纪以来的研究表明，BGK模型的应用范围已经远远超出了Boltzman方程本身，如可以不限于稀疏气体假设。这得益于模型中包含的松弛过程在很多领域和真实的物理过程有着紧密的联系。

在Gas-kinetic的气体动力学描述中，所有相关的流动变量被定义为分布函数的矩。令ρ为密度，(U, V, W)为笛卡尔坐标中的速度矢量，E为总能量。

守恒量q，由下式可得：

$q={\left(\rho ,\rho U,\rho V,\rho W,E\right)}^{\text{T}}$ (1)

可以被写成

$q\left(x,y,z,t\right)={\displaystyle \int \phi f\left(x,y,z,u,v,w,t,\xi \right)\text{d}u\text{d}v\text{d}w\text{d}\xi }$ (2)

其中f是分布函数，u，v，w是相空间速度，n表示内部自由度。在下文中，为了简单起见，我们限制为单原子气体。对于多原子气体，内部自由度在形式上只会造成微小的差异。对于单原子气体，向量φ可以写成：

$\phi ={\left(1,u,v,w,\frac{1}{2}\left(\stackrel{\to }{u}\cdot \stackrel{\to }{u}\right)\right)}^{\text{T}}$ (3)

其中， $\stackrel{\to }{u}={\left(u,v,w\right)}^{\text{T}}$。

分布函数的推演由玻尔兹曼方程给出：

${\partial }_{t}f+\stackrel{\to }{u}\cdot \nabla \stackrel{\to }{f}=Q\left(f,f\right)$ (4)

这里 $Q\left(f,f\right)$ 为碰撞积分，Bhatnagar等人首先介绍了用简单松弛项替换碰撞积分的概念，从而得到BGK方程：

${\partial }_{t}f+\stackrel{\to }{u}\cdot \nabla \stackrel{\to }{f}=-\frac{f-g}{\tau }$ (5)

这里g是由麦克斯韦分布给出的均衡分布：

$g\left(\stackrel{\to }{x},\stackrel{\to }{u}\right)=f_0=\rho \left(\stackrel{\to }{x}\right){\left(\frac{\lambda \left(\stackrel{\to }{x}\right)}{\pi }\right)}^{3/2}{\text{e}}^{-\lambda \left(\stackrel{\to }{x}\right)\left(\left(\stackrel{\to }{U}-\stackrel{\to }{u}\right)·\left(\stackrel{\to }{U}-\stackrel{\to }{u}\right)\right)}$ (6)

其中：

$\lambda =\frac{\rho }{2p}$

p是压力。我们更倾向使用符号 $f_0$ 来表示均匀分布，以强调这个函数实际上是Boltzmann和BGK方程的Chapman-Enskog扩展中的第零阶项。在粗略的合理化中，BGK方程可以被认为是对分布函数的建模，且基于假设流动在一个时间尺度，即碰撞时间上达到局部平衡。

从Chapman-Enskog扩展可以发现，Navier-Stokes方程可以从BGK方程中恢复，粘度系数为：

$\mu =\tau p$ (7)

将碰撞时间和分子碰撞过程直接关联到宏观动力粘度。

我们有兴趣在有限体积方法的框架下利用BGK方程。在本文中使用Gas-Kinetic方法背后的关键思想是，在邻近单元格中的重建分布函数的离散控制体之间的界面处计算通量。考虑控制体 $V_i$ 的积分形式的Navier-Stokes方程：

${\displaystyle {\int }_{V_i}\frac{\text{d}q}{\text{d}t}\text{d}V}+{\displaystyle {\int }_{\partial V_i}\stackrel{\to }{F}\cdot \text{d}\stackrel{\to }{A}}=0$ (8)

F(q)是Navier-Stoke的通量矢量。如果假设局部方向x和对应的相位空间速度u与面(可以在每个面处的旋转坐标上工作)垂直，则可以获得下式：

$\stackrel{\to }{F}\cdot \text{d}\stackrel{\to }{A}={\left(F_{\rho },F_{\rho u},F_{\rho v},F_{\rho w},F_E\right)}^{\text{T}}={\displaystyle \int u\phi f\left({\stackrel{\to }{x}}_l,\stackrel{\to }{u},t\right)\text{d}\Xi \text{d}A}$ (9)

上式中， $\Xi =\text{d}u\text{d}v\text{d}w$ 表示速度空间、 ${\stackrel{\to }{x}}_l={\left(x_l,y_l,z_l\right)}^{\text{T}}$ 为面坐标，推到上式过程中，我们利用了这一事实：玻尔兹曼和BGK方程的右边含有 $\phi$ 的矩消失了。

以这种方式计算的通量将是有效的，只要重组分布函数与精确分布函数一致，至少对于欧拉或Navier-Stokes的精度阶数，取决于求解哪个方程。因此，Gas-Kinetic Scheme的基本任务是构建这样的分布函数。然后可以根据等式计算数值通量，最后可以及时整合以推进数值解。

2.2. 基于笛卡尔网格的自适应加密

本文采用的是基于Cut Cell的笛卡尔网格生成算法，原理类似于Cart3D，在Cut Cell计算中，采用Cohen和Sutherland-Hodgeman剪裁算法的混合方法来求任意几何曲面与立方体的交集。

同时本网格能够进行网格自适应加密。

强间断、边界层以及高温非平衡效应是高超声速流动中比较复杂的几个问题，当这些问题耦合一起进行考虑时，流场在空间和时间上都会出现多尺度化效应，即不同的流场区域需要不同的网格尺度才能保证得到准确的结果。在计算中常采用网格自适应加密方法 [7]，在尺度小的局部区域进行网格加密，在大尺度区域进行网格稀疏，从而提高计算精度，同时提高计算效率。此外，高温非平衡流动的数值计算比较耗时，发展并行化的计算程序，是高超声速数值计算提高计算效率更为有效的途径。

在初始网格的基础上，根据流场剧烈变化程度对网格进行局部加密和稀疏，采用流场物理量的截断误差反映变化的剧烈程度。取流场密度的截断误差，其计算公式如下：

${ϵ}_{\rho }=\max \left(\frac{\left|\stackrel{\to }{l_{ij}}\cdot \left[{\left(\stackrel{\to }{\nabla }\rho \right)}_c-{\left(\stackrel{\to }{\nabla }\rho \right)}_i\right]\right|}{{\alpha }_f{\rho }_c+\left|\stackrel{\to }{l_{ij}}\cdot {\left(\stackrel{\to }{\nabla }\rho \right)}_i\right|},\frac{\left|\stackrel{\to }{l_{ij}}\cdot \left[{\left(\stackrel{\to }{\nabla }\rho \right)}_c-{\left(\stackrel{\to }{\nabla }\rho \right)}_j\right]\right|}{{\alpha }_f{\rho }_c+\left|\stackrel{\to }{l_{ij}}\cdot {\left(\stackrel{\to }{\nabla }\rho \right)}_j\right|}\right)$ (10)

其中， $\stackrel{\to }{l_{ij}}$ 是网格单元i和j之间的方向向量，见图1。图中网格i的所有四个边都利用上式计算截断误差，然后取最大值作为网格j的密度截断误差。 ${\alpha }_f$ 是防止分母出现为零情形的一个小量。

Figure 1. The schematic diagram of the interface AB and its adjacent grids i and j

图1. 界面AB以及相邻网格i和j的示意图

一般来说，取密度的截断误差能够表征流场变化是否剧烈，但在有些区域密度变化较小，而速度变化剧烈，只取密度这一个变量无法保证在速度变化剧烈的区域能够进行网格加密，所以程序中将密度 $\rho$ 和速度u、v的截断误差最大值作为衡量流场剧烈变化程度的特征量。

$ϵ=\max \left({ϵ}_{\rho },{ϵ}_u,{ϵ}_v\right)$ (11)

对于网格单元i，加密和稀疏准则如下：

当 $ϵ>{ϵ}_r$ 时，加密此网格；

当 $ϵ<{ϵ}_c$ 时，稀疏此网格。

${ϵ}_r$ 和 ${ϵ}_c$ 为网格加密和稀疏的判定阀值，一般取 ${ϵ}_r=0.08$， ${ϵ}_c=0.05$， ${\alpha }_f=0.03$。

2.3. 湍流模型

大涡模拟所采用的模型一般都是Smagorisnky给出的最简单的模型。它需要采用比较小的网格，因而也需要比较大的计算机。西德学者采用了一个方程和二个方程的模型，因而可以计算一些更为复杂一些的湍流运动(如圆环中的湍流、后台阶的湍流等)。LES Smagorisnky模型的公式为：

$\frac{\partial }{\partial t}\left(u_i\right)+u_j\frac{\partial }{\partial x_i}\left(u_i\right)=-\frac{\partial }{\partial x_i}\left(\frac{p}{\rho }+\frac{\tau }{3}\right)+\frac{\partial }{\partial x_i}\left[\left(v+v_i\right)\left(\frac{\partial u_i}{\partial x_i}+\frac{\partial u_j}{\partial x_i}\right)\right]$ (12)

其中， $\frac{\partial u_i}{\partial x_i}=0$。

与RANS模型相比，大涡数值模拟的亚格子模型具有较大的普适性。湍流大涡数值模拟方法中需要封闭的量是亚格子应力，它和大尺度脉动的相关微弱。亚格子应力是不可解小尺度脉动和可解尺度之间的动量交换，它和强烈依赖于流动边界的大尺度脉动相关性很小，因此合理的亚格子模型将有较大的普适性。LES Smagorisnky模型可以获得流动的动态特性，而RANS模型只能提供定常的气动力特性。LES Smagorisnky模型的解包含大于过滤尺度的所有脉动，由此可以获得速度谱以及气动力谱等，这些动态气动力特性对于分析高超声速导弹飞行工况下高温非平衡效应是十分重要的。

2.4. 边界条件设定和网格

该算例翼型采用NACA 0012，见图2，不存在扭转。本案例采用粗网格对Caradonna-Tung旋翼案例进行了稳定性测试，如图3所示。本案例中，旋翼展弦比为6，弦长0.625，桨尖马赫数0.439，变距角8˚。最小网格尺寸0.02 ft，网格数19467761，旋翼旋转一周的周期为0.048 s，共计算0.15 s，约3.1圈。

Figure 2. The calculation area and its grid of the NACA0012 airfoil

图2. 计算区域和网格和 NACA0012翼型

Figure 3. Cartesian mesh for airfoil adaptation

图3. 翼型自适应的笛卡尔网格

Table 1. The parameters of the wing

表1. 翼型参数

3. 结果与分析

如表1所示，本数值模拟结果和Caradonna-Tung的试验结果进行对比。Caradonna和Tung发表了一组旋翼压力测量的实验结果，该实验在NASA Ames研究中心完成，测量对象为高速旋转运动的旋翼。由于该实验针对不同的总距及旋转速度(翼尖马赫数)有详尽的实验数据，如今已成为考核各种旋翼数值计算方法的标准算例。

(a) (b) (c)

Figure 4. Comparison between the upper and the lower surface lift of the sections at 0.8, 0.89 and 0.96 stations and the experiment. (a) 0.8 station section; (b) 0.89 station section and (c) 0.96 station section

图4. 0.8、0.89、0.96站位截面上下表面升力与实验对比。(a) 0.8站位截面；(b) 0.89站位截面；(c) 0.96站位截面)

图4为0.8、0.89、0.96站位截面上下表面升力与实验对比。横坐标为弦长，纵坐标为升力。

从0.8站位截面处，整体趋势与实验结果一致，尤其是0.1弦长后，旋翼升力基本稳定，现有的模拟结果可以精确捕捉稳定后的升力。略有区别的是0.05弦长前后，旋翼升力从零迅速上升到最高点，在其最高点和下降的过程中，现有数值模拟结果与实验结果略有误差。但是整体来看，现有的GKS理论能够较为精确的捕捉旋翼升力，在其快速上升到最高点时和迅速下降过程中误差略微偏大一点。

从0.89站位截面处的数值模拟结果来分析，整体趋势与实验结果一致，尤其是吸力面的升力，从零到最大，其数值模拟的结果和实验结果非常接近，整体模拟精确度很高。但是略有不足的是在升力面，数值模拟结果并没有捕捉到升力面上的两个极小值，数值模拟的结果和实验结果的误差交大。但是升力面后面的模拟数据和实验数据非常接近，误差很小。总之，在0.89站位截面处，整体的吸力面和压力面的升力模拟都和实验模拟结果很接近，误差极小，除了在压力面的极小值处，数值模拟结果和实验结果误差交大。

从0.96站位截面处的数值模拟结果来分析，与0.89站位截面处的数值模拟结果较为类似，整体的吸力面和压力面的升力数值模拟结果和实验结果都较为接近，除却压力面的最高升力的小值处误差交大。

综合分析，整体可以看出，可以看到在0.1弦长前，即前两圈，旋翼在相对地面静止的流场中转动时受到较高的升力；在0.1弦长后，即第三圈开始后，旋翼升力基本稳定。整体的数值模拟结果和实验结果非常接近，误差很小。个别的升力极高值附近，尤其是压力面上，存在模拟的精度不足的情况。

如图5所示为不同时刻的速度分布图和尾迹涡图。从右图可以看出，随着旋翼的运动，形成一系列的尾迹涡带，在不同的时刻具有不同的分布。

Figure 5. The velocity distribution

图5. 速度分布图

4. 结论

本文采用GKS理论对旋翼升力进行了非定常的数值模拟，模拟网格采用笛卡尔自适应网格。旋翼采用Caradonna-tung的NACA0012翼型。经过数值模拟结果和实验结果进行对比发现，数值模拟能够较为精确地捕捉压力面和吸力面的升力变化，整体误差较小，精确度很高。误差最大的地方是处于升力变化的极高值上，数值模拟结果相对实验结果较高。初步分析，部分原因是由于笛卡尔网格的关系，由于在边界层处反应边界的形状不够精确造成的。在后面的研究中，希望能够进一步改进网格的划分，能够进行贴体划分，再结合GKS理论进行计算分析。

基金项目

本文来自上海市科学技术委员会科研计划项目，基金号：19511104100。

NOTES

^*通讯作者。