基于格子玻尔兹曼方法计算表面结冰水滴撞击特性的方法

摘要

本发明提出一种基于格子玻尔兹曼方法计算表面结冰水滴撞击特性的方法，以介观尺度的、远小于宏观控制体的流体微团为基础，通过求解带有大涡模拟、分裂碰撞模型、差异格子速度格式的多组分格子玻尔兹曼方程，采用空气组分平衡态退化、引入水气分子聚合数的处理方式使模型更符合物理实际，从而求得空气和水滴的宏观物理量，进而求得水滴撞击特性。本发明可以以流体分子微团为研究对象求解空气水气混合流场，捕捉宏观方法下无法模拟的水滴微观运动；在复杂结冰条件下，有利于解决求解流场时由于复杂几何而面临的网格尺度的约束问题。

说明书

技术领域

本发明涉及计算流体力学技术领域，尤其涉及一种基于格子玻尔兹曼方法计算表面结冰水滴撞击特性的方法。

背景技术

飞机在穿过含有过冷水滴的云层时会发生结冰现象。在数值模拟结冰的过程中，水滴撞击特性的计算是预测结冰冰型的基础。水滴对表面的撞击区、撞击量，以及水滴在撞击区内的分布，统称为水滴对表面的撞击特性，如图1所示。其中水滴撞击极限与局部水滴收集系数是撞击特性的重要参数。

格子玻尔兹曼方法是一种介观的计算流体力学方法，其建立了微观分子运动和宏观流体流动的桥梁，基于流体分子在网格格点的碰撞迁移理论求解速度分布函数(particle distribution function,PDF)，从而获得网格格点上的宏观物理量，比如速度，密度，压力等。图2为二维情况下流体分子的演化示意图，可以看出流体分子在8个方向上进行迁移。格子玻尔兹曼的计算尺度是介观的，以图2为例，可以理解为，通过研究流场中每个格子上具有8个方向运动趋势的流体分子微团来研究整个流场的运动，该微团并非单个流体分子，而是具有相同运动趋势的流体分子组成的微团，该微团比分子大但又远小于宏观模型的控制体。利用微观粒子的速度分布函数演化，可以更加直观的描述流体组分间的相互作用。因此，格子玻尔兹曼方法在多相态、多组分的复杂流动现象方面受到了很大的关注。

现有的技术方案中，欧拉法求解水滴撞击特性的方法是空气-水滴两相流法(参考：林贵平，卜雪琴，申晓斌，郁嘉.飞机结冰与防冰技术[M].北京:北京航空航天大学出版社，2016.)，两相流即气相和液相两种相态。欧拉法将水滴看成连续相，再引入水滴容积分数的概念后，求解水滴运动的连续方程和动量方程，其中水滴容积分数的定义是水滴相在控制体所占的体积与总控制体积之比。欧拉法计算的尺度是宏观的，宏观模型可以理解为，通过研究流场中一个个固定的有限控制体从而研究流场运动，其中有限控制体中包含有一定数量的流体分子以及流体分子的无序运动。在结冰应用中，欧拉法在宏观控制体的前提下计算水滴撞击特性，目标是模拟宏观的结冰冰型即只关注冰的主要特征。因此，欧拉法不考虑分子间作用的影响，也不考虑水滴的破碎现象，并且，欧拉法的水滴收集率算法中，水滴流场与空气流场是相对独立的两套计算体系，空气流场在经过迭代后施加作用力到水滴流场上，进而实现水滴流场的运动与迭代。实际情况下，几微米的水滴和空气之间的相互作用以及水滴的二次碰撞和破碎，会影响水滴撞击特性从而影响结冰，并且，空气、水滴运动的过程涉及相同分子和不同分子之间的相互碰撞和相互作用，通常会出现反常扩散等现象。为了捕捉这些特征，需要从流体分子运动的尺度上研究水滴撞击特性问题。从分子运动的微观尺度计算来说由于计算机的限制，其模拟的分子个数十分有限，并不适合模拟大量分子的运动，介观尺度方法中针对多组分气体运动问题的离散速度法，由于方法的复杂，其只能运用于低维、简单的运动，对于结冰的复杂流动并不适用。

发明内容

本发明针对在流体分子运动的介观尺度上研究表面结冰水滴撞击问题，提供了一种新的基于格子玻尔兹曼方法计算表面结冰水滴撞击特性的方法。该方法可以计算空气-水滴混合绕流流场中空气和水滴的运动状态、计算表面水滴撞击特性。

本发明的技术方案为：

一种基于格子玻尔兹曼方法计算表面结冰水滴撞击特性的方法，包括以下步骤：

步骤1：获取流场网格数据，所述流场为需要计算表面结冰水滴撞击特性的物体周围流场；

步骤2：设定流场初始参数，包括空气的初始速度u

利用设定的流场初始参数计算空气组分平衡态分布初始函数f

利用空气组分平衡态分布函数和水气组分平衡态分布函数分别对应计算空气组分和水气组分的宏观量：空气密度ρ

步骤3：根据空气组分和水气组分的宏观量以及混合流体的宏观量，计算得到空气组分平衡态速度分布函数

步骤4：根据步骤3得到的空气组分平衡态速度分布函数

空气组分碰撞项

水气组分碰撞项

步骤5：利用步骤4得到的空气组分碰撞项

计算碰撞后的速度分布函数，包括空气速度分布函数

其中

再由非网格格点上的水气速度分布函数

步骤6：计算流场边界处的速度分布函数；

步骤7：由步骤5求得的空气组分和水气组分的速度分布函数，计算各组分的密度和宏观速度及压力，以及混合流体的密度、速度和压力；

步骤8：判断是否达到收敛条件，如果达到收敛条件，则输出流场的计算结果，否则返回步骤3；

步骤9：利用步骤8得到的流场计算结果，利用公式

计算局部水滴收集系数β，其中ρ

进一步的，步骤1中获取的流场网格数据包括为实现流场计算而所需的网格特征参数、边界点LBM网格数据文件、流场点LBM网格数据文件、网格单元连接关系数据文件以及网格单元信息文件。

进一步的，步骤2中，空气组分平衡态分布初始函数

其中c

水气组分平衡态速度分布初始函数

其中

进一步的，步骤2中，计算模型速度分布选择D2Q9。

进一步的，步骤3中，

空气组分的平衡态速度分布函数

水气组分平衡态速度分布函数

进一步的，步骤5中，根据公式

由非网格格点上的水气速度分布函数f

进一步的，步骤6中，边界处分布函数的处理方式为：远场边界采用非平衡态外推格式，分层边界采用虚拟分层边界高精度处理格式，物面边界采用曲面边界处理方法。

有益效果

本发明可以解决在流体分子运动的介观尺度上研究表面结冰水滴撞击特性问题；可以以流体分子微团为研究对象求解空气水气混合流场，捕捉宏观方法下无法模拟的水滴微观运动；在复杂结冰条件下，有利于解决求解流场时由于复杂几何而面临的网格尺度的约束问题。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1展示了水滴撞击示意图；

图2展示了二维情况下在一个格子上的流体分子的演化过程；

图3展示了差异格子速度格式下分布函数插值；

图4描述了算法的流程以及各个模型之间的关系；

图5描述了实施例流场计算数据流程；

图6展示了圆柱算例网格示意图；(a)计算域示意图，(b)网格局部示意图；

图7展示了翼型算例网格示意图；(a)计算域示意图，(b)网格局部示意图；

图8展示了针对圆柱表面空气和水滴的运动轨迹；(a)空气(b)水滴；

图9展示了圆柱表面局部水滴收集率的计算结果；

图10展示了针对翼型表面空气和水滴的运动轨迹；(a)空气(b)水滴；

图11展示了翼型表面水滴收集系数的计算结果。

具体实施方式

为了开展在流体分子运动的介观尺度上研究表面结冰水滴撞击特性问题，本发明结合格子玻尔兹曼的算法特点以及结冰应用场景下流场的真实物理特性，提出了流场计算模型假设、并且构建流场计算模型、建立局部水滴收集系数的计算模型，以及相应的提出了针对圆柱表面、翼型表面的基于格子玻尔兹曼方法的水滴撞击特性数值模拟方法。

首先给出流场计算模型假设：

在结冰场景的应用下，空气水滴混合流场是一种云雾场，其中水滴的直径在微米级别，水滴均匀的遍布在整个流场。此时，在分子尺度上追踪气相和液相的界面几乎是不可能的，大量的微小的水滴弥漫在空间中，两者更类似于两种相互混溶的流体。此时就需要引入多组分模型来重现这种现象。同时，格子玻尔兹曼方法具有瞬时特性，其可以将计算步时间步长控制在在极短的物理时间内，以初始马赫数Ma＝0.1732，网格尺寸dx＝0.001，特征长度为1米的计算条件为例，每当LBM计算模型演化一千万步，对应真实物理世界只经过了一秒。在这个极其短暂的时间长度内，如果水滴与来流速度一致，则水滴仅仅运动了0.06微米，该距离为水滴直径的0.6％，也就是说水滴并没有完全离开其原本的位置，而是以其相界面的形变完成了这一过程，从宏观角度上可以近似认为该水滴在原地几乎没有发生位移。但如果此时以流场全局的宏观视角观察这一微观位移过程，该位移过程又过于细小且难以捕捉，水滴团的位移是微观概念，不能在全局流场的宏观尺度上进行准确描述。在以上基础上，为了能够正确描述云雾场中液体的运动状态，基于格子玻尔兹曼方法，本发明提出了，水气和空气单相多组分混合流场的计算模型，通过水气流场近似等效水滴流场来解决以上问题。以下为流场计算模型假设：

(1)气相水趋于发生相变成为液相水；

现实中，流场中的气相水由于自身的温度、密度或是压力影响，处于一种趋向于发生相变变成液相水的状态中，该趋势与化学反应中的平衡状态类似，水在一定的温度、压力、压强下可以互成平衡，形成一种冰-水-汽的三相平衡体系(参考：刘飞.水和甲醇在二维受限下相变的分子动力学模拟研究[D].中国科学技术大学,2011.)。这也是流场中的气相水可以转变成为液相水的根本依据。

(2)流场中每一节点上都有足够的等效水滴分布；

高密度的水汽是难以被观察到的，该状态会短暂存在，是只能在较短的时间内维持的一个特殊的过度状态，但是该状态客观存在在由气相水向液相水转变的物理过程中，在水蒸气冷凝时，蒸汽(气态)变化为水(液态)，同时释放出潜热(参考：鲍鹤鸣,高淑宁,程思源,关欣.水的相变过程中辐射机理研究[J].能源研究与信息,2021,37(01):40-45.DOI:10.13259/j.cnki.eri.2021.01.007.)。因而将此瞬间的流场作为流场数值计算的理论依据，充分利用格子玻尔兹曼方法的时间步长小的特点，建立数值模型进行仿真。

(3)在基于格子玻尔兹曼方法的流程计算过程中等效水滴不发生相变；

在流场计算的过程中，在从一节点向另一节点迁移的过程中，任何等效水滴在当前瞬间都认为是气相，不会发生相变，也就是说被等效的等效水滴不会转化为液相水或者冰。该点假设与欧拉法(参考Lee T,Lin C L.An Eulerian description of thestreaming process in the lattice Boltzmann equation[J].Journal ofComputational Physics,2003,185(2):445-471.)一致，在空气中运动的水滴不会发生升华或者冻结现象进而造成质量、能量损失或交换。

(4)相变过程足够短暂；

为了保证流场的整体流动特性不发生巨大变化，故而假设所有的相变过程都可以在小于一个流场计算时间步的时间长度内完成，这样一来既可以保证水滴在接触物面的瞬间完成相变过程形成结冰所需要的表面水膜，也可以保证流场不会随着时间推进而逐渐在物理性质上发生偏移。

综合假设(3)与假设(4)，在计算过程中的体现就是水的初始相态在流场中始终等于水的平衡相态，并且在接触到物面之前水的相态始终不发生变化，但又维持一个趋向于相态变化的状态。

(5)空气流场对水滴流场的影响是单向的；

由于在结冰过程中水滴直径一般在1-100微米的范围之内(参考：易贤，朱国林.积冰面水滴撞击特性的数值计算.空气动力学研究文集，第15卷，2006)，近似认为它对空气流场没有任何影响，这一思想在欧拉法(参考Durst F,Miloievic D,Schnung B.Eulerianand Lagrangian predictions of particulate two-phase flows:a numerical study[J].Applied Mathematical Modelling,1984,8(2):101-115.)的构建过程中均有体现。

(6)水滴流场与空气流场直接没有热交换；

由于在流动过程中，空气与水滴之间的能量交换相对较少，对于流场整体流动的影响微弱。如果忽略该过程对于数值计算的精度影响也较低，故认为水滴与空气之间不存在热交换(参考Bragg M B.Rime ice accretion and its effect on airfoilperformance/[J].Dissertation Abstracts International,Volume:42-07,Section:B,page:2913.1982.)。

基于以上假设，可以得知水气流场发生相变和格子玻尔兹曼数值模拟过程都足够短暂，水气流场近似等效云雾场的方法，可以通过单相多组分流动中组分间的相互作用来模拟真实的物理运动。

其次构建流场计算模型：

针对格子玻尔兹曼框架下的表面结冰水滴撞击特性计算问题，在研究同相态多组分的流体流动时，对于流场计算模型的搭建，主要考虑如下4个因素：

(1)在大气中，水雾与干燥空气的粘性、质量是不同的，模型必须可以模拟粘度不同的、组分流动可以影响混合流动的流动现象。

(2)飞行过程中的飞行器，往往需要面对较高雷诺数，在这种状态下，必须要选用合理的湍流模型。

(3)结冰应用下云雾场的液态水含量的量级约为1g/m

(4)在结冰场景下的空气水滴混合流场，从宏观上看，水滴的质量和惯性都比空气的要大。而在分子运动的角度，空气的平均分子质量为29，水分子的分子质量为18，因此在求解多组分流动时也要考虑实际物理现象，保证水气的惯性比空气组分的惯性大。

结合上述1、2点因素，对适用于求解结冰应用下的空气-水气混合流场，本发明对动力学模型的搭建选择带有分裂碰撞模型、引入大涡模拟和差异格子速度格式(DifferentLattice Speed Scheme，DLS)的多组分格子玻尔兹曼模型，同时针对第3点因素，本发明提出了对基组分即空气组分的平衡态退化的处理方法。针对第4点因素，本发明引入水气分子聚合数，将数个水气分子组合成流体微团，以达到其惯性比空气大的目的。

对于多组分流动而言，某一组分i的格子玻尔兹曼方程(参考Joshi A S,Peracchio A A,Grew K N,et al.Lattice Boltzmann method for continuum,multi-component mass diffusion in complex 2D geometries[J].Journal of Physics D:Applied Physics,2007,40(9):2961)如下：

公式中

其中c

等式(1)右边的碰撞项展开如下：

其中右边一项表示自体碰撞项，其用BGK近似，该项的计算公式为：

公式中τ是自碰撞松弛时间，

其中以图2模型为例，权重w

根据stefan-Maxwell模型，总密度ρ和宏观速度u满足如下公式：

为了提高计算雷诺数的需求，在此引入大涡模拟。加入Smagorinsky涡粘模型后，涡粘性系数ν

其中S

湍流弛豫时间τ

τ＝τ

等式(4)中的互体碰撞项由下式计算：

对于本文的气液两组分而言，互体碰撞弛豫时间τ

其中M

p＝ρ

考虑到水滴对空气的影响可以忽略不计，本发明对空气组分进行了平衡态退化处理。平衡退化的方法是：空气组分平衡方程中的互体碰撞项为零，其方程退化为单组分流动平衡方程，平衡态函数中混合流体的速度u等于空气组分的速度。

因此，空气组分的平衡态退化的格子玻尔兹曼方程具体表达式为：

平衡态分布函数f

然后建立局部水滴收集系数计算模型：

在格子玻尔兹曼的框架下以水气分子微团为研究对象，本发明中表征水滴撞击特性的局部水滴收集系数β的表达式为：

上式即微元表面的实际水收集量与该微元表面上最大可能的水收集量之比，因此它是表征微元表面的水收集能力的一个参数。其中β为局部水滴收集系数、ρ

综合以上，算法流程以及各个模型之间的关系如图4所示，图中第一步流场初始化可以得出初始的宏观量及速度分布函数；第二步根据初始的宏观量以及速度分布函数计算平衡态分布函数；第三步计算碰撞项，碰撞是速度分布函数变化的原因；第四步分布函数的迁移和插值，此过程中不改变分布函数的值，其只在格子节点上迁移，并对偏离格子节点的组分，通过插值求解节点上的分布函数。

下面以计算圆柱表面和翼型表面水滴撞击特性为例，进一步说明本发明：

实施例1：圆柱表面水滴撞击特性数值模拟：

计算表面水滴撞击特性，首先要获得流场的运动状态。欧拉法中，以宏观尺度的控制体为基础，先求解空气的运动状态，再作为已知条件施加到水滴上，求解水滴的连续方程和动量方程，最后获得水滴容积分布和水滴速度分布，进而求得水滴撞击特性。本发明中，以介观尺度的、远小于宏观控制体的流体微团为基础，通过求解带有大涡模拟、分裂碰撞模型、差异格子速度格式的多组分格子玻尔兹曼方程，采用空气组分平衡态退化、引入水气分子聚合数的处理方式使模型更符合物理实际，从而求得空气和水滴的宏观物理量，进而求得水滴撞击特性。

表1为圆柱算例的计算条件

表1圆柱算例计算条件

该算例参考(陈金瓶,二维圆柱结冰的冰风洞试验研究及水滴撞击特性计算[D].南京航空航天大学,2013)。图5为计算过程的数据流程图，如图所示，图5中虚线代表水气组分数据流动，点线代表空气组分数据流动。在整个流场计算过程中需要外部输入的数据有网格信息与流动的初始分布参数。图6为圆柱算例网格示意图，如图所示，网格采用非均匀的自适应笛卡尔网格，其中分层边界为非均匀网格交界处、物面边界为网格与物面的相交处、远场边界为外流场边界处。网格中流体(空气与水气)的初始信息是同时生成的。结合图4和图5具体的计算步骤如下：

步骤(1)，输入圆柱模型的流场网格数据。网格数据包括为实现流场计算而所需的网格特征参数、边界点LBM网格数据文件、流场点LBM网格数据文件、网格单元连接关系数据文件、全部网格单元信息文件。以上网格信息能够确定网格的节点数、网格的边界信息和几何信息、网格点的类型(边界点还是流场点)。

步骤(2)，流场初始化。流场初始化的参数通过人为给定，其中包括：空气的初始速度u

求解平衡态退化形式的空气组分各速度方向平衡态速度分布函数f

其中初始化时：u

w

水气组分平衡态速度分布函数

其中

通过空气和水气的平衡态分布函数

混合流体的宏观量：密度ρ、速度u、压力p，如下公式求得：

p＝ρ

初始化过程后得到了为后续迭代结算所需的宏观量和平衡态分布函数。

步骤(3)，演化与迭代计算，计算平衡态分布函数。平衡态分布函数计算步骤和公式同步骤(2)，其中所需的宏观量：空气密度ρ

空气组分各速度方向平衡态速度分布函数

水气组分平衡态速度分布函数

步骤(4)，计算碰撞项。由第(3)步得到的平衡态分布函数

对于空气组分的碰撞项

其中：

对于水气组分碰撞项的计算：

其中

步骤(5)，计算空气的速度分布函数

其中

求解网格格点上的水气速度的分布函数

其中

步骤(6)，计算边界处的分布函数。第(5)步的计算得到了所有流场点的速度分布函数，接下来计算边界处的速度分布函数。边界处分布函数的处理：远场边界采用非平衡态外推格式(参考Guo Z L,Zheng C G,Shi B C.Non-equilibrium extrapolation methodfor velocity and boundary conditions in the lattice Boltzmann method.ChinesePhysics,2002,11(4):0366-0374.)；分层边界采用虚拟分层边界高精度处理格式(参考西北工业大学.一种针对树网格格子玻尔兹曼方法中虚拟分层边界的高精度处法:CN202210219613.9[P].2022-06-07.)；物面边界采用曲面边界处理方法(参考西北工业大学.针对格子玻尔兹曼方法数值模拟中曲面边界的处理方法:CN202210219614.3[P].2022-06-07.)。

步骤(7)，计算宏观量。由步骤(5)求得的空气和水气组分的速度分布函数

混合流体的宏观量如下公式求得：

其中ρ是混合流动的密度、u是混合流动的速度、p是混合流动的压力、u

步骤(8)，达到收敛标准后，输出流场的计算数据及结果图像处理。其中包括密度、速度、压力的计算结果。

步骤(9)，计算局部水滴收集系数。局部水滴收集系数在公式上可表示为：

其中β为局部水滴收集系数、ρ

图8给出了圆柱表面水滴撞击特性数值模拟算例的空气(a)和水滴(b)的运动轨迹,从图中可以发现，空气在接近圆柱表面时流线发生剧烈弯曲，气流贴着上下表面流过，表现出明显的绕流特性，在尾缘部分，气流又重新恢复；同时，可以明显的看出水滴的收集特性，由于水滴质量和惯性较大，不容易改变运动状态，因此水滴无法绕过圆柱前缘，从而撞击在壁面上，由于水滴的收集特性，因此在尾缘部分出现了水滴遮蔽区，即此部分没有水滴的存在。

图9给出了圆柱表面局部水滴收集系数的对比结果，可以看出，局部水滴收集系数在0到1的范围内，在前缘驻点处水滴收集系数最大，并且收集系数沿上下表面逐渐减小，直到水滴撞击极限处收集系数为零。本发明计算结果与欧拉法计算较为吻合。

实施例2：翼型表面水滴撞击特性数值模拟：

表2为翼型算例的计算条件

表2翼型算例计算条件

该算例参考(AIAA.Results of an icing test on a NACA0012airfoil in theNASA Lewis Icing Research Tunnel-30th Aerospace Sciences Meeting and Exhibit(AIAA).1992)，在进行翼型表面水滴撞击特性的数值模拟时，所用到的计算模型和计算流程与实施例1相同，其中图7为翼型算例网格示意图。图10给出了翼型表面水滴撞击特性数值模拟算例的空气(a)和水滴(b)的运动轨迹,同样可以发现，水滴在翼型前缘的撞击和收集。图11为翼型算例的局部水滴收集率对比图，可以看出，计算结果与试验的计算结果整体比对较好。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。