0 引言

倾转旋翼飞行器由于特殊的构型和工作条件，流场结构和气动力特征复杂。如悬停工作状态时，旋翼尾迹受到机翼翼面阻挡，在机翼下方诱导出大范围的分离，在机翼上表面中线附近汇集形成喷泉流动；前飞状态旋翼脱出尾迹与机翼相互干扰，此时旋翼作用类似于螺旋桨；前飞和悬停的过渡状态时，旋翼尾迹弯曲严重，非定常气动现象显著。

有关倾转旋翼的计算研究工作已经开展不少，计算方法包括自由尾迹分析^[1-4]和求解Navier-Stokes方程计算^[5-7]等。在采用CFD方法时，各种层次的模型也被广泛使用。由于旋翼流场结构复杂，因而相应计算资源要求也通常较高，选择合适的计算方法，使得计算资源在可承受的范围之内，又能够有效描述气动现象，是一项很重要的课题。求解雷诺平均Navier-Stokes(Reynolds Averaged Navier-Stokes，RANS)方程是被广泛采用的方法，也有一些研究者使用高计算量的大涡模拟(Large Eddy Simulation，LES)^[8]方法进行气动干扰分析，但在工作中，他们把旋翼对流场的作用进行了简化，即采用动量源模型模拟旋翼，使得整体计算的网格数量得到有效控制，这样虽然更高精度地模拟了旋翼下洗流通过机翼时的干扰效应，但旋翼周期运动及桨叶外形影响的敏感性带来的影响被削弱；由于计算量的考虑，一些研究者在旋转坐标系下求解旋翼附近流场区域^[9]，大大节约了计算资源，但由于忽略了非定常效应影响，且惯性系与非惯性系区域间数据交换通常会产生尾迹截断，导致流场数值模拟的精确性下降；一些研究者仅考虑一侧旋翼、短舱及半机身的模型，附加了对称面进行计算^[10]以削减计算量，这样对于双旋翼间相互作用的动态时间影响效果则被忽略。近年来，有关倾转旋翼的数值模拟在国内成为热点，研究者们开展的典型研究工作如过渡态计算方法研究^[11-12]和悬停状态的气动干扰分析等^[13]。

本文主要对倾转旋翼的悬停状态进行数值模拟研究。与先前研究工作不同的是，本文对采用的数值计算方法、计算模型和网格系统进行了综合考虑。首先计算模型真实模拟整个机翼和双旋翼的所有桨叶外形，由于不采用对称面、动量源条件并在地面坐标系下求解，所有气动部件的运动及相互影响被更真实刻画；其次计算网格方面采用了可自适应的直角背景网格，并预先剖分尾迹区域，基于流场变量的周期记录统计，对高物理量梯度区域(特别是旋翼尾迹区域、喷泉流动区及机翼下方的流动分离区域)实现自动网格细分，因而网格数量得到有效控制；最后在控制方程求解方面，采用了DES方法进行计算，对比LES方法对网格数量及计算量的较高要求，DES方法并没有比RANS方法增加过多额外的计算负担，且其计算近壁面分离流动的能力更好，是一种当前可行的研究手段。

针对典型工作状态，采用以上方案，对倾转旋翼双旋翼/机翼干扰模型流场和气动力进行了数值计算，特别是对桨叶/机翼气动干扰、机翼前后缘下方流动分离现象和机翼向下载荷变化等进行了定量分析。通过与相同状态RANS计算结果的比较，分析了采用DES模型时，在模拟旋翼、机翼气动载荷及机翼下方涡演化等流场方面的计算结果差异。

1 模型、计算状态和网格

由于很难查阅到国外倾转旋翼飞行器外形的精确数据，计算模型选用了双旋翼/机翼的干扰类比模型。其中每个旋翼由三片桨叶组成，桨叶半径R为5.8m，有较大和不规则的扭转分布，机翼选择有上反和小前掠的平面矩形机翼模型，半翼展长L为6.2m。

计算状态参考典型倾转旋翼工作情况，选择旋翼桨尖马赫数0.7。桨叶总距角7.5°，选择此总距角的原因在于本文主要讨论不同方法计算的旋翼/机翼间气动干扰差异，此时，旋翼表面气动分离情况不严重，不同方法计算的旋翼气动力差异不会对旋翼/机翼间作用讨论带来更大的干扰影响。

嵌套网格系统由八块网格组成，即分别围绕六片桨叶的贴体网格块，围绕机翼的贴体网格块及包围机翼和桨叶网格块的背景网格块。贴体网格块中由于物体外形并不复杂，均生成纯六面体网格。背景网格采用了自适应直角网格，结合问题实际，进行了网格初始设定和自适应调整：一是初始网格设定，先给定一单元尺度较大的均匀背景网格；二是根据桨叶、机翼贴体网格块上的物体表面空间位置(桨叶做一周期的旋转运动)和网格尺度，寻找背景网格上对应网格区域进行细分和自动网格尺度调整；三是预定尾迹区域网格调整，根据经验给定旋翼尾迹大致范围内区域，按照桨叶弦长预估涡核尺度，细分可能的尾迹区域到合适尺度；四是进行流场初步计算，根据流场物理量反馈，划分变量高梯度区域，生成最终计算网格。经过以上步骤得到的计算网格如图 1所示。

定义向上为Y轴方向。旋翼桨盘向后方向为Z轴方向，X轴方向用右手定则确定。图 2给出了本文用于计算讨论的桨叶方位角位置(俯视图)。

2 数值方法 2.1 控制方程形式

在地面坐标系下，不计体力等产生的源项，采用格心格式的有限体积法，积分形式的雷诺平均Navier-Stokes方程可以写为：

$ \frac{\partial }{{\partial t}}\int\limits_\Omega {\mathit{\boldsymbol{W}}d\Omega + \oint\limits_{\partial \Omega } {\left( {{\mathit{\boldsymbol{F}}_c}-{\mathit{\boldsymbol{F}}_v}} \right)} {\rm{d}}S} = 0 $

(1)

$ \begin{array}{l} {\mathit{\boldsymbol{F}}_c} = \left[{\begin{array}{*{20}{c}} {\rho V}\\ {\rho uV + {n_x}p}\\ {\rho vV + {n_y}p}\\ {\rho wV + {n_z}p}\\ {\rho HV} \end{array}} \right], {\mathit{\boldsymbol{F}}_v} = \left[{\begin{array}{*{20}{c}} 0\\ {{n_x}{\tau _{xx}} + {n_y}{\tau _{xy}} + {n_z}{\tau _{xz}}}\\ {{n_x}{\tau _{yx}} + {n_y}{\tau _{yy}} + {n_z}{\tau _{yz}}}\\ {{n_x}{\tau _{zx}} + {n_y}{\tau _{zy}} + {n_z}{\tau _{zz}}}\\ {{n_x}{\mathit{\Theta }_x} + {n_y}{\mathit{\Theta }_y} + {n_z}{\mathit{\Theta }_z}} \end{array}} \right]\\ {\mathit{\Theta }_x} = u{\tau _{xx}} + v{\tau _{xy}} + w{\tau _{xz}} + K\frac{{\partial T}}{{\partial x}}\\ {\mathit{\Theta }_y} = u{\tau _{yx}} + v{\tau _{yy}} + w{\tau _{yz}} + K\frac{{\partial T}}{{\partial y}}\\ {\mathit{\Theta }_z} = u{\tau _{zx}} + v{\tau _{zy}} + w{\tau _{zz}} + K\frac{{\partial T}}{{\partial z}} \end{array} $

其中，Ω为控制体的体积，S为积分面面积，W为守恒变量，F_c和F_v分别为对流和粘性通量项，ρ为密度，p为压强，V为垂直网格交接面的对流速度，u、v、w为速度，n_x、n_y、n_z为交接面单位外法矢，H为总焓，τ为应力项，T为温度，K为传热系数。

涡粘性采用一方程SA模型^[14]及其对应DES方法^[15]进行，积分形式的方程可以写为：

$ \frac{\partial }{{\partial t}}\int\limits_\mathit{\Omega } {\tilde \upsilon {\rm{d}}\mathit{\Omega } + \oint\limits_{\partial \mathit{\Omega }} {\left( {{F_{c, T}}-{F_{\upsilon, T}}} \right)} {\rm{d}}S} = \int\limits_\mathit{\Omega } {{Q_T}{\rm{d}}\mathit{\Omega }} $

(2)

其中${\tilde \upsilon } $为粘性系数，F_{c, T}为对流通量项，F_{υ, T}为粘性通量项，Q_T为源项，其中包含生成项和破坏项等。

DES计算时，不采用最近壁面距离，而是采用混合尺度，特征长度计算采用下式

$ \tilde d = \min \left( {d, {C_{{\rm{DES}}}}\Delta } \right) $

(3)

d为计算单元到所有壁面的最近距离。因计算采用格心格式，Δ为当前计算单元中心到其所有邻居单元中心的距离最大值，C_DES为可调节系数。由于采用了混合尺度，该方法综合了大涡模拟方法和RANS方程的优点，在近壁面区域，计算出的尺度与RANS方法相同，在远离壁面区域，湍流封闭模式是亚格子模型。

需要注意的是，在旋翼流场计算中，随桨叶运动，固定的背景网格上计算点到壁面的距离不断变化，因而需要在每个物理时间，重新计算特征长度。

2.2 求解方法

无粘通量计算时，空间方向采用重构方法^[16]构建二阶ROE^[17]格式计算交接面上值，粘性通量计算采用中心格式。主控和湍流方程计算采用隐式时间推进^[18]，为提高非定常计算的效率，还耦合使用了双时间方法^[19]。物理时间步长为桨叶旋转过方位角1/4°的时间，每个物理时间内迭代20次，计算9个周期后，得到相对稳定的计算解。

3 算例验证 3.1 NACA0015翼型较大迎角状态气动力计算

为验证本文DES方法的计算能力，进行了NACA0015翼型在较大迎角(17°)状态下的气动力计算。计算的翼型弦长是0.3048m，来流马赫数是0.29。图 3是计算得到的翼型表面压强与试验值^[20]的对比。可以发现对于较大迎角计算状态，翼型上表面附近发生大范围的分离流动时，相对RANS计算，采用DES方法可以更好地预测翼型上表面的压强分布。

3.2 旋翼悬停气动力计算

为验证本文方法在旋翼气动力计算方面的能力，选用了被研究者们广泛采用的有试验结果可供对比的“Caradonna & Tung旋翼”^[21]为研究模型。该旋翼桨叶的翼型为NACA0012，展长R是1.143m，展弦比为6，无扭转尖削。为体现出DES方法与RANS方法的计算差异，计算状态选择为桨尖马赫数0.794，桨距角为12°。在该状态下，桨叶靠近桨尖的部分截面上出现激波，激波和附面层相互干扰，形成桨叶上表面的分离流动。图 4给出了采用DES和RANS方法计算的桨尖附近的桨叶上表面压强分布。观察DES结果，可以发现随着桨叶旋转(方位角变化)，桨叶上表面激波后方形成的气流分离区域逐渐向桨叶尾缘推移(反映在桨尖附近的桨叶表面压强随方位角变化)，到300°方位角时几乎消失。而采用RANS方法时，各个物理时刻计算出的表面压强基本相同，表示其未能捕捉到激波后方时变的分离流动现象。

图 5为计算得到的桨尖附近截面(r=0.96R)表面压强系数随桨叶方位角变化结果。可以发现，在该截面，DES方法计算的桨叶翼型上表面压强随方位角变化明显。采用RANS计算的结果几乎随方位角不变，为清晰起见，只给出其0°方位角时的结果。

图 6为计算得到的各个桨叶截面表面压强系数与试验值的对比(从旋翼旋转第10个周期开始进行3周的气动力平均结果)。在大部分截面上，计算和试验结果都有较好的一致性。在非常靠近桨尖的截面(r=0.96R)，计算和试验结果的激波位置存在一定差异。比较RANS和DES计算结果，发现其差别体现在靠近桨尖侧的截面(r=0.89R，r=0.96R)，对于远离桨尖的截面，两者差异消失。究其原因在于，对应于该计算状态，桨叶表面的分离仅发生在靠近桨尖的激波后方，其对远离桨叶尖部的桨叶表面附近流动影响有限(r=0.68R，r=0.8R)。在靠近桨尖附近的区域，DES计算得到的激波位置略偏前缘，而激波后的桨叶上表面压强高于对应的RANS结果。

4 旋翼/机翼干扰状态计算结果分析 4.1 流动现象比较分析

在倾转旋翼的悬停状态，流场中典型的流动状态即喷泉流动、旋翼下洗及旋翼尾流和机翼干扰后的分离流动。其中标志性的喷泉流动是由于两个旋翼的尾流冲击机翼翼面后沿翼面向双旋翼中线处汇集并卷起后形成。本文计算出的结果如图 7所示(左为DES结果，右为RANS结果)。由图可以看出，采用DES方法计算出的流场无论是机翼上方中线附近的喷泉流动区域还是机翼下方的诱导流动，都体现出了左右不对称的时变特征，而采用RANS方法计算出的流场则对称性相对较好。

图 8中给出了桨叶方位角0°时，计算得到的旋翼尾流受机翼阻挡，在机翼下方不同翼展平面诱导出的分离流动情况(左为DES结果，右为RANS结果)。由图可见，旋翼诱导的尾流冲击机翼，受机翼阻挡，形成机翼下方的旋转分离流动。采用RANS计算的结果一般是预测出靠近机翼前缘和尾缘下方的两个较大范围的涡，而采用DES方法则计算出较多的向下发展的连续小涡。DES结果给出了更细致的机翼下方流动分离情况的描述。在靠近机翼尖部位置(X/L=0.8)，两种计算方法得到的流动形式相接近，这可能是两个原因导致的，一是该位置是旋翼桨根的下方，诱导的下洗速度相对较小，二是机翼尖部的三维横向流动效应削弱了诱导分离的影响。

图 9给出了采用DES方法计算得到的不同桨叶方位角时的特征空间截面(对应机翼位置X=0.4L)上的涡量，可以发现本文计算较为清晰地捕捉到了旋翼尾迹及机翼下方脱落涡的时间发展变化历程。

4.2 气动力比较分析

表 1给出了本文计算得到的旋翼(一个)和机翼整体气动力结果。由对比得知，采用DES方法计算得到的旋翼拉力和扭矩都大于RANS结果，但旋翼拉力的差别极小。最大的差别在于机翼向下载荷结果，采用DES方法计算出的机翼向下载荷较小，比RANS结果小2.3%左右。导致此项差别的原因可能是采用不同方法时，计算得到的机翼下方诱导流动(特别是机翼前后缘的涡脱落)形式不同。由计算还可以得知，采用DES计算的机翼下洗载荷占整体旋翼拉力的8.34%左右，小于RANS方法对应结果。

图 10给出了一片桨叶在一个旋翼旋转周期内沿截面段的拉力分布变化。可以发现在大部分桨叶方位角，两种计算方法差异较小，两者主要差异是桨叶240°到300°方位角之间(正好对应于桨叶通过机翼上方位置)，采用DES方法计算得到了更为剧烈的拉力变化。由图还可以得知，在桨叶通过机翼上方时，拉力最小，在远离机翼的方位角时，桨叶拉力得到恢复。

图 11给出了几个特征方位角时的桨叶拉力系数比较，可以发现在240°到300°方位角之间，桨叶拉力下降主要来源于桨尖部分的气动力损失。发生这种现象的原因可能是该时刻桨叶位于机翼上方，两侧旋翼桨叶靠近且均受到喷泉流动形成的向下气流影响，桨叶靠近桨尖截面的有效迎角减小。

图 12给出机翼向下载荷的周期变化，由于每个旋翼都由三片桨叶组成，因而旋翼下洗作用在机翼上，形成了周期性的机翼向下载荷。DES方法计算出的机翼向下载荷相对较小，且出现峰值的相位也略靠前。

图 13给出了机翼上表面的压强等值图(左为DES结果，右为RANS结果)，可以看出在不同桨叶方位角时的机翼上表面高压区域的变化，反映出旋翼对机翼的下洗作用影响。

5 结论

在嵌套网格系统下，基于DES方法，进行了倾转旋翼双旋翼/机翼干扰模型的流场和气动力计算，并与RANS结果进行了对比研究，结果表明：

1) 由于计算的桨叶总距角较小，DES方法得到的旋翼整体气动力与RANS方法相差不大，但在机翼向下载荷计算方面，DES方法计算得到的向下载荷较小；

2) 在桨叶通过机翼上方时，拉力出现周期内最小值，原因可能是向下的喷泉流动减小了桨尖部分截面的迎角。采用DES方法计算时，桨叶通过机翼上方过程的拉力变化更剧烈；

3) 采用DES方法计算的流场左右对称性相对RANS方法要差些，非定常效应现象更加显著。机翼下方靠近前缘和后缘的涡结构形式描述也更为细致，表明了该方法在壁面附近存在较大分离流动时的模拟效果更好。