2. 中国直升机设计研究所, 景德镇 333000
2. China Helicopter Research and Development Institute, Jingdezhen 333000, China
“沙盲(Brownout)”是指当直升机在干燥，含有大量沙尘颗粒、雪花等离散相的近地面飞行时，沙尘被卷起形成沙尘云，严重阻碍飞行员视线、进而诱发严重事故的特殊现象[1-4](图 1)。
|
| 图 1 直升机悬停时沙尘云中沙尘颗粒的宏观分布 Fig. 1 Macroscopic distribution of dust particlesin dust cloud during helicopter hover |
当直升机在地效悬停或近地低速机动时，其桨尖涡脱落并经过一段距离后抵达地面，然后与地面作用产生反向的上洗并将地面沙尘颗粒从旋翼附近卷起[5-6]。这些被卷起的沙尘颗粒到达旋翼上方后又被卷入旋翼下洗流中，这个过程循环往复，最终形成遮天蔽日的沙尘云，严重侵害飞行安全，且这些沙尘颗粒对桨叶和机身还会造成磨损并减少使用寿命，同时还会导致机载设备失效，甚至对周围起降的飞机和人员安全造成严重威胁。因此，深入研究沙盲产生的机理、形成过程并分析沙尘颗粒的运动状态和分布规律具有重要意义。
早期研究沙盲的方法主要是基于PIV(Particle Image Velocimetry)技术对桨尖涡与地面的相互作用进行观察。为了测量地面涡以及周围的流场速度信息，Ganesh和Komerath使用该技术对前进比0.05~0.08范围内的缩比模型所产生的地面涡进行了观察，并对沙盲进行了可视化研究[7]；Lee等基于PIV对不同悬停高度下的缩比模型进行了测量并发现在中间悬停高度下，涡的伸展要比涡的耗散现象更为明显，但是随着悬停高度的增加，这个趋势则逆转了[8]。Milluzzo使用该技术对固定在1R高度上的4种不同桨尖外形的旋翼缩比模型进行了研究，分别测量了涡核尺寸以及800°方位角时刻的涡剖面速度[9]；与先前的研究不同，Hance则引入了机身，以研究沙盲环境下，机身对旋翼桨尖涡的影响[10]。Continuum Dynamics, Inc.(CDI)的Whitehouse等采用缩比旋翼，其桨尖转速达到200 ft/s时(弦长设为1 in)能匹配全尺寸的旋翼所产生的涡核强度。在对旋翼进行沙尘试验发现，通过在桨叶后缘处布置小翼可以加速桨尖涡的耗散[11]，阻碍桨尖涡在地面处的叠加[12]，从而减弱沙尘云的浓度，并认为改变桨叶气动外形会直接影响沙盲严重程度。随着计算机性能的快速发展，数值模拟方法在对地效流场和参数设定方面的优势逐渐体现并被应用于沙盲研究之中，Syal和Leishman采用基于拉格朗日框架下的时间精确自由涡尾迹方法，并结合应用于每个沙尘颗粒上的力学方程对沙尘云分布进行了求解，得到了可视化的沙盲结果[3]。Thomas等通过求解雷诺平均Navier-Stokes(RANS)方程，模拟了旋翼的地效流场，并将结果推广到沙尘颗粒的影响研究中[13]，这些对桨尖涡与地面的相互作用以及沙尘颗粒的力学分析进行的研究为深入了解沙盲机理奠定了基础。但全尺寸的直升机沙盲试验由于难以参数化，代价太高且存在相当的危险性，而在室内或者风洞中的缩比模型进行的沙盲研究虽然能取得和全尺寸直升机接近的涡核强度以及流场细节，但仍然在不少方面与真实环境存在差异，例如：真实环境中的沙尘在颗粒粒径、粒径范围和颗粒数量等参数上与实验室中不尽相同，因此，采用缩比模型的试验仍难以准确反映出全尺寸直升机沙盲现象的状态，再者，之前的数值方法没有考虑沙尘颗粒之间以及沙尘颗粒(由沙尘颗粒构成的地面)与桨叶、机身发生的碰撞，而它们也会对沙盲结果产生重要影响。因此，针对真实尺寸的直升机地效流场进行分析并将沙尘颗粒物理属性以及碰撞影响等因素考虑进耦合计算当中十分必要。
鉴于此，本研究拟构建一种基于计算流体力学-耦合离散单元法(CFD-DEM)用来模拟直升机沙盲现象。这是一种基于Euler-Lagrange模型求解流体、离散相以及耦合方程获得颗粒动量、质量和能量信息的方法，该方法通过迭代计算确定出每一个时间步长颗粒的受力和位移，从而追踪单个颗粒的微观运动并得到离散相的宏观运动规律。Yu等将这种耦合方法应用于煤矿开采行业中，模拟掘进中产生的粉尘运动规律，并得到了不同粒径的沙尘颗粒在不同时刻的分布特征[14]；Hou等模拟了在火星大气环境中的非球形沙尘颗粒在探测器上的分布，得到了不同流场速度下沙尘颗粒的速度分布规律[15]，这些研究均取得了和试验接近的结果。说明采用CFD-DEM对在流场中受到流场影响的微小粒径的沙尘颗粒运动轨迹和分布规律等进行分析是有效的。
基于以上研究成果，本文应用CFD-DEM方法研究沙盲时，采用了3种不同外形的非球形颗粒加入耦合计算当中，用来模拟更加真实的沙尘颗粒外形并同时赋予其真实物理属性，同时添加了两类碰撞接触模型。首先通过与已有的试验数据对比验证了耦合计算的准确性，之后，针对直升机不同悬停高度，从沙尘云中沙尘颗粒速度和分布规律两方面，深入分析了流场中沙尘颗粒运动的机理，得出了一些有实际意义的结论。
1 CFD-DEM数值计算方法 1.1 CFD求解器和湍流模型采用基于求解RANS方程的求解器STAR CCM+对悬停流场进行求解，对于差分表面积为da的任意控制体积V，密度基的耦合流体求解器采用笛卡尔积分形式的控制方程组可以写为[16]
| $ \mathit{\pmb{{\Gamma}}} \frac{\partial}{\partial t} \int\limits_\boldsymbol{V} \boldsymbol{Q} \mathrm{d} \boldsymbol{V}+\oint(\boldsymbol{F}-\boldsymbol{G}) \cdot \mathrm{d} \boldsymbol{a}=\int\limits_{\boldsymbol{V}} \boldsymbol{H} \mathrm{d} \boldsymbol{V} $ | （1） |
式中：Q为Navier-Stokes方程中的原始变量；F和G分别为无黏通量项和黏性通量项；H代表体积力；Γ表示预处理矩阵。将方程(1)应用于网格单元0的网格单元居中控制体积时，得到以下离散方程组：
| $ \boldsymbol{V}_{0} \mathit{\pmb{{\Gamma}}}_{0} \frac{\partial \boldsymbol{Q}_{0}}{\partial t}+\sum\limits_{f}\left(\boldsymbol{f}_{f}-\boldsymbol{g}_{f}\right) \cdot \boldsymbol{a}=\boldsymbol{h} \boldsymbol{V}_{0} $ | （2） |
式中：ff和gf分别为通过面f的无黏和黏性通量；V0为网格单元0的体积；h为单元0的体积力；Γ0为在网格单元0中计算的预处理矩阵。连续性、质量和动量守恒方程通过耦合流体求解器采用隐式积分格式求解，流场中气体设置为可压缩的理想气体，采用两方程的Menter剪切应力输运(SST)k-ω湍流模型求解壁面黏性流动，采用隐式双时间方法[17]以进行时间推进，在每个伪时间步上CFL数设为30以达到计算效率和鲁棒性的平衡，物理时间步长则设定为每旋转1°所需的时间[18-20]。
1.2 网格系统和物理模型本研究采用重叠网格方法将整体计算域划分为3块区域，包含机身的背景区域含有约700万网格，两块相同的桨叶区域分别含有约150万网格，如图 2所示，所有区域均采用正六面体为核心的切割体网格以保证计算精度，地面边界条件为无滑移壁面，四周边界距离旋翼中心距离为20R(R为旋翼半径)，上边界距离旋翼中心为15R，边界条件均为压力出口。
|
| 图 2 重叠网格系统 Fig. 2 Overset grid system |
全尺寸直升机采用一副无扭转的矩形C-T(Caradonna Tung)旋翼作为旋翼系统[21]，桨尖马赫数固定为0.7，根切为20%，展弦比为6，参考弦长c=1.5 m。机身采用删除了垂尾和平尾的简化模型以便保留其对流场主要影响的同时降低网格数量。
为了分析沙尘云中沙尘颗粒的运动状态，对计算域按图 3所示的方式共划分出8块长方体区域，其中，区域1、区域5、区域3和区域7这4个矩形区域位于桨盘平面下方，旋转轴周向间隔均为90°，长宽高分别为6.67R、1.5c和1H(H为悬停高度)，区域1位于机头正前方；区域2、区域6、区域4和区域8则分别位于区域1、区域5、区域3和区域7的正上方，高度为3.5R，宽度和长度不变，此8个区域在轴向和径向范围涵盖了整个沙尘云范围。为了进一步分析沙尘云中沙尘颗粒在桨盘平面上、下层的颗粒粒径和数量分布以及在Z轴方向上的速度分布，将计算域以桨盘平面为中心划分为同样大小，长宽高分别为9R、9R和2.7R的上层和下层两块长方体区域，如图 4所示。
|
| 图 3 区域划分 Fig. 3 Regional division |
|
| 图 4 上、下层区域划分 Fig. 4 Division of upper and lower regions |
为了模拟从形态学上更加接近自然界中受到岩石磨损而形成的沙尘颗粒的真实外形，在EDEM软件中通过“多球法”[22]以球形为基本构型创建了3种不同外形的沙尘颗粒：标准的球形、有部分尖角的非球形和尖锐的长条形颗粒。这些非球形颗粒通过多个标准球形按照一定的空间排布叠加而成，其外形和轮廓投影列在表 1中，沙尘颗粒、机身和桨叶的物理参数列在表 2中，其中，恢复系数是指碰撞过程中动能恢复的度量，可以用碰撞前法向相对速度与碰撞后法向相对速度的比值表示，静摩擦系数是最大静摩擦力与接触面间面积的比值，滚动摩擦系数是滚动摩擦力矩与支持力的比值[14]。在本研究中，所有的沙尘颗粒按照构成它们的标准球形直径(5、10、50、100 μm)的不同而分为4组，这3种不同外形的沙尘颗粒数量在每一组分中都是均分的。同时，为了模拟自然界中较为真实沉积状态的沙尘厚度，总计数量为20 000的沙尘颗粒全部从位于旋翼正下方地面上的高为10 cm、长和宽均为6.67R的长方体区域在最初时刻同时静态生成，位置随机分布。
| 材料 | 颗粒密度/(kg·m-3) | 泊松比 | 剪切模量/Pa | 恢复系数 | 静摩擦系数 | 滚动摩擦系数 |
| 沙尘颗粒 | 1 400 | 0.5 | 1×108 | 0.5 | 0.4 | 0.05 |
| 桨叶，机身 | 2 700 | 0.3 | 7×1010 |
在本研究中存在着两种不同的接触类型：①沙尘颗粒与沙尘颗粒之间的碰撞接触；②沙尘颗粒与地面以及旋翼和机身的碰撞接触。为了更加准确地描述这些接触，沙尘颗粒与颗粒之间的接触采用对于离散相计算高效且准确的Hertz-Mindlin (No Slip)模型求解，接触力合力由以下方程给出：
| $ F_{\text {contact }}=F_{\mathrm{n}}+F_{\mathrm{t}} $ | （3） |
式中：Fn和Ft分别为法向分量和切向分量，即
| $ F_{\mathrm{n}}=-K_{\mathrm{n}} d_{\mathrm{n}}-N_{\mathrm{n}} v_{\mathrm{n}} $ | （4） |
| $ F_{\mathrm{t}}=\left\{\begin{array}{ll} -K_{\mathrm{t}} d_{\mathrm{t}}-N_{\mathrm{t}} v_{\mathrm{t}} & \text { if }\left|K_{\mathrm{t}} d_{\mathrm{t}}\right|<\left|K_{\mathrm{n}} d_{\mathrm{n}}\right| C_{f{\mathrm{s}}} \\ \frac{\left|K_{\mathrm{t}} d_{\mathrm{t}}\right| C_{f_{\mathrm{s}}} d_{\mathrm{t}}}{\left|d_{\mathrm{t}}\right|} \end{array}\right. $ | （5） |
式中：Kn和Kt分别为法向和切向刚度；Nn和Nt表示法向和切向阻尼分量；dn和dt分别为颗粒在法向和切向重叠量；vn和vt是颗粒在法向和切向上的速度分量；Cfs为静摩擦系数。
而沙尘颗粒与物面的接触模型可以看作捕获模式[14]，采用Hertz-Mindlin与JKR(Johnson-Kendall-Roberts)凝聚力接触模型，即
| $ F_{\mathrm{JKR}}=-4 \sqrt{\pi \gamma E^{*}} \alpha^{3 / 2}+\frac{4 E^{*}}{3 R^{*}} \alpha^{3} $ | （6） |
| $ \delta=\frac{\alpha^{2}}{R^{*}}-\sqrt{4 \pi \gamma \alpha / E^{*}} $ | （7） |
式中：γ为表面能；E*为当量杨氏模量；α为接触半径；R*为当量半径；δ为重叠量。
沙尘颗粒在流场中的运动受到诸多力的影响，在本研究中只考虑了主要的重力、曳力和浮力的影响，其中曳力公式采用Free Stream曳力模型，方程为[23]
| $ F_{\mathrm{drag}}=0.5 C_{\mathrm{d}} \rho_{\mathrm{g}} A\left|v_{\mathrm{p}-\mathrm{g}}^{\mathrm{rel}}\right| v_{\mathrm{p}-\mathrm{g}}^{\mathrm{rel}} $ | （8） |
| $ C_\text{d}=\\ \left\{\begin{array}{ll} 24 / R e & R e \leqslant 0.5 \\ 24\left(1.0+0.15 R e^{0.687}\right) / R e & 0.5<R e \leqslant 1000 \\ 0.44 & R e>1000 \end{array}\right. $ | （9） |
| $ R e=\rho v L / \mu $ | （10） |
式中：Cd为曳力系数；ρg为空气密度；A为颗粒投影到流场方向的面积；vp-grel为流场和颗粒间的相对速度；Re为雷诺数；ρ为流体密度；μ为动力黏性系数；v为流场的特征速度；L为特征长度。
由于沙尘颗粒在整个流场计算域的体积分数(Volume Fraction)小于10-6，这意味着流场能够对沙尘颗粒的运动产生影响，而反过来，沙尘颗粒几乎不对流场产生影响，因而可以通过API (Application Programming Interface)采用单向的Lagrange-Euler耦合模型进行耦合计算，该API采用C++语言编写，负责在流场与沙尘颗粒之间传输动量、质量和能量信息(本研究中只考虑传递动量信息)，随着信息的交换，沙尘颗粒在每一个瑞利时间步(Rayleigh Time Step)上获得动量并产生位移，最终由全部的沙尘颗粒在空间不同位置上的分布而形成宏观的沙尘云全貌。流场与沙尘颗粒之间通过API传递信息的过程如图 5[22]所示，固相和流体之间的耦合可以分为动量、能量和质量耦合3类，能量耦合主要用来衡量流体与固相间的热传递和能量传递速率，质量耦合主要用来衡量颗粒表面物质的吸附与解离，或者颗粒上、颗粒内的化学反应情况。动量耦合是用来衡量颗粒受到的曳力和流体的对流动量通量的情况(通常可以用动量耦合系数来表示)，而本文主要研究流场中颗粒受到流场曳力作用后的运动和分布规律，不考虑能量和质量变化，因而只需考虑动量耦合，能量和质量耦合则不在本研究考虑范围之内。
对于CFD-DEM耦合计算沙盲方法的验证可分为两部分，第1部分是基于Light[24]在有地效(IGE)和无地效(OGE)状态的多个悬停高度下的拉力系数(CT)试验数据，对本研究中用于IGE流场计算所采用的数值方法准确性进行验证(耦合算例中的网格结构和边界条件等均与验证算例中保持一致)。如图 6所示，无地效状态下的CT计算值和试验值对比接近，只有-2%~-4%的偏差，符合预期，但是有地面效应状态下(不同地效悬停高度状态列在图中, h为桨盘平面距离地面高度)的CT与试验值偏差增大至约10%，这是因为数值模拟中忽略桨毂会导致桨根处产生喷泉效应导致的，因而这些偏差是能够接受的[25]。为进一步验证数值方法对桨尖涡位置的捕捉能力，当悬停高度为0.52R时将桨尖脱落涡位置与试验值进行对比，如图 7所示，桨尖脱落涡在不同方位角处轴向和径向的位置和试验值吻合较好，因而可以认为该数值方法能够为CFD-DEM耦合计算提供准确的流场信息。
|
| 图 6 拉力系数的试验值和计算值对比 Fig. 6 Comparison of experimental and numericalthrust coefficients |
|
| 图 7 桨尖涡位置的试验值和计算值对比 Fig. 7 Comparison of experimental and numericalblade tip vortex location |
第2部分是基于Wong和Tanner[4]采用摄影测量法观察到的沙尘云外观与Leishman采用自由尾迹和沙尘动力学所计算得到的沙尘颗粒分布来对本研究所采用CFD-DEM耦合方法得到的结果进行的对比验证。在试验中，他们使用6台摄影机拍下了当EH-60L滑行通过着陆区产生沙尘云的整体外观图，并对其外观轮廓进行了分析。在文献[3]中则采用基于Lagrangian数值方法求解流场(自由尾迹方法)耦合基于Lagrangian数值方法求解沙尘颗粒的算法，得到了沙尘云中颗粒分布。本研究中的验证算法中采用了与EH-60L尺寸近似的简化模型，将CFD-DEM耦合计算得到的沙尘云外观与试验和该文献的沙尘云轮廓进行了对比，结果如图 8所示。需要说明的是在图 8中文献的数据是考虑了载重约7 400 kg，配平后的沙尘云数据，另一组数据由于被遮挡无法完全辨识而未列出。但是，考虑到试验中的真实环境和CFD-DEM耦合算法中的简化环境无法完全一致。例如，试验中EH-60L在飞行轨迹上存在着横向和纵向的加速度，且飞行高度也随时间发生变化；试验中土壤中的沙尘颗粒粒径连续分布，跨度极大，且被卷入到空中较轻的沙尘颗粒占绝大多数[2]；摄像头拍摄的照片由于像素所限，捕捉到的沙尘云较实际范围偏小。同时，与本研究的算法不同，文献[3]所采用的耦合算法并未考虑沙尘颗粒的体积和密度等属性，而是对位于地面非黏性流场中较小的沙尘颗粒进行半经验的方式处理，这些沙尘颗粒受到大于重力和颗粒内部力之和的非稳态力的作用从而提升，导致沙尘云外观与本研究有所差异。排除以上因素，试验结果中沙尘云轮廓和CFD-DEM耦合计算得到的沙尘云外观轮廓接近，并且从图 8(b)可以看出，沿飞行方向受到上洗作用影响而提升的沙尘颗粒初始位置接近，因而可以认为该耦合算法是有效的。
|
| 图 8 沙尘云的试验测量和耦合数值模拟对比(t=10 s) Fig. 8 Comparison of dust cloud for experimentalmeasurements and coupling numericalsimulation (t=10 s) |
图 9中给出了本研究计算的地效流场中涡量图以及X、Y、Z方向的速度云图。与直升机OGE悬停相比，直升机IGE悬停时，由于旋翼下洗流在地面附近受到阻挡无法向下运动，而趋于向外扩散，流场及气动特性变得复杂。从图 9(a)中可以看到，最初脱落的桨尖涡外形规则且光顺，但随着其越来越接近地面而变得粗糙且不稳定，导致这个现象的主要原因有两点：一是桨尖涡在脱落过程中撞击到机身上而破碎，从而不能维持初始外形；二是接近地面的桨尖涡受到地面流场的挤压而与后一个桨尖涡发生融合，导致外形扭曲。分析图 9(a)和图 9(d)中位于桨根部位出现的“喷泉效应”[11]，可以看出桨根涡与桨尖涡沿Z方向速度相反，这是由于数值计算中未考虑桨毂对流场的影响，从而导致下洗流在桨根处受到机身的挤压而向上运动形成的。分析对比图 9(b)和图 9(c)可知，在有地效流场中，位于地面附近的流场速度主要沿径向向外，而切向速度较小，且由于流场受到机身，尤其是尾部的影响，导致尾部区域速度大小不稳定且方向各异，从图 9(d)中可以看出，位于桨根处Z方向速度向上，桨尖处Z方向速度向下，而且在桨盘上方能看出两处明显的沿Z方向向下的速度，这是由于旋翼上方流场受到旋翼影响而被吸入导致的。
|
| 图 9 有地效流场中涡量图和X、Y、Z方向速度云图 Fig. 9 Vorticity structure colored by Z-velocity andX, Y, Z directional velocity contours withIGE flow field |
沙尘云中沙尘颗粒的分布按高度可分为位于桨盘平面上方(上层)的沙尘颗粒和桨盘平面下方(下层)这两块区域，而主要导致阻碍飞行员视线的是位于上层的沙尘颗粒，为了分析上、下层沙尘颗粒的运动状态，划分出如图 3所示的8个区域。当H=0.67R时，从图 10(a)和图 10(b)可以看出1、3、5和7区域的沙尘颗粒在耦合计算时间t≈5 s前，它们的径向速度大于切向速度。这表明位于地面附近，最初受到桨尖涡干扰的沙尘颗粒主要受到流场的影响而沿径向向外运动，同时也存在沿旋翼旋转的切线的分速度，这点可以从图 9(b)和图 9(c)中得到验证。同样的结果可以从H=1R的图 11(a)和图 11(b)与H=1.5R的图 12(a)和图 12(b)中看出，但随着悬停高度的增加，切向速度逐渐减弱，径向速度变化不大。这是由于悬停高度越低，旋翼和机身离地面越近，地面附近流场变得更为复杂导致的，因而沙尘颗粒初始运动状态也更加不稳定。但是，随着悬停高度的增加，由于桨尖涡需要经过更长距离脱落到地面而易于耗散，因而切向运动变得更加不明显。当耦合计算大约在5~15 s之间，沙尘颗粒的径向和切向速度均逐渐开始减小，这是由于在这段时间内，沙尘云范围不断扩大，远处流场对沙尘颗粒的影响减弱从而导致沙尘颗粒运动速度减小而造成的。但是随着沙尘云发展至15~20 s之间时，沙尘颗粒这两个方向的速度又变得更加不稳定，这是因为随着时间的推移，更多位于上层的速度不稳定的沙尘颗粒被卷入进旋翼下洗流中，导致下层颗粒的平均速度也变得不稳定。从图 11(a)和图 11(b)和图 12(a)和图 12(b)的数据来看，随着悬停高度的增加，上述现象同样存在，但是切向和径向的平均速度较H=0.67R时更加稳定。从图 10(c)和图 10(d)中可以看出，对比下层区域，当H=0.67R时，1 s时刻开始有沙尘颗粒进入上层区域，2 s时刻上层沙尘颗粒速度达到最大，但速度方向不一致，在切向和径向的平均速度并未像下层区域差别明显，而是更加相近，随着时间的推移，这个差异逐渐减小，且在各个方向上的平均速度均存在较大波动，这是由于提升至上层的沙尘颗粒受到旋翼诱导速度的影响而被卷入到下洗流当中，而这部分的流场高度不稳定从而导致颗粒受到的曳力急剧变化造成沙尘颗粒速度方向各异。但是，随着悬停高度的增加，上层沙尘颗粒在切向和径向速度并未减小，这是由于更高的悬停高度导致地面上洗速度减小从而沙尘颗粒受到的曳力亦减小，只有粒径更小且数量更少的沙尘颗粒被提升至上层(详见后文分析)，从而更易受到流场影响并具有较高的速度。
|
| 图 10 H=0.67R时不同区域中沙尘颗粒在X和Y方向的平均速度 Fig. 10 Average velocities in X and Y directions of dust particles in different regions at H=0.67R |
|
| 图 11 H=1R时不同区域中沙尘颗粒在X和Y方向的平均速度 Fig. 11 Average velocities in X and Y directions of dust particles in different regions at H=1R |
|
| 图 12 H=1.5R时不同区域中沙尘颗粒在X和Y方向的平均速度 Fig. 12 Average velocities in X and Y directions of dust particles in different regions at H=1.5R |
沙尘颗粒在Z方向上的速度决定了沙尘云的高度分布，下层的沙尘颗粒按照+Z方向速度的大小不同可分为3类：①速度几乎为0，在地面附近难以受到流场影响的沙尘颗粒；②速度较小，只在下层运动而无法运动至上层的沙尘颗粒；③速度较大，能够受流场影响并运动至上层的沙尘颗粒。而第3类沙尘颗粒则是造成沙盲的主因。从图 13中可以看出，当H=0.67R时，上层沙尘颗粒在+Z方向的速度在1 s时刻达到最大，同时-Z方向的速度几乎为0。这表明此时刻并未有沙尘颗粒被卷入到下洗流当中，而在之后的5 s内，+Z方向的平均速度开始急速下降，-Z方向的速度逐渐增加。这是由于在这段时间内，第一，沙尘颗粒是瞬间同时生成的，从而导致进入上层的沙尘颗粒后续难以得到持续补充；其次，这些已经到达上层的沙尘颗粒开始逐渐被卷入到下洗流当中，产生了向下的速度，而在这之后的时间，随着沙尘云范围的逐步扩大，流场对沙尘颗粒的影响逐渐减弱造成+Z方向的速度逐渐减小。下层粒径(沙尘颗粒直径)较大的沙尘颗粒大部分集中在地表附近，难以受到流场影响而被提+Z和-Z方向平均速度均很小。随着悬停高度的升高，上层沙尘颗粒+Z方向速度达到最大的时刻也随之推迟，这表明更高悬停高度导致沙尘颗粒需要运动更长时间才能到达上层并最终形成沙盲。
|
| 图 13 不同悬停高度时沙尘颗粒在Z方向上的平均速度 Fig. 13 Average velocity of dust particles in Z direction at different hovering heights |
位于上层的沙尘颗粒的数量按照粒径的不同而差别巨大。从图 14中可以看出，当H=0.67R时，达到上层粒径为5 μm的沙尘颗粒数量比粒径为10 μm的沙尘颗粒数量要多。这是因为粒径更小的沙尘颗粒更易在上洗区域被输运至上层，但数据表明位于上层区域这两种粒径的沙尘颗粒数量接近，流场对于这两种尺寸的沙尘颗粒作用差异并不敏感。但是，在该区域内，50 μm和100 μm粒径的沙尘颗粒数量与之相比差异较大，这表明重量较重的沙尘颗粒难以在上层聚集导致浓度较低。从图 14中还能看出，沙盲中绝大多数沙尘颗粒并不能到达旋翼上方，只是在地表附近做平面运动。这说明即使悬停高度较低，对于粒径为50 μm和100 μm的沙尘颗粒而言，也只有少部分能够形成沙盲。当悬停高度增加至1R时，如图 15所示，位于上层的沙尘颗粒数量和H=0.67R时的接近。这可能是由于在该悬停高度下，沙尘云范围比同时刻H=0.67R时候的更小，离旋翼更近的流场能为沙尘颗粒提供的曳力更大，导致即使悬停高度更高，上层的沙尘颗粒浓度却与H=0.67R时的接近。但是，当H=1.5R时，如图 16所示，桨尖涡在流场中需要经过更长的距离才能到达地面，从而对沙尘颗粒的曳力减小，导致难以输运更多沙尘颗粒至上层，从而使直升机逐渐脱离沙盲。
|
| 图 14 H=0.67R时不同粒径沙尘颗粒数量分布 Fig. 14 Quantity distribution of dust particles withdifferent sizes at H=0.67R |
|
| 图 15 H=1R时不同粒径沙尘颗粒数量分布 Fig. 15 Quantity distribution of dust particles withdifferent sizes at H=1R |
|
| 图 16 H=1.5R时不同粒径沙尘颗粒数量分布 Fig. 16 Quantity distribution of dust particles withdifferent sizes at H=1.5R |
从图 17~图 19中可以通过观察每一粒沙尘颗粒在沙尘云中的微观状态，直观地分析沙尘云发展至1 s和20 s时刻的宏观全貌，当H=0.67R时的1 s时刻，地面附近的沙尘颗粒从旋翼中心向外呈现辐射状扩散，集中在速度高的环形区域内，速度大小分布从内到外逐渐减小。从图 17(a)的前视图和侧视图中观察到，在1 s时刻距旋翼中心约3R处，沙尘颗粒开始被提升且有少量沙尘颗粒达到上层。随着时间的发展，沙尘云范围不断扩大，上层沙尘颗粒数量逐渐增多，当沙尘云发展至20 s时刻，从图 17(b)中可以明显看出位于旋翼中心附近部分上层的沙尘颗粒被卷入下洗流当中，从而造成沙尘颗粒在上、下层之间循环往复的运动。随着悬停高度增加至1R，在同样的1 s时刻，从图 18(a)中的前视和侧视图观察到，沙尘颗粒开始被提升的位置离旋翼中心约2R处，离旋翼中心距离比H=0.67R时更近，在此时刻，当H=1.5R时，地表附近的沙尘颗粒尚未明显聚集，但是速度由内到外逐渐减小。在20 s时刻，对比图 17(b)、图 18(b)和图 19(b)可以看出：随着悬停高度的增加，沙尘云范围以及在外层的沙尘颗粒速度也均逐渐减小，速度较大的沙尘颗粒受到下洗流影响而集中在旋翼附近。这表明随着时间的发展，沙尘云中外层的沙尘颗粒会趋于悬浮在空中。
|
| 图 17 H=0.67R时沙尘颗粒的宏观分布图 Fig. 17 Macroscopic distribution of dust particles at H=0.67R |
|
| 图 18 H=1R时沙尘颗粒的宏观分布图 Fig. 18 Macroscopic distribution of dust particles at H=1R |
|
| 图 19 H=1.5R时沙尘颗粒的宏观分布图 Fig. 19 Macroscopic distribution of dust particles atH=1.5R |
此外，需要说明的是，由于沙尘颗粒在流场中受到重力和曳力等的影响，且沙尘颗粒粒径十分微小，会导致有少量的沙尘颗粒在单独的时间步内，穿过物面边界而脱离计算域，从而出现图 17~图 19中部分沙尘颗粒在地面之下的情形。相比位于流场中的沙尘颗粒而言，脱离计算域(地面之下)的沙尘颗粒数量很少，并不会影响计算结果而可以忽略。
4 结论本研究构建了用于直升机地效飞行沙盲现象模拟的CFD-DEM方法。采用该方法针对全尺寸直升机在3个不同高度下地效悬停并发生沙盲时，所产生的沙尘云中的沙尘颗粒进行了受力分析。求解了它们在悬停流场中不同时刻的运动状态和分布规律，得出了以下结论：
1) 沙尘云中大部分沙尘颗粒并未受流场的影响而被提升至上层区域，只在地表附近扩散。沙尘云径向范围随着悬停高度的升高而逐渐减小，但轴向范围变化不大。
2) 沙盲现象发生时，沙尘云中下层的沙尘颗粒主要做径向运动，而切向运动较弱，这种趋势随着悬停高度的增加变得更加显著。上层区域沙尘颗粒在桨盘平面X和Y方向上的平均速度波动较大，且无明显规律可循。
3) 沙尘云中初始到达上层区域的沙尘颗粒在轴向只有+Z方向的速度，随后受到下洗流的影响出现-Z方向的速度，随着时间的发展，+Z方向速度急剧减小，同时-Z方向速度先增加而后逐渐稳定。
4) 沙尘云形成初期阶段，地面处的沙尘颗粒由内向外扩散，随着时间的发展，沙尘云内层的沙尘颗粒仍然具有较大速度，而外层的沙尘颗粒速度逐渐减小，并随着沙尘云范围逐渐扩大而趋于悬浮在上层区域。
5) 位于沙尘云内层空间(旋翼及机身周围)的沙尘浓度较低，位于沙尘云外层空间的沙尘浓度较高。
| [1] | POLZIN J W, GUNTUPALLI K, RAJAGOPALAN R G. Discrete blade model for rotorcraft brownout[C]//29th AIAA Applied Aerodynamics Conference. Reston, VA: AIAA, 2011. |
| [2] | SYAL M, LEISHMAN J G. Comparisons of predicted brownout dust clouds with photogrammetry measurements[C]//67th Annual Forum of the American Helicopter Society, 2011. |
| [3] | SYAL M, LEISHMAN J G. Modeling of bombardment ejections in the rotorcraft brownout problem[J]. AIAA Journal, 2013, 51(4): 849-866. |
| Click to display the text | |
| [4] | WONG O D, TANNER P E. Photogrammetric measurements of an EH-60L brownout cloud[C]//American Helicopter Society 66th Annual Forum Proceedings, 2010. |
| [5] |
叶靓, 招启军, 徐国华. 基于非结构嵌套网格方法的旋翼地面效应数值模拟[J]. 航空学报, 2009, 30(5): 780-786. YE L, ZHAO Q J, XU G H. Numerical simulation on flowfield of rotor in ground effect based on unstructured embedded grid method[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(5): 780-786. (in Chinese) |
| Cited By in Cnki (13) | Click to display the text | |
| [6] |
朱明勇, 招启军, 王博. 基于CFD和混合配平算法的直升机旋翼地面效应模拟[J]. 航空学报, 2016, 37(8): 2539-2551. ZHU M Y, ZHAO Q J, WANG B. Simulation of helicopter rotor in ground effect based on CFD method and hybrid trim algorithm[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2539-2551. (in Chinese) |
| Cited By in Cnki (7) | Click to display the text | |
| [7] | GANESH B, KOMERATH N. Study of ground vortex structure of rotorcraft in ground effect at low advance ratios[C]//24th AIAA Applied Aerodynamics Conference. Reston, VA: AIAA, 2006. |
| [8] | LEE T E, LEISHMAN J G, RAMASAMY M. Fluid dynamics of interacting blade tip vortices with a ground plane[J]. Journal of the American Helicopter Society, 2010, 55(2): 022005. |
| Click to display the text | |
| [9] | MILLUZZO J I. Effects of blade tip shape on rotor in-ground-effect aerodynamics[D]. Maryland: University of Maryland, 2012: 39-83. |
| [10] | HANCE B T. Effects of body shapes on rotor in-ground-effect aerodynamics[D]. Maryland: University of Maryland, 2012: 42-66. |
| [11] | WHITEHOUSE G R, WACHSPRESS D A, QUACKENBUSH T R, et al. Exploring aerodynamic methods for mitigating brownout[C]//American Helicopter Society 65th Annual Forum, 2009. |
| [12] | WHITEHOUSE G R, WACHSPRESS D A, QUACKENBUSH T R. Aerodynamic design of helicopter rotors for reduced brownout[C]//International Powered Lift Conference, 2010. |
| [13] | THOMAS S, AMIRAUX M, BAEDER J D. Modeling the two-phase flowfield beneath a hovering rotor on graphics processing units using a FVM-RANS hybrid methodology[C]//21st AIAA Computational Fluid Dynamics Conference. Reston, VA: AIAA, 2013. |
| [14] | YU H M, CHENG W M, WU L R, et al. Mechanisms of dust diffuse pollution under forced-exhaust ventilation in fully-mechanized excavation faces by CFD-DEM[J]. Powder Technology, 2017, 317: 31-47. |
| Click to display the text | |
| [15] | HOU X Y, DING T X, DENG Z Q, et al. Study of the creeping of irregularly shaped Martian dust particles based on DEM-CFD[J]. Powder Technology, 2018, 328: 184-198. |
| Click to display the text | |
| [16] |
魏可可, 高霄鹏. 基于STAR-CCM+对5415船模的阻力预报[J]. 兵器装备工程学报, 2016, 37(9): 157-161. WEI K K, GAO X P. Resistance prediction of 5415 ship model based on STAR-CCM+[J]. Journal of Ordnance Equipment Engineering, 2016, 37(9): 157-161. (in Chinese) |
| Cited By in Cnki (5) | Click to display the text | |
| [17] | DUBUC L, CANTARITI F, WOODGATE M, et al. Solution of the unsteady Euler equations using an implicit dual-time method[J]. AIAA Journal, 1998, 36(8): 1417-1424. |
| Click to display the text | |
| [18] | QU Q L, WANG W, LIU P Q, et al. Airfoil aerodynamics in ground effect for wide range of angles of attack[J]. AIAA Journal, 2015, 53(4): 1048-1061. |
| Click to display the text | |
| [19] | QIN Y P, LIU P Q, QU Q L, et al. Numerical study of aerodynamic forces and flow physics of a delta wing in dynamic ground effect[J]. Aerospace Science & Technology, 2016, 51: 152-221. |
| Click to display the text | |
| [20] | QIN Y P, LIU P Q, QU Q L, et al. Wing/canard interference of a close-coupled canard configuration in static ground effect[J]. Aerospace Science & Technology, 2017, 69: 60-75. |
| Click to display the text | |
| [21] | KUMAR M, MURTHY V. Analysis of flow around multibladed rotor using CFD in the frequency region[C]//25th AIAA Applied Aerodynamics Conference. Reston, VA: AIAA, 2007. |
| [22] | NOROUZI H R, ZARGHAMI R, SOTUDEH-GHAREBAGH R, et al. Coupled CFD-DEM modeling (formulation, implementation and application to multiphase flows)CFD-DEM applications to multiphase flow[M]. . |
| [23] |
杜俊.基于CFD-DEM方法的希相气力输送数值模拟研究[D].武汉: 武汉大学, 2015: 20-21. DU J. Simulation of dilute pneumatic conveying by CFD-DEM[D]. Wuhan: Wuhan University, 2015: 20-21 (in Chinese). |
| [24] | LIGHT J S. Tip vortex geometry of a hovering helicopter rotor in ground effect[J]. Journal of the American Helicopter Society, 1993, 38(2): 34-42. |
| Click to display the text | |
| [25] | KUTZ B M, KOWARSCH U, KEßLER M, et al. Numerical investigation of helicopter rotors in ground effect[C]//30th AIAA Applied Aerodynamics Conference. Reston, VA: AIAA, 2012. |