对于分布式电推进技术这种存在分布式动力单元的布局形式，如何高效且精确地模拟多动力单元同机体耦合的流场是需要解决的问题之一。

旋翼/螺旋桨等动力单元的数值模拟方法主要有多重参考系法^[1](multiple reference frames, MRF)及基于滑移网格^[2]和嵌套网格^[3]的非定常数值模拟。但采用真实模型的数值模拟均占用较多的计算资源和时间，不适用于多动力的计算。因此人们对复杂动力进行简化处理，激励盘法和动量源法(momentum source method, MSM)被提出。

激励盘方法是将旋翼等效成一个无厚度的圆盘，通过改变将盘前后的压力及切向速度来模拟真实流动。国外方面，德国宇航院^[4]应用激励盘方法实现对直升机前飞状态的模拟；Veldhuis等^[5-6]将其应用于螺旋桨的模拟，研究了螺旋桨滑流对机翼展向载荷分布的影响及螺旋桨处于不同相对位置时的作用。国内方面，左岁寒、李博、曹栋等^[7-9]分别应用激励盘法对螺旋桨滑流及旋翼尾流进行了模拟。

与激励盘法不同，动量源法是将桨叶的旋转效果时均等效到一个薄体圆盘上，根据叶素理论求得桨盘力分布，再将其转化为动量源项加入Navier-Stokes(N-S)方程求解。国外Rajagopalan、Zori、Chaffin等^[10-12]对该方法进行了研究和发展，O′Brien^[13]则研究了桨盘不同力分布模型对直升机前飞状态模拟的影响，并指出在求解不可压N-S方程时，动量源法在数值上更加稳定。国内方面，康宁、童自立等^[14-15]将其用于旋翼的流场求解，宋长虹^[16]采用动量源法模拟了涵道螺旋桨流场及气动特性，成宝峰^[17]则将其用于模拟倾转旋翼不同运动状态下飞行器各部件间的干扰。

动量源法多采用定常计算的方式，但定常计算并不能模拟出桨叶真实旋转下的非定常特征。为解决这一问题，并发挥动量源法计算高效的特点，非定常动量源法(unsteady momentum source method, UMSM)被提出。

国内外学者对非定常动量源法的研究不多，主要有Kim^[18-19]、Guntupalli^[20-21]及南航招启军^[22]、朱秋娴等^[23]。Kim以局部网格区域代替真实桨叶，采用动态入流模型及局部平均获取桨盘速度，结合叶素理论建立非定常动量源法。而Guntupalli提出的方法是将动量源添加到桨叶单位时间步所扫过的扇形网格区域中，扇形区域与真实桨叶无关。2种方法的作用网格区域均随时间的推进而改变，以此模拟真实桨叶的旋转。

国内南航的研究者则采用嵌套网格思想，从背景网格中提取速度，将沿弦向的动量分布，添加到对应的网格单元中，得到的流场结果更加准确。但其仍涉及网格运动，且在求解二维翼型流场及处理嵌套边界信息上需耗费一定时间，相比于前2种方法计算效率降低。

本文给定力分布作为输入，以旋翼悬停状态为例，对定常动量源法进行验证，并提出一种能够较好地模拟桨尖涡发展的非定常动量源法，探究动量源法中影响尾流模拟的主要因素，为进一步研究多动力耦合机体的流场计算做好准备。

1 基于真实桨叶几何的数值模拟 1.1 旋翼模型

计算采用Caradonna-Tung旋翼模型^[24]。该旋翼由2片等弦长、无扭转的桨叶组成，桨盘半径R=1.143 m, 弦长c₀=0.190 5 m, 桨叶剖面为NACA0012翼型。计算忽略桨毂和连接杆的干扰, 忽略桨叶挥舞效果, 无周期变距。

1.2 计算状态

来流压力P=103 027 Pa, 温度T₀=289.75 K, 旋翼转速ω=1 250 r/min。分别采用多重参考系法及嵌套网格非定常模拟2种方法进行计算, 计算基于理想气体、S-A湍流模型。其中, 非定常时间步大小取为Δt=0.000 266 7 s, 即180个时间步为一个旋转周期, 子迭代步数取20。旋转10个周期后, 取3个周期的时均结果作为最终结果进行分析。

1.3 结果分析

旋翼拉力、扭矩系数定义如下, ρ为来流密度

(1)

(2)

从表 1可看出, 基于真实模型的2种模拟方法对旋翼拉力的计算有着较好的精度, 并且均能够得到较为准确的桨叶力分布, 如图 2所示。

2 基于定常动量源法的数值模拟 2.1 定常动量源法简介

传统的动量源法将叶素理论与CFD求解结合, 基于定常计算。基本步骤为:从CFD计算流场中获取叶素截面的当地来流速度, 根据已知叶素截面的几何信息, 求得该截面上的力, 再将其转化为动量源参与CFD求解, 重复进行这个过程, 直至流场收敛。

设单个桨叶在一个周期内旋转用时为t₀

(3)

对于桨盘中的网格单元, 假设其所占圆周角度为dϕ, 则单个桨叶旋转dϕ用时为t_dϕ

(4)

将对应叶素截面作用力dF在一个周期内的作用效果等效到该网格单元上

(5)

对于N片桨叶, 经过等效处理, 便得到动量源S, 式中V_c为网格单元的体积

(6)

通过坐标变换得到S_x, S_y, S_z, 最后将其加入到对应网格单元控制体的动量方程中, 参与迭代求解。

叶素理论在整个过程中实际扮演求取桨叶径向力分布的角色, 如若有准确的力分布作为输入, 那么直接将其转化为动量源参与求解, 可省去力分布求解过程, 加速计算, 前人的研究^{[4, 13]}表明该方法可取得较好的效果。

对于来流均匀的轴流状态(悬停属于该状态的特殊情况), 桨盘不同相位处的力分布相同, 因此只需输入任意相位下的力分布结果即可。本文将拉力、扭矩分布作为给定输入建立定常动量源法, 命名为MSM_g, 计算流程如图 3所示。

2.2 计算状态

计算状态同上。采用薄体网格圆盘等效真实旋翼, 厚度d=0.01 m, 如图 4所示。桨盘周向网格采用均布处理, 并对尾流区域进行加密, 计算网格数量44万。

2.3 桨叶拉力分布对比

动量源法的关键便是取得桨叶径向力分布, 获取桨叶力分布的前处理可采用求解效率较高的MRF方法或叶素理论。对比2种方法得到的力分布结果, 如图 5所示。两者求解得到的力分布趋势相同, 只是其值在靠近桨根处偏小, 而桨尖处偏大。两者的差距是因为, 叶素理论通过xfoil快速计算得到的翼型气动数据进行叶素气动力的计算, 且采用经验公式处理奖尖附近的三维流动, 同CFD求解存在一定差距。

2.4 尾流速度对比

将采用MRF求解作为前处理的方法, 命名为MSM_g1, 采用叶素理论求解的方法, 命名为MSM_g2。而嵌套网格非定常模拟能够较好地反映流场非定常特征, 因此以其结果对动量源法进行校验。

对比桨盘后0.175R距离处尾流速度分布, 如图 6所示。2种定常动量源法所得结果与嵌套网格非定常模拟时均结果在变化趋势上基本相同, 而相比之下MSM_g1结果符合更好。

而从图 6a)中可以看出, MSM_g2相比于MSM_g1所得轴向速度在桨根处偏小, 而在桨尖处偏大, 恰好与两者拉力分布的差距相对应。这便说明了准确的桨叶力分布输入是精确模拟尾流的关键。

此外, 在模拟出的尾流中出现了较为准确的径向流动, 如图 6c)所示, 这是因为动量源法在CFD计算的流场中获取速度信息进而计算气动力, 而流场是连续且相互影响的, 更加接近真实状态。而叶素理论将叶素单独处理, 忽略了叶素间的相互影响。这说明动量源方法能够弥补叶素理论的不足。

从图 7中可以看出, 其能够反映出桨盘后气流的收缩和旋转, 对旋翼悬停流场的模拟取得了较好的效果。

2.5 计算效率对比

计算设备采用Inter Core i7 2.8 GHz四核处理器, 32G内存。对比不同方法的计算效率, 如表 2所示, 不考虑前处理时, 定常动量源法相比于真实模型数值模拟计算效率更高, 而以叶素理论作为前处理的方法则能够真正发挥动量源法的高效优势。

此处2种基于真实模型的计算与表 1相同, 且两者网格量的差距是由于多重参考系法需要划分桨叶所在的柱状旋转域的网格, 而嵌套网格只需单独划分桨叶网格。因此在桨叶网格布点数目基本一致的情况下, 嵌套网格方法所用的网格量更少。

3 基于非定常动量源法的数值模拟

由于定常动量源法不能模拟旋翼真实的非定常特性, 故其在研究动力单元与机翼等部件耦合的流场时, 并不能模拟出涡系之间的相互影响, 精度不及真实模型非定常模拟。为解决这一问题, 并发挥动量源法本身的计算优势, 本文在定常动量源法的基础上提出一种非定常动量源法。

3.1 方法介绍

Kim及Guntupalli提出的非定常动量源法均耦合叶素理论求解桨叶力分布, 且均采用固定的桨盘网格, 如图 8所示。

每个时间步添加动量源的网格区域随着时间推进而变化, 以此模拟真实的旋转运动。其中, Kim的方法是以与桨叶大小相接近的网格区域代替真实桨叶, 如图 8a)所示, 并考虑力沿弦向的变化, 将力分布转化为动量源分布添加到剖面对应的网格中。但其选取的网格区域边界不光滑, 在桨盘一周布置网格较少的情况下尤为突出, 且此时桨叶剖面对应的弦向网格个数较少, 不能很好地体现力分布的作用。因此, Kim采用的方法需在周向布置大量网格, 这样便降低了计算效率。

Guntupalli提出的非定常动量源法是将叶素上的力通过(7)式转化为动量源, 添加到桨叶Δt时间扫过的ω·Δt大小的扇形中, 添加动量源的网格区域大小不固定, 如图 8b)所示。

(7)

基于这种思想, 如果将时间步取为单个桨叶旋转一周的时间, 即ω·Δt=2π, 则上式与定常动量源(4)式在桨叶数为1时相同, 故上文中定常动量源方法可看作是该非定常动量源法的一种特殊状态。但该方法前后2个时间步选取的网格区域不重叠, 流场连续性及收敛性较差, 且时间步不易取大, 理想状态为每个时间步扫过一个网格单元。

以上2种方法均需保证每个时间步推进整数倍的网格单元, 不应在添加动量源的网格区域中存在网格单元不完整的现象。虽出发点不同, 但都保证了一个周期内添加到网格单元的总动量相同, 殊途同归。

本文基于Kim的思想, 但采用固定大小的扇形区域来等效真实桨叶, 避免了网格区域边界的不光滑, 并保持了流场连续性及收敛性。在图 3定常动量源流程基础上, 建立非定常动量源法。计算流程如图 9所示。

3.2 桨叶几何因素影响分析

以MRF取得的力分布作为给定输入。用扇形网格区域代替真实桨叶, 则需考虑桨叶几何因素带来的影响:桨叶大小决定扇形区域的大小; 需考虑弦向力的分布形式。

对于真实的桨叶, 桨叶截面力沿弦向的分布不均, 主要体现在升力之上。根据薄翼理论^[25], 对于对称翼型, 压力系数沿弦向的分布为

(8)

升力系数

(9)

由于压差是升力产生的主要原因, 因此沿弦向的升力分布函数可写作

(10)

采用下式(其中c为所处剖面弦长)对弦向力分布进行处理

(11)

将分布函数f应用于拉力分布, 扭矩分布采用平均处理, 则网格单元的受力如(12)式所示

(12)

(13)

式中，x_c1, x_c2为网格单元在当地截面弦长中所处的始末位置。得到对应网格单元上的作用力dF后, 将其转化为非定常动量源项

(14)

为分析代替桨叶的扇形面积大小变化及不同弦向分布对尾流的影响, 以MRF模拟所得力分布作为给定输入, 采用以下3种方法:

1) 取周向5个网格的扇形区域代替桨叶, 拉力采用f分布, 扭矩均分, 命名为UMSM_f;

2) 取周向5个网格的扇形区域代替桨叶, 拉力、扭矩均采用均分方式, 命名为UMSM_5;

3) 取周向2个网格的扇形区域代替桨叶, 拉力、扭矩均采用均分方式, 命名为UMSM_2。

时间步取为桨叶扫过一个网格单元所用时间, 由于周向网格数目为180, 因此非定常时间步大小亦取为Δt=0.000 266 7 s, 由于气动力给定不需要迭代求解, 收敛相对较快, 因此子迭代步数取10, 旋转5个周期后, 取3个周期的时均结果作为最终结果进行分析。非定常模拟基于MSM_g1的定常收敛结果进行。

对比图 10真实模型非定常模拟与图 11 UMSM计算结果, UMSM所得轴向速度分布整体趋势上与真实模型非定常模拟结果符合较好。两者间的差距主要在桨叶附近, 而在尾流区中随着远离桨盘两者的差异逐渐减小。这说明了桨叶几何外形所带来的影响主要集中在桨盘附近。这也间接验证了动量源法简化真实模型的可行性。

从图 11中可以看出, 3种UMSM结果间的差异亦表现在桨盘附近, 在桨盘后0.087R处三者结果已趋于一致。

进一步对比尾流时均轴向速度分布, 如图 12所示, 三者结果几乎相同, 均与真实模型非定常模拟结果符合较好。这便说明了桨叶几何外形对尾流的影响主要是通过决定桨叶力分布实现, 在拥有准确力分布的前提下, 替代桨叶的扇形大小及弦向力分布形式等由桨叶几何外形所决定的参数对整个尾流速度影响较小, 主要集中在桨盘后小范围内。

建议扇形区域的选取以桨叶平均弦长为参考, 使两者平均弦长接近; 为保证流场连续性较好, 扇形区域周向至少布置2个网格, 时间步的选取建议一个时间步推进一个网格; 而为提高适用性, 弦向拉力、力矩分布均采用均布处理。

3.3 2种非定常动量源法结果对比

采用叶素理论求得的力分布作为输入, 同样取周向2个网格的扇形区域代替桨叶, 拉力、扭矩均采用均分方式, 命名为UMSM 1。将其与UMSM 2进行对比, 分析不同前处理方式应用于非定常动量源法中的特点。

两者盘后轴向速度分布如图 13所示, 其与上文定常动量源法的对比结果基本相同。

从图 14可以看出, 本文所建非定常动量源法能够模拟桨尖涡的生成和发展。进一步对比桨尖涡下滑位置变化, 如图 15所示, 2种方法均取得了较好的模拟结果。

3.4 不同转速下方法验证

从图 16桨尖涡对比中可以看出, 不同转速下本文所建非定常动量源法对桨尖涡的模拟依然有较好的模拟精度。

3.5 计算效率对比

对比不同方法的计算效率, 如下表所示。非定常动量源法计算所用网格数量远小于真实模型非定常模拟, 即使考虑采用MRF进行前处理, 用时也只有其1/8左右, 具有更高的计算效率。

4 结论

本文以桨叶力分布作为输入，建立了高效求解旋翼悬停流场的定常及非定常动量源法，并验证了这种动力单元简化方式的可行性，得到的结论如下：

1) 对于动量源法，准确的桨叶力分布是精确模拟尾流的关键，桨叶几何模型对尾流的影响主要体现在桨叶力分布之上。

2) 对比定常、非定常2种动量源法，文中建立的非定常动量源法对桨盘后的尾流模拟精度更高，且能够较好地模拟桨尖涡的生成和发展。

3) 采用MRF获取力分布的非定常动量源法，获得的流场更加准确，可用于较为精确的流场模拟；虽前处理需消耗一定的时间，但效率依然高于嵌套网格非定常模拟，且在求解多动力相关流场时，只需计算单独动力单元，能够显著提高求解效率。

4) 采用叶素理论求解力分布的非定常动量源法，可省去真实模型模拟的过程，效率更高，可用于耦合多动力流场快速及初步的计算。

文中以旋翼悬停为例，进行了非定常动量源方法的验证，而旋翼悬停属于均匀轴流状态，但在非轴流状态，桨叶力分布是变化的，需获得不同相位的力分布，才能使该方法正常进行，此时力分布不能通过MRF较快地给出。因此，如何处理非轴流状态的计算是接下来的研究重点。