一种基于通量限制器的水滴流场及水滴收集率计算方法

摘要

本发明适用于计算流体动力学领域，提供了一种基于通量限制器的水滴流场及水滴收集率计算方法，在计算流出控制体的对流净通量的过程中引入通量限制器，通过合理构造通量限制器的使用模型，在提高水滴流场计算精度的同时，也能解决其求解过程中发散的问题，非常适用于飞机结冰计算中水滴流场的高阶精度计算。

说明书

技术领域

本发明涉及计算流体动力学领域，尤其是涉及一种基于通量限制器的水滴流场及水滴收集率计算方法。

背景技术

飞机穿过含有过冷水滴的云层时，水滴撞击到飞机表面可能凝结成冰，结冰会改变飞机升力部件的外形并影响其气动性能，导致飞机的操纵性和稳定性恶化；其他部位结冰也可能导致传感器失灵、操纵面卡死、发动机失火等不良影响，严重时直接导致空难事故。

数值模拟是研究飞机结冰的重要手段，飞机结冰是水滴在空气绕流作用下，撞击到飞机表面发生凝结的现象。其中，水滴收集率计算是整个结冰计算的基础，是其他模块的输入依据，同时，水滴收集率也是防冰系统设计的重要输入依据。

水滴收集率的计算方法有拉格朗日法和欧拉法两种。对于欧拉法，一阶格式的精度不够，而单纯的高阶格式又容易引起数值振荡和密度脉冲。

目前用欧拉法高阶格式求解有两种解决方法，一是在水滴流场控制方程中加入扩散项，如Tong等分析了欧拉法数值不稳定的原因，在原始的水滴控制方程基础上添加了扩散项用来消除迭代过程中的奇异性，招启军等建立了遮蔽区扩散模型，解决尾流等区域的密度脉冲现象引起的稳定性和收敛性问题。该方法的稳定性较好，但是需要修改水滴的控制方程，增加了编程的复杂度。

二是简单的使用带有通量限制器的高阶格式，如张大林等使用欧拉法MUSCL离散格式计算了二维NACA0012翼型的水滴收集率，易贤等使用欧拉法发展了适合于三维、复杂构型的水滴收集率计算程序，其对流项的离散是采用迎风插值和线性插值相结合的方法。使用带单一简单的通量限制器的格式可能会无法解决复杂构型的情况。

综上所述，现有技术对水滴流场的模拟计算存在着计算过程中发散导致的计算精度低，或者计算量大，计算复杂的技术问题。

发明内容

为了解决以上技术问题，本发明提供一种基于通量限制器的水滴流场及水滴收集率计算方法。

现有技术中，使用高阶格式求解水滴流场时，一般有两个原因导致迭代求解过程的发散：

其一，水滴流场中，虽然不存在空气流场中的激波现象，但是由于物面遮挡效应所导致的无水区边缘的液态水含量的梯度很大，使用高阶精度格式计算时也会出现数值振荡。

其二，由于欧拉法中水滴速度的单值性引发水滴流动交汇处出现奇异值。在物体后部一般会有水滴交汇的现象出现，对于拉格朗日法，每个水滴是独立运动的，其轨迹可以发生交叉和重叠，此处液态水含量的值为两股交汇水流的液态水含量之和；但对于欧拉法，由于水滴被视为是连续相，流线不能出现交叉或重叠，两股水滴流动分别从上至下和从下至上汇聚成一股向同一方向流动，导致了此处的液态水含量异常增大，若离散格式缺乏耗散项及时将此处的异常流动抹平，则极其容易导致计算的发散。

因此需要引入限制准则，在对流项的离散格式中添加通量限制器来使其达到总变差减小的条件，以此来保证数值解的单调性和有界性。

通量限制器能够在物理场的梯度不同处，组合使用不同的离散格式(一阶迎风格式、二阶迎风格式、中心差分格式等)。因此可以用来抑制水滴流场间断处的波动和抹平奇异值。

本发明在不添加扩散项的前提下，通过在二阶精度格式中引入系列通量限制器来解决收敛性的问题，以发展飞机结冰计算中水滴流场的高阶精度计算格式。

本发明提供一种基于通量限制器的水滴流场计算方法，其特征在于，包括以下步骤：

S10.确定工况参数、网格类型、边界信息及空气流场信息；

S20.给定初始值，根据水滴流场控制方程及边界条件计算流出控制体的对流净通量和控制体的源项通量；

S21.计算流出控制体的对流净通量

S211.引入通量限制器进行插值计算控制体表面的守恒变量Q

i,j，k分别为控制体在x,y,z方向的编号，ξ，η，ζ分别为曲线坐标系下的三个方向；

ψ(r)为通量限制器函数，通量限制器函数满足以下限制准则：

S212.根据步骤S211计算的Q

F

令

V

将F

Q

将F

S22.计算流出控制体的源项通量

其中，

S30.进行时间推进

根据以下公式计算下一时间步的Q值：

其中，

S40.将步骤S30计算得到的Q值代入到步骤S20-S30；

当迭代步数小于设定值，继续进行迭代计算；

当迭代步数等于设定值，结束计算。

本发明还提供一种基于通量限制器的水滴收集率计算方法，包括以下步骤：

S100.采用如前所述的一种基于通量限制器的水滴流场计算方法进行流场计算，得到Q值；

S200.根据Q值的表达式计算得到水滴速度向量u；

S300.根据以下公式计算水滴收集率β：

其中，α

相对于现有技术，本发明至少具有以下有益效果：

本发明使用通量限制器构造的TVD格式(The total-variation diminishing，TVD，总变差减小)来离散水滴流场的对流项，在提高计算精度的同时，也能解决其求解过程中发散的问题。因为一阶格式的精度低，而单纯的使用二阶格式又会存在间断处的波动或者水流交汇处奇异值的问题，而通量限制器能够根据物理场的梯度来自动调节来抑制数值波动和奇异值。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对本发明实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面所描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例的一种基于通量限制器的水滴流场计算方法的流程图；

图2是本发明实施例的一阶迎风格式和五种带不同限制器的TVD格式计算得到的水滴收集率与实验数据的对比。

具体实施方式

以下的说明提供了许多不同的实施例、或是例子，用来实施本发明的不同特征。以下特定例子所描述的元件和排列方式，仅用来精简的表达本发明，其仅作为例子，而并非用以限制本发明。

一种基于通量限制器的水滴流场计算方法，如图1所示，包括以下步骤：

S10.确定工况参数、网格类型、边界信息及空气流场信息；

本领域技术人员知晓，本申请的水滴流场计算方法基于有限元分析来计算，按照本领域常规的方法进行网格划分和边界信息输入即可，具体网格如何划分，边界如何设定不作为对发明的限制；

S20.给定初始值，根据水滴流场控制方程计算对流通量、空气相对水滴相的外力相；

水滴流场控制方程的积分形式：

其中，u

τ

f为阻力函数：

由于常见的飞机飞行时的雷诺数超出了斯托克斯公式计算水滴粘性阻力的理论允许范围，因此采用修正的阻力函数：

Re

μ

将水滴流场控制方程的积分形式半离散化：

i,j，k分别为控制体在x,y,z方向的编号，Vol

S21.计算流出控制体的对流净通量

S211.引入通量限制器进行插值计算控制体表面的守恒变量

在非均匀结构网格中水滴流场控制体表面左右两侧的守恒变量可以写为：

公式中，ψ(r)为通量限制器函数，通量限制器函数需要满足Sweby’s TVD限制准则：

Sweby’s TVD限制准则不依赖于CFL数，对于求解半离散对流方程的恒定状态解来说，不依赖CFL数的限制准则更加适用。而在空气流场恒定的情况下，水滴流场就是定常的，其求解过程是使用“伪时间步”推进的。“k-格式”系列所包括的区域与Sweby’s TVD限制区域重合的部分，既可以保证总残差减小又能保证空间二阶精度。

不同的通量限制器函数决定了对流项离散格式的精度和收敛性，具体地，作为优选，可选用以下任意通量限制器函数：

Minmod：ψ(r)＝max(0,min(r,1))

Minmod是一种最简单的限制器，使用时，实际上是只用迎风梯度或中心梯度来修正一阶迎风量。当在极值处时r≤0，离散格式直接降为一阶迎风格式，而远离极值处时，离散格式会在中心差分和二阶迎风格式中切换。

MUSCL：

MUSCL限制器在光滑区域具有Fromm格式的精度和收敛性；

Superbee：ψ(r)＝max(0,min(2r,1),min(r,2))

Superbee限制器在y-r图中是二阶Sweby’s TVD限制区域的上界，此格式是被设计用来详细刻画流场间断的，但是由于在高曲率区域大量使用一阶顺风格式导致了光滑区域的过度压缩。和Minmod限制器一样，Superbee限制器在光滑区域也是使用中心差分与二阶迎风格式的组合，不过判断条件与其相反。

Harmonic：

Harmonic格式实质上是用中心梯度和迎风梯度的调和平均值来修正一阶迎风项。其也是基于Fromm格式的斜率限制器，但只在r＝1处具有和Fromm格式相同的斜率。

Van Albada：

Van Albada限制器在整个定义域上都是光滑、连续且可微的，这为离散格式带来了良好的收敛性。与其他限制器不同的是，此限制器在r≤0时并不为0，即完全不使用一阶迎风格式。

S212.根据步骤S211计算的Q

令F

Q

由于水滴流场控制方程的Jacobian系数矩阵有4重特征值V

V

令：

也就是说，先计算V

值得说明的是，式中

S30.计算流出控制体的源项通量

源项通量按下式计算：

S40.进行时间推进

对于水滴流场的半离散式，令：

略去下标i,j,k，则可以写为：

根据步骤S20和S30计算的

当迭代步数超过设定值，则结束计算，迭代步数小于设定值，继续进行迭代计算，并输出每步计算的Q值，作为优选，可以每间隔一定的步数输出Q值。

作为优选，本发明采用四步Runge-Kutta发来进行时间推进，即根据第n步得到的Q值来计算第n+1步的Q值的具体计算过程：

Q

Q

Q

Q

Q

Q

其中α

进一步地，可以根据计算得到的水滴流场信息计算水滴收集率：

经过前述计算得到了Q值，根据Q的矩阵式可以得到水滴速度向量u，

根据以下公式计算水滴收集率β：

其中，α

图2为一阶迎风格式和五种带不同限制器的TVD格式计算得到的水滴收集率与实验数据的对比，可以看出，相比于一阶迎风格式，具有限制器的TVD格式均能更好的计算水滴收集率，计算精度都有较好的提升，特别是在驻点附近，与实验值的耦合度较好。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。