WO2013078628A1

WO2013078628A1 - 直升机旋翼飞行结冰的数值模拟方法

Info

Publication number: WO2013078628A1
Application number: PCT/CN2011/083189
Authority: WO
Inventors: 路明
Original assignee: 天津空中代码工程应用软件开发有限公司
Priority date: 2011-11-30
Filing date: 2011-11-30
Publication date: 2013-06-06
Also published as: US20140257770A1

Abstract

一种用于模拟直升机旋翼的飞行结冰的数值方法，其包括在旋翼尾迹中、旋转坐标系下在用于模拟空气-超冷水滴两相流的旋涡运动的单流体模型中的动量方程和能量方程中加入涡量补偿力的算法、在超冷水滴的滑移速度方程中加入离心力和科里奥利力的算法；考虑离心力和科里奥利力的水膜运动和结冰模型；利用上述算法和模型进行直升机旋翼飞行结冰数值模拟计算的流程。其中模拟空气-超冷水滴的运动的单流体两相流流动的模型是一组描述旋翼尾迹流场中空气和超冷水滴二元混合物运动的偏微分方程组，加上一个两相流滑移速度模型。旋翼表面上的水膜运动和结冰过程的模型中考虑了因旋翼的旋转运动使水膜和冰层内产生的质量和能量的变化，以获得旋翼表面的结冰形状。

Description

直升机旋翼飞行结冰的数值模拟方法技术领域

本发明涉及航空工程领域，是计算流体力学在航空工程领域的应用。具体地讲，是一种用于模拟直升机旋翼的飞行结冰的数值方法。这种数值模拟技术可以用计算机高级程序语言实现，并通过计算机运行计算机程序来模拟直升机旋翼在空中飞行期间遭遇结冰时的状态。

说背景技术

直升机以其良好的机动性和灵活性，在军事、书民用等各个领域有着广泛的应用。直升机的飞行高度多处于 6000米以下，在此高度范围内，大气中的超冷液态水滴（温度在冰点以下，以液态水滴形式存在）会碰撞在飞行器表面的部件上形成水膜，如机翼、机身、驾驶舱、尾翼、发动机进气口等部件的表面都容易形成大量积聚的液态水膜。如果超冷液态水含量（单位体积内的超冷液态水滴的总重量，量纲 kg/m³) 较高，积聚的水膜会在飞行器件表面上结成冰层。这种现象被称为飞行结冰。和其它任何类型的飞行器一样，表面的结冰极大影响飞行器的空气动力学特性，严重的会引发事故和灾难。就直升机来讲，直升机的主旋翼具有多种用途,它在结构上同时起到机翼、升降舵、方向舵以及推进装置的作用。旋翼表面结冰将直接威胁直升机的安全。首先，旋翼上的结冰会使旋翼上的负载增加，引起输出扭矩的增加。结冰后，旋翼形状的改变会导致机翼失速特性的恶化。旋翼上冰块的脱落又会使负载变得不均勾，引起直升机更大的振动和噪声。这些不利因素都表明直升机的防冰尤为重要。而主旋翼的防冰是直升机防冰最核心问题，也是最复杂的问题。与固定的机翼相比，旋转的机翼要与前一个机翼的尾迹有干涉现象。此外，旋转的机翼会产生离心力，这些复杂的空气动力学现象在不稳定运动（例如，机翼的震颤）情况下会变得更加复杂。

现代的飞行器都装有除冰系统，而且在仪表规则 IFR ( Instrument Fl ight Rules ) 条件下使用，利用电加热的方法对飞行器在飞行期间进行除冰和防冰。电加热设备的启动是通过系统对安置在飞行器外表面各处的结冰传感器的反馈控制完成的。

为保证飞行器在遭遇结冰条件下飞行的安全性，飞行器制造者必须表明飞行器能在各种结冰飞行条件下满足飞行包线以取得适航证。适航取证的过程是通过空中飞行测试、风洞测试、计算机数值模拟三种手段共同实现的。空中飞行是最直接的测试手段，完全在自然条件下进行。但是这种方法不仅成本昂贵，而且自然条件无法完全达到飞行包线上的所有条件，不可能进行逐个条件下的验证。在陆地上唯一替代空中飞行测试方法就是倚重结冰风洞，在风洞中产生高空飞行时遇到的超冷水含量和水滴尺度等结冰条件。结冰风洞制造成本昂贵，而且风洞中无法产生的大气中超冷水滴尺度的分布。由于风洞中采用缩小的飞行器几何模型，所以流动雷诺数不准确，以致难以准确预报真实的结冰状态。

自从上世纪 80年代，国际上逐步开始了飞行器飞行结冰过程的数值模拟的研究，不仅在理论上对这一复杂的现象有了深刻的认识，而且已经将数值模拟的结果应用到工程上。飞行结冰的计算机数值模拟已经成为一种新产品开发设计和适航取证的支持工具。先进的数值模拟技术还可以弥补结冰风洞实验中的误差和缺陷。飞行结冰过程的数值模拟技术经过 20多年的发展后，出现了具有代表性的、具有工业应用价值的相关产品。例如，来自美国的 LEWICE软件、来自法国的 0NERA软件、来自加拿大的 Fensap软件。这些软件集中体现了飞行器飞行结冰过程的数值模拟技术现状。

飞行器飞行结冰的现代的数值模拟技术中所采用的流程都是基于三个主模块的调用。三个主模块及其各自主要功能是分别为- 空气流动模块：用于求解飞行器外部流场（包括飞行器表面）信息，即求解空气运动的控制方程；

超冷水滴运动模块：用于求解大气中的超冷水滴与飞行器表面的碰撞过程，以获得飞行器表面的液态水的状态（用液态水收集率表示），即求解水滴运动方程；

结冰状态模块：用于求解液态水膜的结冰过程，获得结冰后的几何形状。目前，这些软件中使用的数值模拟技术都是为固定机翼的飞行器开发的，没有考虑直升机旋翼运动的复杂性以及对结冰过程的影响。例如，直升机旋翼的尾迹是以空气-超冷水滴的两相流旋涡运动为主（vortex-dominated ) 的可压缩流场，从前一个旋翼脱落的旋涡与后一个旋翼发生干涉，这种流动现象与固定机翼周围的流场不同。流场的旋涡运动的特点是弱不连续，跨过不连续界面，密度和切向速度不连续，但没有质量传递、径向速度和压力连续。对于这种接触不连续的数值解需要依靠数值方法中固有的数值耗散或者加入人工数值耗散获得，但是数值耗散即使很小也会使弱不连续界面变得模糊，降低数值解对流场的预测精度。此外，超冷水滴在与旋转运动的机翼表面碰撞后生成的水膜会受到离心力和科里奥利效应（Coriol is effect ) 的影响，会在旋翼翼展的方向发生流动，因而其结冰过程也与固定机翼的结冰过程不同。所以，基于直升机旋翼周围流场流动的复杂性，需要设计特殊的数值计算方法，才能对直升机旋翼的飞行结冰过程进行精确的模拟，以达到分析、预测这一复杂物理现象的功能。

首先，由于直升机旋翼旋转运动特点，数值计算采取的计算空间需要按照旋翼角速度建立旋转坐标系，因而会使流体运动的控制方程的原始形态发生变化。图 1 给出了用于直升机旋翼流场计算的坐标系以及速度和力的矢量图。图中表示，空间任意一点 Λ 其位置用 ί=表示，围绕固定旋转轴 Ζ-轴，以角速度 _∞ =[0,0,ω_ζ。]转动。在这个旋转坐标系下必须考虑离心力和科里奥利力的作用。单位质量离心力和单位质量科里奥利力 I 分别定义为：

(1)

-2(ώ_ζο xV), (2) 其中， V是流体在 ⁵点在固定直角坐标系下的运动速度（V = w' + v' + ^，其中 Μ,ν, 是三维直角坐标系的人 j、 A方向上的分量），也是在旋转坐标系下的相对速度。式中符号 X代表叉乘运算。在旋转坐标系下的流体运动控制方程中，将这两个单位质量力乘以流体密度 p，并作为源项加入到固定坐标下的动量和能量方程中。

此外，大气中的超冷水滴的温度低于冰点，统计平均尺度一般在 50 μιη以下，其密度在飞行结冰条件下可以用超冷液态水含量（Z T)表示，一般在 10—³¾ ₃的量级，比较均勾地分布

/ m

在空气中，和空气一起做对流运动，成为一种空气-超冷水滴的两相混合流体。空气-超冷水滴的两相相互作用程度可以通过下面的方法进行比较：

超冷水滴体积分数 α₂定义为 Z T与超冷水滴密度之比， ¾约为水的密度，在飞行器结冰表面，由于超冷水滴的集聚， L¥C变大，极限情况下等于空气的密度，所以有

密度比 ξ定义为超冷水滴密度 ρ₂和空气密度 Α之比，

斯托克斯数（Stokes number) 定义为超冷水滴和空气的豫迟时间之比，实验数据表明，超冷水滴豫迟时间 ₂约[1。-⁵ ~10] 、空气的豫迟时间约[10-³ ^ιο-¹]^ , 所以，

空气受干扰强度定义为

超冷水滴受干扰强度 /,定义为

从以上比较分析的结果可以认为在远离飞行器的区域一远场，空气运动受超冷水滴的影响可以忽略不计（因为《₂较小）；在距离飞行器较近的区域一近场，两相间的相互影响增大，需要考虑空气受超冷水滴对的影响，而超冷水滴受空气的影响始终存在与整个流场。此外，在旋翼的尾迹中除了需要考虑两相的相互影响，还要考虑旋涡运动的特点。根据以上分析，在直升机旋翼飞行结冰的数值模拟中可以将整体计算域分成三个部分：远场、近场和尾迹。图 2给出直升机旋翼数值模拟使用的流场分区策略，按照各个区域的流动特点分别运用不同的流体运动控制方程和物理模型描述空气-超冷水滴的两相运动。远场

空气和超冷水滴运动的控制方程可以分别建立。空气运动因为不受水滴运动的影响，其控制方程组包括旋转坐标系下的连续、动量、能量方程和湍流模型，均为公知的。通过求解空气运动控制方程可以获得远场中的空气速度、密度、压力、温度、动力粘度及其湍流特性。超冷水滴运动的控制方程可以通过颗粒动力学的牛顿定律的原理获得，需要考虑旋转坐标系下的空气的粘性阻力、超冷水滴的重力、离心力和科里奥利力的作用，可以获得超冷水滴在远场的运动速度，其方法是公知的。近场

飞行器外部流场中空气-超冷水滴组成的二元混合物，混合均勾，在一个很小的流体微团内可视为连续介质，就相当于某种等效的单一流体混合物的流动，可以用单流体两相流流动的处理方法建立起流动控制方程。其中，等效混合物的物理性质，如密度、比热、粘性系数、导热系数等都可从两种组分中的对应参数按照体积分数进行加权平均；而混合物的速度按照质量分数进行加权平均。为了体现两相间的相互影响，还需要加入体积分数扩散方程和空气 -超冷水滴的滑移速度方程，目的是为了在必要时将混合物流动变量分解成两相各自独立的流动变量。两相流中的滑移速度是指两相流动速度之差，即超冷水滴速度为与空气运动速度之差。空气-超冷水滴的单流体两相流流动控制方程组包括：连续、动量、能量方程、两相体积分数的扩散方程、超冷水滴滑移速度方程，如果来流是湍流，还需要增加单流体两相流的湍流模型方程。以上方程组在固定坐标系下的形式均是公知的，在旋转坐标系下需要改写，与远场的处理方法一样。

空气-超冷水滴的单流体两相流流动的模拟结果给出近场的流动信息，即各个计算网格点上（或是网格单元内部）的空气和超冷水滴的密度、压力、速度，以及由上述独立变量推导出的温度、动力粘度、湍流特征等信息。近场的结果作为边界条件将被用于尾迹计算、水膜运动和结冰模型中。

发明内容

本发明是一种关于直升机旋翼飞行结冰的数值模拟方法，这种数值模拟技术可以用计算机高级程序语言实现，并通过计算机运行程序来模拟直升机在空中飞行期间遭遇结冰时的状态。该数值模拟方法特点是更接近真实的飞行条件，兼顾计算精度、效率和功能。本发明提出的方法的主要特征是包括在旋翼尾迹中、旋转坐标系下在用于模拟空气 -超冷水滴两相流的旋涡运动的单流体模型中的动量方程和能量方程中加入涡量补偿力的算法、在超冷水滴的滑移速度方程中加入离心力和科里奧利力的算法; 考虑离心力和科里奧利力的水膜运动和结冰模型; 利用上述算法和模型进行直升机旋翼飞行结冰数值模拟计算的流程。其中模拟空气-超冷水滴的运动的单流体两相流流动的模型是一组描述旋翼尾迹流场中空气和超冷水滴二元混合物运动的偏微分方程组，加上一个两相流滑移速度模型，表明在旋涡为主的流动区域需要精确捕捉捕捉旋涡结构的算法。旋翼表面上的水膜运动和结冰过程的模型中需要考虑因旋翼的旋转运动使水膜和冰层内产生的质量和能量的变化，以获得旋翼表面的结冰形状。计算流程需要即保证计算精度，又兼顾计算效率，达到工业实用性的目的。尾迹

两个相邻旋翼之间的区域，其中最大的特征是以旋涡为主的空气-超冷水滴的两相流动。尾迹中两相流旋涡和下一个旋翼由干涉现象，是直升机噪声和振动的主要原因。而旋涡中的超冷水滴和下一个旋翼的碰撞过程直接影响旋翼上的结冰过程。和如前所述，数值方法中固有的数值耗散和人工数值耗散会削弱数值解中对于旋涡一类的弱不连续间断面的捕捉精度。一个用来描述旋流场的变量是涡量（vorticity), 数值耗散会使旋流场的涡量有明显的损失，致使旋涡运动界面的数值解变得模糊。本发明提出一个提高可压缩两相流旋流场的数值模拟精度的方法，其基本原理是人为地补充被耗散的涡量，被补充的涡量作用于每一个计算网格单元，以两相流流体微元的体积力形式体现，这个体积力被称为涡量补偿力。

直角坐标系下的空气 -超冷水滴两相流混合流体运动速度¹^ = uj + vj + w_mk，其中

^，！^，^^是三维直角坐标系的人 j、 A方向上的分量，首先按照涡量的定义，直角坐标系下空气-超冷水滴两相流涡量 _m = ¾^ + ^' + ¾ ^是速度的旋度， 0_m表示为

涡量补偿力是由混合流体的密度 ^乘以单位质量涡量补偿力 i_ffl，定义为

(9) 式中，符号 X表示叉乘运算； fi_ffl代表涡量最大梯度变化的方向；代表涡运动补偿粘性系数，与流体的运动粘度有相同的量纲，是一个标量；代表补偿涡量的特征半径。公式（9) 中， fi„,代表涡量最大梯度变化的方向，即 ν

( 10) 其中， ^是涡量的模；是涡量的模的梯度； |ν^|是涡量的模的梯度的模，分别定义为

+ ,「 δφ

+ ( 13) dx dy

涡量的模的梯度的模表示了最大涡量的变化程度。涡量补偿的空间方向应该是沿着涡量梯度变化的最大方向与涡量耗散的方向的叉乘方向。涡运动补偿粘性系数 ν„,定义为

R

( 14) 其中，二维和三维空间的补偿涡量的特征半径分别定义为

R =丄0 禾 Β R =丄0 ， ( 15) c 2 ^c 2

其中 Ω在二维情况下代表计算网格面积；三维情况下代表计算网格体积尾迹中空气 -超冷水滴运动的数值模拟计算所需要的假设条件为：

1. 空气运动不计空气重力；

2. 空气运动是湍流；

3. 空气运动受超冷水滴运动的影响；

4. 超冷水滴运动受空气阻力的影响；

5. 超冷水滴运动受重力、离心力和科里奥利力、涡量补偿力的影响，不计其他体积力;

6. 不考虑超冷水的湍流脉动；

7. 超冷水滴运动中没有聚合、碰撞、分裂现象发生；

8. 空气和超冷水滴运动的两相间没有相变发生。根据以上分析和假设，直升机旋翼飞行结冰的数值模拟方法中需要在旋转坐标系下建立模拟空气-超冷水滴两相流的旋涡运动的单流体模型，共包括五个偏微分方程组，其中，连续方程和组分扩散方程是公知的，本发明提出的算法是：在动量方程和能量方程将涡量补偿力加入到源项的位置；在超冷水滴的滑移速度方程中加入离心力和科里奥利力。在无粘流的动量方程、能量方程和超冷水滴的滑移速度方程的具体形式如下：

动量方程 dp_my_m

+ V · P_m V_m V_m = -Vp_m - V · \p_mc₂ (1 - c₂ )V,V, + p_m{g + f_centr + f_corio + ( 16) dt

+ · (p_mH_r Y_m ) =— V · (p_mc₂ (1 - c₂ )v,v, )v_m + p_m (i + f_centr + f_∞,₀ + ΐ_ω )v„ ( 17)

超冷水滴的滑移速度方禾 '王

V, ( 18)

在动量方程和能量方程中，涡量补偿力 A„f_ffl被作为源项处理。公式（14) 一 ( 16 ) 中的变量的下标表示混合物、 7表示空气、 2表示超冷水滴；密度、三维速度矢量、压力、温度、时间，分别用 p,V,p,7 表示；速度矢量在笛卡尔坐标方向上的分量分别为 Μ, ν, ; 物性参数,动力粘度系数、定容比热、定压比热、导热系数，分别用 /,C_V , Q , 表示。重力向量表示为 g = gj + gyj + g_Zk；由于使用旋转坐标系，能量方程需要用 MI (roergy) 和转焓 H„ (rothalpy) ， u_m + v_m + w_m ω_οζ x + y

E =C T + (19) u_m +v_m + w_m ω_οζ \x + y

C T + (20) 元混合物中，空气的体积分数为《,、超冷水滴的体积分数为《₂，密度 p_m、物性参数，如混合物的粘性系数 /_m，按照体积分数进行加权平均， p_m - a_xp_x + a₂p_: (21) μ_η = «1 Α + «2^2 (22) 混合物和速度 „按照质量分数进行加权平均

V_m =—— + ₂p₂Y, (23)

P_m ' 并定义 c,为超冷水滴的质量分数

公式（18) 中表示超冷水滴的滑移速度；是托动系数，

其中，超冷水滴单滴的雷诺数

P. V d

Re (26)

离心力科向: .和里奥利力向: 的形式分别是已经由公式（1) 和（2) 给出，需要将^带入其中，在动量方程、能量方程中被处理为源项。此外，如果考虑湍流影响，上述方程采用雷诺平均的方法，并加入湍流模型以封闭方程组。

求解上述方程组可以获得旋翼尾迹中空气、超冷水滴的速度、密度、压力、温度、粘度、湍流特性等流动信息。水膜运动

水膜是超冷水滴碰撞到物体壁面或结冰表面形成的，结冰将首先在水膜和飞行器干净的壁面之间发生，之后将在水膜和已经结成的冰层之间发生。图 3表示旋转坐标系下旋翼上水膜流动的示意图。图中表示，在旋转坐标系下的机翼上建立当地三维正交坐标系坐标 (₀, ,ζ",/7)，其中 ?7-轴方向是壁面法向方向（用向量，也是水膜运动的法向方向和水膜厚度方向、 -轴方向是翼展方向、轴方向是弦长方向。水膜流动的边界包括水膜 /空气 -超冷水滴两相混合物、水膜 /壁面、水膜 /冰层。水膜运动速度是^ =M^ + A^' + ^，其中， ^,^, ^是 ^分别在 ,ζ·,/7)方向上的分量。首先对水膜运动做如下的假设条件-

1. 水膜流动为不可压缩层流流动，并可按照薄膜理论作简化处理；

2. 在水膜-空气界面不考虑水膜对空气的作用，只考虑空气对水膜流动的影响；

3. 不考虑超冷水滴撞击对水膜流动的影响。机翼表面的水膜很薄，水膜运动速度在其流动的法线方向是一个小量，即 <<^, 。同时，在流动方向的变化率是一个小量，即^ 》^^^^。在建立水膜运动的动量方程时，可以假设^ =0，同时消去^^,^^项。此外，为进一步对水膜运动方程进行合理简化，可以用一组典型数据进行水膜运动中各种受力情况的量级分析。直升机旋翼飞行结冰经常遇到的一种情况是：旋翼转速《 = 1200rw7；水膜流动速度 ^ =\mls；水膜厚度 = ΙίΓ⁴^；水膜密度 p_w= Kglmⁱ; 水膜特征长度（机翼弦长） / = 0.1m；机翼翼展长度 = 2m _; 水膜运动粘度 v_w=10- ^{6 2}/s ;水膜表面张力系数 = lN/ 。各种受力的量级分析以粘性力为参考，用各种受力与之进行比较。

表面张力 /粘性力

重力比 /粘性力

ωΒ

离心力 /粘性力 -4x10 ；

zV_f ν,, V,

2p_wcoV_f ~ 2ωδ

科里奥利力 /粘性力: 4x10—'。

2 / ,

按照受力的重要性排序，实际计算中不计与粘性力之比在 10— ³以下的力。根据以上分析，可进一步给出其他假设条件，

4. 不考虑水膜内的压力梯度、惯性力；

5. 不考虑水膜在和空气界面上的表面张力；

6. 水膜流动受重力、离心力和科里奥利力的影响，其他体积力不计。根据以上假设，可以获， ρ_κΐΛΥ ξ,ς,_η θ , (27)

其中， ^ ,ζ·,/;)是水膜的运动速度； /„和 ^分别是水膜的动力粘和密度，是水膜温度 7；的函数。 ϊ^^,^,ζ·,/;)是水膜所受的体积力的合力，是单位质量重力、离心力和科里奥利里之和，并投影到当地坐标系，即

（²⁸) 其中，单位质量重力和离心力 ^^的合力在当地坐标系下的投影 ¾ ，即

(29) 式中代表原始坐标系 (ο,χ, ,ζ)和当地坐标系 (ο,^,ζ·,/;)之间进行坐标变换的张量矩阵，

- cos a sin Θ cos a cos Θ sin a (30)

其中，角度 ^和《见图 3中的定义。同理，可将旋翼在 (ο,χ,^,ζ)坐标系下的旋转角速度 ώ_ΰζ投

ω。ζ,/ = e * ω。_ζ (31) 在当地 (₀, ,^77；)坐标系下的单位质量科里奥利力 ϊ ，按照公式（2) 给出的定义有

(32)

如果水膜厚度是，将公式（28) 至（32) 带入（27)，对其在 |7,/^的范围内进行积分，并利用水膜剪切应力的定义 ϊ„ = A 和假设条件 5 iT = u 和 ^分别

是水膜 /空气-超冷水滴两相混合物接触面上（当 7； = /^时），在两个方向上的切应力向和动力粘度系数），可以求解水膜内高度为 ?7的平面上的速度向量即

上式中，是科里奥利力系数张量矩阵，反应了科里奥利力的影响， κ k₃ 1 -k (34)

其中的元素 [^,^Λ]是向量在 ( ,ζ",/7)方向上的三个分量，向量 S为

OZ (35) μ、

可以定义水膜的平均速度为，

V/ {ξ,ς,'η) =— f Y_f (ξ, ς,η)άη。 (36) 水膜运动的速度或平均速度是水膜运动的动量方程的解，建立动量方程是使用的假设条件在建立水膜结冰过程的模型过程中不再使用。结冰模型

当水膜发生结冰的时候，在水膜 /空气-超冷水滴两相混合物界面上有超冷水滴的收集、水膜的蒸发现象；在水膜 /冰层界面上、水 /壁面上有冰的生成现象，都导致水膜的质量传递。此外，由于旋翼的旋转运动，水膜会受到离心力和科里奥利力的作用，由此也将产生质量的传递。在一定时间内，水膜内部达到质量平衡。同时，与之接触的冰层之间有相变的产生，冰层内部的温度分布也会随之变化。图 4 (a) 表示旋转坐标系下旋翼上水膜结冰过程中的质量平衡的示意图。图中， ^是超冷水滴碰撞到水膜所加入的质量流率、 ₂是水膜表面蒸发 /汽化流出的质量流率； ^是水膜与冰层或者壁面因结冰流出的质量流率。按照质量守恒的原则，得到水膜结冰的连续方程

+ \p_wY_†h_† )= m_l - ( 37) 其中， _V = ^^ + ₊ ；是超冷水滴的体积分数《₂、速度 ₂和水膜温度的函数； '₂ οξ ος οη 是水膜温度 Γ„的函数，其表达式是公知的；图 4 (b) 表示旋转坐标系下旋翼上水膜结冰过程中的能量平衡的示意图。图中， G是超冷水滴碰撞到水膜所加入的能量、是水膜表面蒸发

/冰层升华流出的能量； έ是水膜与冰层或者壁面因结冰流出的能量； ά是壁面流入的热通量。按照能量守恒的原则，得到水膜结冰的能量方程

其中，在旋转坐标系下超冷水滴碰撞到水膜所加入的能 J

其中， c_pw是谁的定压比热、 r_2∞是超冷水滴在运出的温度、是集水率。 ρ₂、 ^和^均是 r_w 的函数，其表达式的形式是公知的。

公式（37)、 ( 38 ) 联立，并加上水膜结冰的正定约束条件（其表达式的形式是公知的），可求得、 Γ„和 ₃，再将转化为冰层厚度，则完成了直升机旋翼飞行结冰在某一个时间间隔内的模拟完成一个时间间隔内的飞行结冰计算后，按照有新的结冰界面形成的计算区域重新生成计算网格，进入下一个时间间隔的计算。图 5表示了本发明提出的直升机旋翼结冰的数值模拟方法计算流程的详细步骤，具体是

(1) 未结冰的旋翼周围生成旋转坐标系和计算网格；

(2) 划分远场、近场和尾迹三个计算区域；

(3) 规定计算的开始时刻 t;

(4) 进行旋翼远场空气-超冷水滴的运动的单流体两相流流动模拟，给出空气运动的流动变量，速度、密度、压力、温度、湍流特性和超冷水滴的运动速度；

(5) 进行旋翼近场和尾迹空气-超冷水滴的运动的单流体两相流流动模拟，给出两相混合流体的压力 p_m、速度 Ϋ_Μ、温度 T_m、动力粘度系数 /_m、剪切力 r_{m ;}

(6) 求旋翼壁面水膜运动的平均速度^

(7) 求水膜厚度/ ? _;

(8) 求结冰量和结冰厚度；

(9) 获得新的冰层形状；

(10) 重新划分计算网格；

(11) 回到步骤（2) ，进行下一个时刻 (t + At)的流体两相流单流体流动模拟计算。

本发明提出的数值方法的流程的特点是：

在旋翼流场的远场，不考虑超冷水滴对空气的影响，分别求解空气和水滴的运动方程；在近场和尾迹中考虑两相间的相互影响，求解空气-超冷水滴的运动的单流体两相流流动的运动方程和两相流滑移速度模型，可以精确地捕捉旋涡结构。旋翼表面上的水膜运动和结冰过程的模型中考虑了因旋翼的旋转运动对水膜的影响，以获得旋翼表面的结冰形状。计算流程即保证了计算精度，又兼顾计算效率，达到工业实用性的目的。

附图说明

图 1 用于直升机旋翼流场计算的坐标系以及速度和力的矢量图

图中， 1:P点的空间位置向量、 2: 旋转轴角速度向量； 3: 卷吸速度、 4: 离心力向量、

5：相对速度向量、 6: 科里奥利力向量。

图 2 直升机旋翼数值模拟使用的流场分区策略图中， 1 : 远场、 2: 尾迹、 3: 近场、 4: 尾迹、 5: 机翼、 6: 近场。

图 3 旋转坐标系下水膜流动示意图和水膜运动坐标系

图中， 1 : 机翼、 2: 水膜。

图 4 (a) 旋翼上水膜结冰过程中的质量平衡示意图。

图 4 (b) 旋翼上水膜结冰过程中的能量平衡示意图。

图 5本发明提出的飞行器飞行结冰的数值模拟方法计算流程图。

图 6剖面为 NACA0012的三维机翼和表面计算网格。

图 7 计算结束时的结冰状态。

具体实施方式

以下以一个具体实施方式进一步说明本发明提出的一种用于直升机旋翼飞行结冰的数值模拟方法，这种数值模拟技术可以用 Fortr_an90计算机高级程序语言实现，并通过计算机运行来模拟横截面为 NACA0012形状的三维旋翼在空中飞行期间遭遇结冰时的状态。

具体实施方案中，直升机采用四个旋转的机翼。展翼比为 7，二维剖面是 NACA0012机翼。计算域取四分之一象限，包含一个机翼， Z-轴为旋转轴。计算网格采用多块结构化六面体网格，机翼周围采用 0-型网格，围绕机翼剖面方向 192个；法线方向 48个；翼展方向 32个。如图 6 所示，剖面为 NACA0012 的三维机翼和表面计算网格。数值模拟的空间离散采用有限体积法 (Finite Volume Method), 每个计算网格成为控制单元，流动变量在控制单元的中心处。

数值模拟的飞行条件是：

来流速度： 57m/s

Z-轴为旋转角速度： 1200rpm

来流温度： 243K

大气压力： lOOkPa

攻角： 0°

液态水含量： 2. 58g/m³

超冷水滴直径： 50 μ ιη

冰的密度： 917kg/m³ 首先按照计算流程将计算域划分为远场、近场和尾迹三个部分。图 6是围绕三维旋翼的计算空间；距离机翼弦长十倍远之外的区域是远场；其余部分为近场和尾迹。尾迹在机翼的上游和下游部分， z-轴方向高度为五倍机翼弦长，被包含在近场之内。

远场

数值模拟计算所需要的假设条件为：

1. 空气运动不受超冷水滴运动的影响；

2. 空气运动不计空气重力；

3. 空气运动是湍流；

4. 超冷水滴运动受空气阻力、重力、离心力和科里奥利力的影响，不计其他体积力

5. 不考虑超冷水的湍流脉动；

6. 超冷水滴运动中没有聚合、碰撞、分裂现象发生；

7. 空气和超冷水滴运动的两相间没有相变发生。按照以上假设，空气和超冷水滴运动的控制方程可以分别建立。空气运动控制方程组包括连续、动量、能量方程，湍流模型采用 Baldwin-Lomax模型，均为公知的。其中，能量方程中需要使用的转内能 E 和转焓 H 的定义为

u_x + j + w₁ ω, + y

H„ oz x

( 41 ) 2

通过求解空气运动控制方程可以获得远场中的空气速度、密度、压力、温度、动力粘度及其湍流特性；通过求解超冷水滴的运动方程可以获得远场中超冷水滴的运动速度。

近场

距离机翼弦长十倍以内的计算域。数值模拟计算所需要的与近场的假设条件不同的是： 1. 空气运动受超冷水滴运动的影响。

按照以上假设，空气运动的控制方程，连续、动量、能量方程和湍流模型均需要重新制定。如前所述，可以用空气-超冷水滴的两相单流体流流动的处理方法建立起流动控制方程。为了体现两相间的相互影响，需要加入超冷水滴的体积分数扩散方程和空气-超冷水滴的滑移速度方程，目的是为了将混合物流动变量分解成两相各自独立的流动变量。两相流中的滑移速度 ^ 是指两相流动速度之差，即超冷水滴速度为¹ ^与空气运动速度之差，

(42) 在求解混合流体的控制方程之后，按照上式可获得滑移速度 Ϋ_Ϋ，之后按照下式求取超冷水滴与混合流体的相对速度 Ϋ

相对速度定义为单相流体速度与混合流体速度之差，所以超冷水滴的速度 ^即可获得 v₂ = v_m2 + v„ (44)

空气-超冷水滴的单流体两相流流动控制方程组包括：连续、动量、能量方程、两相体积分数的扩散方程、湍流模型、空气-超冷水滴的滑移速度方程，模拟结果给出旋翼外部近场的流动信息，即各个计算网格点上（或是网格单元内部）的空气和超冷水滴的密度、压力、速度，以及由上述独立变量推导出的温度、动力粘度、湍流特征等信息。

尾迹

尾迹包含在近场中，所需要的数值模拟假设条件与近场的相同。尾迹是以空气-超冷水滴两相旋涡流动为主要特征的，将描述其运动规律的单流体模型包括组分扩散方程、连续方程，以及本发明给出的动量方程、能量方程联立的积分守恒形式为 i^ + i。(U = i^ ，（⁴⁵) 式中守恒变量 _m、对流项向量 _m、粘性项向量 _m、滑移速度弓 I起的对流项向量、源项 _m，分别为，

0

Pm^C2 (l ~ ^C2

0

其中，混合流体的速度不变量和滑移速度 ^ = uj + vj + w_sk的速度不变量分别定义为

(47)

V,, - U K + VM„ + WM (48)

+

混合流体的转内能 E™和转焓 H™的表达形式分别与公式（40) 和（41) 相同，差别在于将混合流体速度替换其中的空气速度。公式（46) 中， ί代表体积力，即单位质量重力和涡

(49)

此外，本发明中加入了旋转坐标系下的超冷水滴滑移速度方程，见公式（18)，其具体的维速度分量形式是

水膜运动

如前所述，旋转轴为 Z-轴，机翼旋转角速度向: :ώ_ζο =[θ,0,ω_ζο], 在 (ο, ,ζ",77)坐标系下，公式（29) 和 Z-轴旋转角速度分别表示为 ω_ζ ² (x cos 6* + j sin θ)

gcf g_z sin α-ω_ζ ²(χ cos αύηθ-y cos a cos θ) (51) g_z cos a + ί¾?_ζ ²( sin αύηθ-y sin a cos θ)

得上式带入公式（35) 中，并消去二阶小量；；²，可得 β在 (^ς·,//)方向上的三个分

进而，公式（34) 中，代表科里奥利力影响系数张量矩阵 K

求解水膜的平均速度向量时，可以消去科里奥利力影响张量系数矩阵与重力、离心力向: 的点乘运算中的二阶小量 /；²，于是水膜速度公式（36) 变成

其中， Ϊ是单位对角阵；剪切应力向量 ¾_m=r_m + r_mi ， ^和^^分别为 ^: 在 ^,ζ")两个方向上的分量。如果^，^ 是^分别在 ^,ζ·,//)方向上的分量，贝 IJ

其中，向量 Α， Β, έ分别为

A (57)

按照公式（36)，水膜的平均速度为，

按照公式（37) 重新

其中， m_x是超冷水滴的体积分数《₂和超冷水滴速度¹^的函数； m₂和^是水膜温度 7；的函数，其表达式是公知的； m_x是超冷水滴碰撞到水膜所加入的质量流率

m₀是水膜表面蒸发 /汽化流出的质量流率

其中，是气体常数；是水的传热系数；是刘易斯数，此处为 1-、和压力，是温度 Γ„的函数，工程上常以多项式拟合的形式给出，即

(T_w ) = a₀+a, (T_w - 273) + ₂ (T_w - 273)² + ₃ (T_w - 273)³ + ₄ (T_w - 273)⁴， (62) 其中，常数 a。=611.01、 a, = 44.4816 α₂ =1.4188 α₃ - 0.0239 α₄ = 0.0002；式中，相对湿度定义为

Ρ,

Φ. (63)

Ρ_: 按照公式（38) 重新写出水膜结冰的能量方程，并假设其中 =0，则有

其中， Q是超冷水滴碰撞到水膜所加入的能量，由公式（39) 给出;

d是水膜表面蒸发 /冰层升华流出的能量

(65) Q₃ =(L_f -C_piT_w)m, (66) 其中，和£ ^是冰点的蒸发潜热和升华潜热、是冰融化潜热。

在本实施方案中，假设水的定压比热 ^、冰的定压比热 „、水的密度 p_w、水的运动粘度// _w均为常数，水膜结冰过程的质量守恒和能量守恒方程（59)、 (64) 联立，其积分的守恒形式为

aw

·+ GdS = Q, (67) 其中，守恒变

W (68) h_fT_w 对流项向量 0 = f + F 在 ^方向上的分量分别为 ρ„ω_ζ ²(χ COS0 + ysinO) ₊ ^ ₄2 ?_wo_zr_mc cosa

3

F (69)

²A, ¾ cos a

p_w (g_z sin α-ω_ζ ²(χ cos α ηθ-y cos a cos 6¹)) ₄2 ?_wo_zr_mi cos

F

_T

(70) 源项表示为

求解方程组（67) 的约束条件为 h_f >0; mice > 0； h_f -T >0; m_ice-T <0, (72) 可求得/ ^、 T_w、和 _W3，再将 w₃除以冰的密度 _A.，转化为冰层厚度，即 h.=^。 (73) 上述方程（67)求解采用有限体积法，变量存储在计算网格中心，方程组的空间离散二阶 Roe方法、时间离散采用 LU-SGS方法，具体细节均是公知的。图 7给出积分时间为 7分钟计算结束时的直升机旋翼的结冰状态，表示为在不同的与翼展方向垂直的剖面上的结冰形状。图中，

平面 1: 位于距离机翼的翼尖端部 2倍弦长、与翼展方向垂直的剖面；

平面 2: 位于距离机翼的翼尖端部 1倍弦长、与翼展方向垂直的剖面；

平面 3: 位于距离机翼的翼尖端部 0.5倍弦长、与翼展方向垂直的剖面。在不同平面上，机翼前缘部位的结冰形状不同。例如，平面 1在前缘上沿来流方向较厚，结冰厚度大约 0.02倍机翼弦长。平面 2上的结冰形状表现为在机翼上下方向较厚，在前缘上沿来流方向较薄，说明离心力和科里奥利力的共同作用，这也是旋转机翼与固定机翼表面结冰状态的区别。

Claims

权利要求

1. 一种用于模拟直升机旋翼的飞行结冰的数值方法，其特征为包括在旋翼尾迹中、旋转坐标系下在用于模拟空气 -超冷水滴两相流的旋涡运动的单流体模型中的动量方程和能量方程中加入涡量补偿力的算法、在超冷水滴的滑移速度方程中加入离心力和科里奥利力的算法；考虑离心力和科里奥利效应的水膜运动和结冰模型；利用上述算法和模型进行直升机旋翼飞行结冰数值模拟计算的流程。

2. 根据权利要求 1所述的一种用于模拟直升机旋翼的飞行结冰的数值方法，其特征在于所述的在旋翼尾迹中、旋转坐标系下在用于模拟空气 -超冷水滴两相流的旋涡运动的单流体模型中的动量方程和能量方程中的源项位置加入的涡量补偿力等于混合流体的密度 ^乘以单位质量涡量补偿力 i_ffl，其中的定义为

其中， _m代表直角坐标系下空气-超冷水滴两相流涡量；符号 X表示叉乘运算； fi_ffl代表涡量最大梯度变化的方向；代表涡运动补偿粘性系数。 3. 根据权利要求 2所述的涡运动补偿粘性系数，其特征在于，其表现形式为 v_ffl = ^，其中， ^代表补偿涡量的特征半径，二维和三维空间分别为^ = ^Ω^/2 和 ^ ΩΑ ,

2 2 其中， Ω在二维情况下代表计算网格面积、三维情况下代表计算网格体积。

4. 根据权利要求 1所述的在超冷水滴的滑移速度方程中加入离心力和科里奥利力的算法，其特征在于，所述的离心力和科里奥利力以向量的形式。

5. 根据权利要求 1所述的一种用于模拟直升机旋翼的飞行结冰的数值方法，其特征在于所述的考虑离心力和科里奥利效应的水膜运动模型中，在水膜运动的当地坐标系 (//, 下，水膜内高度为 /；的平面 ζ")上的速度向量为

其中， /„和 ^分别是水膜的动力粘和密度、 ^和 _m分别是水膜 /空气-超冷水滴两相混合物接触面上在（ ζ")两个方向上的切应力向量和动力粘度、是水膜高度、 i_ge 是单位质量重力向量和离心力向量 iL_fr的合力在当地坐标系 ζ·)下的投影、是科里奥利力系数张量矩阵。

6. 根据权利要求 5所述的科里奥利力系数张量矩阵，其特征在于，其表现形式为

1 — ^3 k₂

h 1 -k

1 其中的元素 [^,^Λ]是向量在 ( ,ζ",77)方向上的三个分量，向量 S为「 τ -2pjj{h_f -η) ^

k = Ί 3」 = ^ω0Ζ,/，其中， ώ。_ζ; 是直升机旋翼的旋转角速度 ώ。_ζ 在坐标系 (/；， ^下的投影。

7. 根据权利要求 1所述的一种用于模拟直升机旋翼的飞行结冰的数值方法，其特征在于所述的考虑离心力和科里奥利效应的结冰模型中，使用权利要求 5所述的水膜平均速度。

8. 根据权利要求 1 所述的一种用于模拟直升机旋翼的飞行结冰的数值方法，其特征在于所述的利用上述算法和模型进行直升机旋翼飞行结冰数值模拟计算的流程，具体是

(1) 未结冰的旋翼周围生成旋转坐标系和计算网格；

(2) 划分远场、近场和尾迹三个计算区域；

(3) 规定计算的开始时刻 t;

(5) 进行旋翼近场和尾迹空气-超冷水滴的运动的单流体两相流流动模拟，给出两相混合流体的压力、速度、温度、动力粘度系数、剪切力；

( 6 ) 求旋翼壁面水膜运动的平均速度

(7) 求水膜厚度；

(8) 求结冰量和结冰厚度；

(9) 获得新的冰层形状；

(10) 重新划分计算网格；

(11) 回到步骤（2) ，进行下一个时刻的流体两相流单流体流动模拟计算。