0 引 言

当飞行器靠近地面、水面或冰面等无孔介质时^[1]，受到地面影响靠近地面的流体将沿着平行于地面的方向流动^[2]，这就是地面效应。飞机起飞着陆过程中离地高度很近，不可避免地会受到地面效应的影响。飞机起降时要尽可能降低飞行速度，缩短滑跑距离，常采用多段翼型来满足起飞时高升阻比或着陆时高升力的要求。研究表明多段翼型的地面效应与单段翼型有所不同^[3]：随着离地高度的减小，单段翼型的升力增加，而多段翼型的升力可能减小。因此，为保证飞机起降气动性能和操纵性能，有必要开展多段翼型地面效应研究。

目前，关于翼型地面效应的研究大都是固定飞行高度的定常状态下单段或多段翼型地面效应的研究。Ahmed等^[4]和Justin等^[2]分别对NACA4412翼型在不同迎角下的地面效应和低雷诺数下的地面效应进行了实验研究。Zhang等^[5]对反转的两段翼型的地面效应进行了实验研究。肖涛^[6]等对面效应下加装Gurney襟翼的NACA0012机翼进行风洞实验，测量其气动力，并采用PIV技术观察襟翼在地面效应下速度场的分布。

秦绪国等^[7]使用商用软件对固定离地高度下多段机翼的定常地面效应进行了数值模拟。近年来，已有对单段翼型非定常地面效应的研究。Ye-Hoon Im等^[8]对NACA6409靠近非定常波形物面时的流动进行了计算分析。Wu等^[9]应用浸入边界格子玻尔兹曼方法对层流中离地高度固定情况下扑翼前飞时受地面的影响进行了数值研究。对多段翼型非定常地面效应的研究则很少，王忠清等^[10]使用重叠网格方法研究了翼型在某离地高度下迎角变化对非定常流动的影响。但鲜有关于翼型下落靠近地面的非定常地面效应研究。本文模拟了翼型的着陆过程，研究其离地越来越近时非定常地面效应对翼型周围流动的影响。

针对着陆过程存在的翼型大位移运动问题，直接应用常规动网格方法会导致翼面附近网格密度分布差异较大而产生较大的计算误差，提出尽量保持翼型周围和尾迹区域网格跟随翼型做刚体运动的动网格策略，并采用一种简单的网格重构方法以保证翼型下降过程中具有合理的网格分布。据此计算了NACA0012翼型和GAW-(1)两段翼型匀速下降问题，非定常升力计算结果均和相应准定常状态相吻合。进一步建立起多段翼型非定常地面效应N-S方程计算方法，计算GAW-(1)匀速下降过程的非定常地面效应，并与下降速度折合迎角加至翼型的准定常计算结果进行比较，揭示两者之间的差异。

直角坐标系下基于任意拉格朗日欧拉描述法的二维非定常N-S方程为^[11]：

点击浏览公式

采用有限体积空间离散法、双时间推进法求解非定常N-S方程。远场采用无反射边界条件。

为了避免网格变形引入的误差，Thomas和Lombard^[12]提出了几何守恒律(GCL)。由质量守恒方程可以推出其积分形式的表达式为：

点击浏览公式

基于RANS模拟湍流，湍流模型采用一方程Spalart-Allmaras模型，该模型比较适合多段翼型的粘性绕流数值模拟^[13]。SA模型的控制方程的积分形式为^[11]：

点击浏览公式

快速生成高质量的动网格是动边界非定常问题数值计算的基础。分块结构网格是复杂组合体和非定常流动数值模拟的理想网格形式，其动网格通常在保持分块网格拓扑结构和对接关系不变的前提下，重新分布网格点，以避免搭接网格的再分区过程引入插值误差。

利用基于径向基函数(RBF)插值和超限插值(TFI)的动网格方法^[14]。该方法分为两步：第一步，选择变形已知的各网格块顶点作为中心点，并通过径向基函数插值算出相应网格块边界的变形。径向基函数插值时使用插值函数f定义物理空间内的位移，表示为径向基函数的线性组合：

点击浏览公式

插值系数a_i、b₀、b₁、b₂可通过附加约束求得：

点击浏览公式

第二步，在各网格块内根据边界的变形进一步通过超限插值得到其余网格点的坐标。二维超限插值公式如下：

点击浏览公式

为保证初始网格的分布规律，使用基于弦长的超限插值方法，i、j方向控制函数的定义如下：

点击浏览公式

点击浏览公式

多段翼型流场存在尾迹流动的相互干扰、流动分离等复杂现象^[15]，数值模拟存在较大困难。以GA(W)-1两段翼型(图 1)为例，验证本文N-S方程计算方法数值模拟多段翼型定常流动的可靠性。翼型缝道参数见表 1，其中c为计算参考弦长，取为1 m。计算马赫数Ma=0.2，雷诺数Re=2.2×10⁶。

图 1 GA(W)-1两段翼型示意图 Fig. 1 Schematic of GA(W)-1 airfoil

图选项

图 2是GA(W)-1两段翼型的本文定常升力系数计算结果与实验值^[16]的对比，两者吻合度很高，说明采用本文N-S方程计算方法数值模拟多段翼型定常流动是可靠的。

图 2 GA(W)-1定常升力系数 Fig. 2 Steady lift coefficient of GA(W)-1

图选项

然后以非定常运动NACA0012翼型为例，验证本文计算方法模拟非定常问题的可靠性。计算马赫数Ma=0.755，雷诺数Re=5.5×10⁶，参考弦长取为1 m。翼型的运动规律表示为迎角随着时间变化的函数α=0.016+2.51sin(ωt)，其中ω是角频率。本文计算结果与参考文献计算结果及实验值^[17]吻合得很好，并更接近实验值(图 3)，说明采用本文计算方法数值模拟翼型小振幅非定常运动问题是可靠的。

图 3 非定常运动NACA0012翼型法向力系数 Fig. 3 Comparison of normal force

图选项

翼型匀速下降问题理论上等效为一个翼型固定不动且将下降速度考虑成折合迎角的准定常问题。第2节所采用的非定常运动NACA0012翼型是验证2D非定常算法的一个典型算例，很多文献给出了计算值和实验值比较，计算与实验之间均存在如图 3所示相当的误差，有计算和实验两方面原因。因此，采用非定常算法计算翼型匀速下降问题并与等效准定常问题进行比较，是发展和完善非定常计算方法的有效手段，也是本文研究多段翼型下降过程中非定常地面效应的基础。

NACA0012翼型匀速下落。计算马赫数Ma=0.2，雷诺数Re=2.88×10⁶，来流迎角为0°，参考弦长取为1 m。翼型以w₀ =3.5637 m/s匀速下落，折合迎角为3°。为尽量减小下降过程中远场离翼型太近引起的误差^[18]，初始静态网格的远场取为50倍参考弦长。如图 4所示，采用C型网格，网格数1.1万左右。

图 4 NACA0012翼型初始网格 Fig. 4 Schematic of computational grid

图选项

考虑到翼型不断下降会导致翼型大位移移动，本文采用两种动网格方法，研究动网格对计算结果的影响。图 5给出了t=1.5 m时两种动网格。(1)代表使用基于径向基函数(RBF)插值和超限插值(TFI)的动网格方法；(2)代表翼型周围及尾迹区网格随翼型一起做刚性运动，在此基础上应用扰动衰减规律生成其余部分动网格^[19]。

图 5 t=1.5 s时刻NACA0012翼型两种动网格 Fig. 5 Dynamic grid of t=1.5 s

图选项

经物理时间步长/虚拟定常内迭代步数组合与收敛性验证后，非定常计算时间步长均取为1.0×10^-3 s，内迭步数代取为50步。图 6为两种动网格下算出的升力系数随时间的变化过程，表 2给出t=1.5 s两种动网格下的升力系数及准定常计算结果。比较分析表明：基于动网格方法(1)计算出的升力相比于准定常值误差很大，达到了-2.267%；而基于动网格方法(2)计算出的升力与准定常值相吻合，仅相差0.504%。对于翼型下降问题，直接应用动网格方法导致一部分网格被压缩，而一部分网格被拉伸，随高度降低出现网格大变形情况，此时翼型上下表面网格高度出现较大甚至很大差异，引起计算误差。因此，动网格方法的直接应用不能很好地保证物面附近及尾迹区的网格质量，为提高非定常数值计算精度，建议采用保持该区域网格与物体一起做刚体运动的方法生成高质量动网格。

图 6 匀速下落NACA0012翼型升力系数历程 Fig. 6 Lift coefficient evolution of NACA0012

图选项

计算分析GA(W)-1两段翼型匀速下落的非定常地面效应问题。计算马赫数Ma=0.2，雷诺数Re=2.2×10⁶，参考弦长取为1 m，来流迎角取为0°。初始远场为50倍弦长、翼型迎角为4.8°，匀速下落速度为w₀ =3.5637 m/s，折合迎角为3°。经物理时间步长/虚拟定常内迭代步数组合与收敛性验证后，非定常计算时间步长取为1.0×10^-3 s，内迭步数代取为400步。

根据第3节结论，非定常计算应尽量保持物面附近及尾迹区网格跟物体做刚体运动，为便于翼型下落过程中动网格生成其余区域采用H型网格。计算时使用C-H-O型分块结构网格拓扑，共10个区(图 7)，图 8是相应的分块网格图，网格总数约9.3万。非定常计算时，翼型周围的网格区5~9(见图 7)随翼型做刚体运动；1~4区仍然使用基于扰动衰减的动弹性网格方法^[19]进行变形；网格区10为最下方的地面，计算过程中保持不变。

图 7 GA(W)-1两段翼型分块网格拓扑结构图 Fig. 7 Schematic of grid topology

图选项

图 8 GA(W)-1两段翼型分块网格图 Fig. 8 Schematic of computational grid

图选项

随着翼型向下运动，1区的网格越来越稀，2区的网格越来越密。因此数值模拟多段翼型的着陆过程必然涉及到网格的大变形。为保证合理的网格分布，本文在上述动网格基础上采用一种简单易行的网格重构方法。二维情况下，以垂直I方向为例，将初始网格的第i层平均高度记为h_i，动网格变化后第i层平均高度记为h′_i。当第i层网格被拉伸至

点击浏览公式

类似的，当第i层网格被压缩至

点击浏览公式

图 9给出了GA(W)-1两段翼型运动到地面附近时的动网格。

图 9 GA(W)-1两段翼型运动至地面附近时的动网格 Fig. 9 Schematic of grid as GA(W)-1 approaches ground

图选项

2 s时翼型离地高度h(定义为襟翼偏转后翼型后缘点到地面的距离)约为42.9c，此时地面对翼型作用可以忽略不计，而翼型匀速下落状态已趋于稳定。而后对2 s时刻的动网格进行定常计算，得出准定常气动力结果。非定常计算稳定后的升力系数为C_{L unsteady}=3.411，与相应准定常状态的升力系数C_{L steady}=3.396相比，大了约0.44%。图 10是非定常、准定常算法计算得到的翼型压强系数分布，两者相吻合。

图 10 GA(W)-1两段翼型压强分布 Fig. 10 Pressure distribution of GA(W)-1

图选项

为研究非定常地面效应，以下计算翼型继续以w₀ =3.5637 m/s匀速下落至物面附近过程中的气动力。同时使用翼型运动后的网格，将非定常匀速下落等效为来流迎角为3°加至翼型进行定常计算。将两者结果进行对比。计算时物面边界条件根据实际情况进行选取，翼型表面满足：

点击浏览公式

图 11是不同离地高度下准定常情况(下降速度折合迎角3°加至翼型)和非定常计算得到的升力系数比较图。

图 11 GA(W)-1两段翼型升力系数对比 Fig. 11 Lift coefficient comparison of GA(W)-1

图选项

(1) 两种情况下的升力系数均随离地高度的减小而减小，文献^[7]给出了地面效应导致多段翼型升力减小的原因；

(2) 随着翼型的下落，非定常地面效应对气动力的影响发生改变(表 3)。在h/c<h₀的时候，非定常地面效应小于准定常地面效应，非定常升力的减小小于准定常升力的减小；在h/c>h₀的时候，翼型的下落运动和非定常地面效应的综合作用越来越明显，导致非定常升力略小于准定常升力(h₀≈3)。

本文建立起多段翼型非定常地面效应N-S方程数值计算方法，能应用于单段或多段翼型起飞、着陆过程中地面效应的实时数值模拟。相比于传统的不考虑翼型起降运动的准定常地面效应算法，本文首次得到更为真实的模拟翼型的地面效应，其差量是实际应用中应尽可能加以考虑的。本文结果可为增升装置设计以及飞机起降气动性能和操纵性能评估提供参考。