飞行器积冰冰形模拟方法、装置、计算机设备及存储介质与流程

1.本发明属于计算流体力学领域，具体涉及一种飞行器积冰冰形模拟方法、装置、计算机设备及存储介质。


背景技术：

2.由于积冰导致的飞行事故多发，带来了人身安全和社会财产方面的巨大损失。在飞行器飞行过程中，机翼前缘、发动机进气口、机头挡风玻璃以及各种传感器检测仪器的迎风面都是易于形成冰层的区域。飞行器本身是由机翼表面上下压力差产生的升力来实现飞行的。然而由于积冰的存在使得飞行器表面气动外形发生了改变，影响了飞行器的气动特性。研究表明，干净构型气动特性随雷诺数变化明显，雷诺数越大，其最大升力系数、失速越大；但当飞行器积冰后其气动特性几乎不受雷诺数效应的影响。
3.此外，发动机进气口处结冰会减小发动机进气量，导致发动机功率不足，甚至造成空中停车的严重后果。而传感器附近积冰会影响采集数据数值，挡风玻璃处积冰会影响飞行员视野对安全驾驶造成隐患。
4.因此针对积冰冰形的研究具有重要意义，目前针对飞行器积冰的实验成本高昂且难度大，因此数值模拟成为研究积冰问题的重要方法
1.。针对飞行器积冰问题国外学者很早之前就展开了相关研究，已推出多款成熟商业软件产品如美国国家航空航天局开发的lewice,加拿大开发的fensap
‑
ice,意大利开发的cira和ciramil，英国开发的dra和法国开发的onera等，其中最为著名的是美国航空航天局开发的lewice软件。lewice软件以及法国航空研究中心的计算软件被美国联邦航空局和欧洲联合航空局所采用，作为确认飞机是否满足结冰适航的工具。而国内方面的研究起步较晚，空气动力研究与发展中心易贤等、南京航空航天大学孙志国等、西北工业大学周志宏等、北京航空航天大学杨胜华等和常士楠等均采用数值模拟的手段针对飞行器积冰问题进行了相关研究。而且已有研究成果主要为二维问题，针对三维问题的研究正在逐步开展。已开发软件如nuaa
‑
ice3d，数量较少与国外还存在一定差距。
5.运用工程估算的方法主要是使用一些经验总结出来的公式和数据针对飞行器积冰结果进行估算，这种方法比较粗糙，而且只能粗略估算出积冰量，无法模拟出积冰的增长过程和冰形样貌。
6.利用实验方法针对飞行器积冰进行研究主要分为实际环境的积冰实验和冰风洞试验两种。飞行试验对飞行员要求高且存在一定的危险性，成本高昂，而冰风洞方法难点在于完全模拟自然条件下的飞行环境。
7.在数值计算方面。目前国内部分研究人员针对飞行器表面积冰问题采用稳态方法进行模拟，然而飞行器表面积冰是一个动态发展的过程，冰形会随着积冰时间不断发生变化，因而稳态方法与实际积冰情况存在差距。另外，飞行器表面积冰冰形复杂，不规则，且不同冰形样貌对于飞行器表面气动特性影响很大。而部分现有研究中，将冰形归纳为诸如钝头体、双角体和尖头体等类型以对问题进行简化，简化后冰形类型的组合无法涵盖实际复
杂冰形样貌，对计算结果准确度造成影响。与此同时，还存在利用nurbs样条曲线等拟合积冰冰形样貌的方法，利用样条曲线只能光滑近似整个区域的冰层形貌，而无法实现更细致的局部冰形样貌的精细模拟，例如楔形冰表面的毛刺状积冰等。最后，上述采用方法多借助商业软件如fluent等实现，学习成本高，虽然功能完备但是体量庞大，相较于更具有针对性的自编程方法效率不高。


技术实现要素：

8.针对现有技术存在的问题，本发明提供一种飞行器积冰冰形模拟方法、装置、计算机设备及存储介质。
9.为实现上述技术目的，本发明采用的技术方案如下：
10.飞行器积冰冰形模拟方法，包括：
11.输入飞行器翼型，确定翼型上的积冰研究区域，采用正交网格进行网格划分；
12.根据当前积冰研究区域边界区分固体区域和流体区域；
13.通过求解流动方程组，获取当前积冰研究区域中流体区域的速度场分布；
14.输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程；
15.对水滴收集过程中得到的已收集水滴进行积冰生长，积冰沿着当前积冰研究区域边界法向增长，计算积冰研究区域不同位置的积冰生长高度，对积冰研究区域积冰后的边界进行更新，得到积冰生长后新的积冰研究区域边界。
16.进一步地，本发明可以根据需要设定飞行器上的积冰研究区域。如设定飞行器上的积冰研究区域为机翼，输入飞行器翼型以机翼轮廓线作为初始积冰研究区域边界，在机翼轮廓线以内的区域为固体区域，机翼轮廓线以外的区域为流体区域。
17.进一步地，本发明中，对于积冰研究区域的流体区域内的任一网格单元，构建流动方程组，如下：
[0018][0019]
其中u、v为网格单元其x、y方向的速度，r
e
为雷诺数，p为网格单元上的压力，q为添加的质量源项，f
x
、f
y
均为添加的体积力源项。
[0020]
进一步地，本发明对于第i行、第j列的网格单元v
i,j
，网格单元v
i,j
其中心点p
i,j
位于当前积冰研究区域边界内部，其x、y方向的速度分别为u
g
,v
g
，体积力源项的表达式如下：
[0021]
[0022][0023]
其中，δt为时间步，为已知量；
[0024]
设当前积冰研究区域边界与网格单元v
i,j
的两个顶点相交，(u
x,in
,v
y,out
)为当前积冰研究区域边界外部与网格单元v
i,j
边界垂直的速度分量；(u
x,out
,v
y,in
)为当前积冰研究区域边界内部与网格单元v
i,j
边界垂直的速度分量，q
i,j
为p
i,j
点处的质量源项，由质量守恒得：
[0025][0026]
其中，
[0027]
结合连续方程得到：
[0028][0029]
由此得到：
[0030][0031]
差分格式为：
[0032][0033]
其中
ɑ
为网格单元边界被当前积冰研究区域边界所占据的比例，当网格单元边界全部在当前积冰研究区域边界内部时；
ɑ
＝1；反之，当网格单元边界全部位于当前积冰研究区域边界外部时；
ɑ
＝0；
[0034]
将体积力源项、质量源项代入到流动方程组中，联立求解得到积冰研究区域的流体区域内的任一网格单元的u、v，得到积冰研究区域的流体区域内的任一网格单元的气流速度进而得到流体区域的速度场分布。
[0035]
进一步地，本发明对水滴进行受力分析，仅考虑重力、浮力以及空气阻力的影响，由牛顿第二定律，采用拉格朗日方法得到水滴运动轨迹方程，进而获取水滴移动轨迹。
[0036]
进一步地，本发明中的水滴运动轨迹方程为：
[0037][0038]
上式中，m
d
为水滴质量，为水滴位移，ρ
a
、ρ
d
分别代表空气密度和水滴密度，v
d
为水滴体积，为当地重力加速度，a
d
代表水滴迎风面积，c
d
为阻力系数，为阻力系数，分别代表气流速度和水滴速度。
[0039]
进一步地，本发明中对于当前积冰研究区域边界上的任一网格单元，其将凝结成
冰的水滴质量等于流入该网格单元的水滴质量减去流出该网格单元的水滴质量；
[0040]
则第i个时间步内，当前积冰研究区域边界上的任一网格单元内的积冰生长高度为：
[0041][0042]
其中ρ
ice
表示冰的密度，a
s
代表网格单元面积。
[0043]
进一步地，本发明提供一种飞行器积冰冰形模拟装置，包括：
[0044]
第一模块，用于输入飞行器翼型，确定翼形上的积冰研究区域，采用正交网格进行网格划分；
[0045]
第二模块，用于根据当前积冰研究区域边界区分固体区域和流体区域；
[0046]
第三模块，用于通过求解流动方程组，获取当前积冰研究区域边界中流体区域的速度场分布；
[0047]
第四模块，用于输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程；
[0048]
第五模块，用于对水滴收集过程中得到的已收集水滴进行积冰生长，积冰沿着当前积冰研究区域边界法向增长，计算积冰研究区域不同位置的积冰生长高度，对积冰研究区域积冰后的边界进行更新，得到积冰生长后新的积冰研究区域边界。
[0049]
本发明提供一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
[0050]
输入飞行器翼型，确定翼型上的积冰研究区域，采用正交网格进行网格划分；
[0051]
根据当前积冰研究区域边界区分固体区域和流体区域；
[0052]
通过求解流动方程组，获取当前积冰研究区域中流体区域的速度场分布；
[0053]
输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程；
[0054]
对水滴收集过程中得到的已收集水滴进行积冰生长，积冰沿着当前积冰研究区域边界法向增长，计算积冰研究区域不同位置的积冰生长高度，对积冰研究区域积冰后的边界进行更新，得到积冰生长后新的积冰研究区域边界。
[0055]
本发明提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
[0056]
输入飞行器翼型，确定翼型上的积冰研究区域，采用正交网格进行网格划分；
[0057]
根据当前积冰研究区域边界区分固体区域和流体区域；
[0058]
通过求解流动方程组，获取当前积冰研究区域中流体区域的速度场分布；
[0059]
输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程；
[0060]
对水滴收集过程中得到的已收集水滴进行积冰生长，积冰沿着当前积冰研究区域
边界法向增长，计算积冰研究区域不同位置的积冰生长高度，对积冰研究区域积冰后的边界进行更新，得到积冰生长后新的积冰研究区域边界。
[0061]
本发明能够达到的有益技术效果是：
[0062]
本发明提出的结合浸润边界法的针对积冰问题冰形的模拟相较已有的冰形组合方法对于冰形样貌的模拟更加精细。另外由于浸入边界法采用正交网格，自始至终针对边界进行变化，因而不存在网格畸变等问题，报错概率低。
附图说明
[0063]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
[0064]
图1是本发明一实施例的流程图；
[0065]
图2是本发明一实施例中输入的飞行器模型的示意图；
[0066]
图3是本发明一实施例中翼型示意图；
[0067]
图4是对图3所示翼型进行网格划分的示意图；
[0068]
图5是图4机翼部分网格局部放大图；
[0069]
图6是识别区域类型示意图；
[0070]
图7是本发明一实施例中体积力源项求解网格示意图；
[0071]
图8是本发明一实施例中质量源项求解网格示意图。
[0072]
本发明目的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0073]
为了使本公开发明的目的、技术方案和优点更加清楚明白，下面结合具体实施例，并根据附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，未描述的内容以及部分英文简写为所属技术领域中普通技术人员所熟知的内容。以下实施例中给定的一些特定参数仅作为示范，在不同的实施方式中该值可以相应地改变为合适的值。
[0074]
参照图1，为了快速、准确描绘出积冰冰形样貌，本发明一实施例中提出一种飞行器积冰冰形模拟方法，包括以下步骤：
[0075]
s1，输入飞行器翼型，确定翼型上的积冰研究区域，采用正交网格进行网格划分；
[0076]
s2，根据当前积冰研究区域边界区分固体区域和流体区域；
[0077]
s3，通过求解流动方程组，获取当前积冰研究区域中流体区域的速度场分布；
[0078]
s4，输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程；
[0079]
s5，对水滴收集过程中得到的已收集水滴进行积冰生长，积冰沿着当前积冰研究区域边界法向增长，计算积冰研究区域不同位置的积冰生长高度，对积冰研究区域积冰后的边界进行更新，得到积冰生长后新的积冰研究区域边界。
[0080]
将得到的新的积冰研究区域边界用于下一次迭代。整个迭代过程结束后，基于最
终得到的积冰研究区域边界信息，可以精细化描绘出积冰冰形样貌，以及通过与初始积冰研究区域边界对比得到积冰研究区域不同位置处的积冰生长高度。基于积冰冰形样貌以及积冰生长高度，再结合积冰前后飞行器气动特性的变化，还可以分析得到积冰对飞行器气动特性的影响规律，进而指导后续的飞行器设计以及飞行器除冰问题的研究工作。
[0081]
浸入边界法(immersed boundary method)基本思路是：在cartesian结构网格中求解流动方程组时，在流动方程组中添加源项。这些源项反映了流动边界和流体的相互作用，通过求解带源项的流动方程组以数值求解随时间运动的物面和界面的流动问题。
[0082]
本发明一实施例的s1中，对积冰研究区域建立正交背景网格，采用正交网格进行网格划分。输入飞行器模型后，将其飞行器翼型作为积冰研究区域，为了提高后续识别过程的精确性可以针对积冰研究区域的正交网格进行适当加密。然后设置初始来流条件，如来流速度、密度等。参照图2，为本发明一实施例中输入的飞行器模型，图3为飞行器翼型。图4是对图3所示翼型进行网格划分的示意图；图5是图4机翼部分网格局部放大图。
[0083]
本发明一实施例的s2中，以机翼轮廓线作为初始积冰研究区域边界，在机翼轮廓线以内的区域为固体区域，机翼轮廓线以外的区域为流体区域。
[0084]
进一步地，在本发明一实施例中，s2中可以根据对初始积冰研究区域中各个网格单元的网格坐标点的属性进行判断，进而更精细的划分固体区域和流体区域。参照图6，以机翼轮廓线作为初始积冰研究区域边界，按照顺序依次对初始积冰研究区域中各个网格单元的网格坐标点的属性进行判断，方法如下：
[0085]
首先当前积冰研究区域边界以内的网格坐标点均为固体点；
[0086]
对于处于当前积冰研究区域边界上的网格单元，其位于当前积冰研究区域边界以外且最靠近当前积冰研究区域边界的网格坐标点，标记为力源项点；
[0087]
位于当前积冰研究区域边界以外的所有非力源项点的网格坐标点，均标记为流体点。
[0088]
此时所有网格点被分为三种类型，即固体点(标记为
‑
1)，流体点(标记为0)和力源项点(标记为1)，划分区域完成。
[0089]
在流动方程组添加源项有两种方式，第一种是在动量方程中添加力源项，第二种是在连续方程中添加质量源项。进而浸入边界法衍生出了两种算法：一类被称为连续力算法，离散前源项就有解析表达式，这种方法在生物流和多相流问题已得到成功运用；另一类被称为离散力算法，源项必须通过求解离散方程组得到，并且一般不是解析表达式，这种方法适合处理刚性壁面问题。
[0090]
本发明一实施例的s3中采用体积力源项和质量源项这两种源项结合起来的浸入边界法，即在流动方程组的动量方程中添加体积力源项，在流动方程组的连续方程中添加质量源项。
[0091]
对于积冰研究区域的流体区域内的任一网格单元，构建流动方程组，如下：流动方程组为：
[0092][0093]
其中u、v为网格单元其x、y方向的速度，r
e
为雷诺数，p为网格单元上的压力，q为添加的质量源项，f
x
、f
y
均为添加的体积力源项。
[0094]
求解流动方程组得到n步速度场
[0095]
动量方程组通过虚拟网格法求解，如图7所示，对于第i行、第j列的网格单元v
i,j
，网格单元v
i,j
其中心点p
i,j
位于当前积冰研究区域边界内部，其x、y方向的速度分别为u
g
,v
g
。根据无滑移边界条件，体积力源项的表达式如下：
[0096][0097][0098]
其中，δt为时间步，为已知量。时间步取前后两次迭代所花费的时间差，比如前一次迭代是1秒，当前次迭代是1.5秒，那δt取0.5秒。
[0099]
连续方程的质量源项求解过程如图8所示。设当前积冰研究区域边界与网格单元v
i,j
的两个顶点相交，(u
x,in
,v
y,out
)为当前积冰研究区域边界外部与网格单元v
i,j
边界垂直的速度分量；(u
x,out
,v
y,in
)为当前积冰研究区域边界内部与网格单元v
i,j
边界垂直的速度分量，q
i,j
为p
i,j
点处的质量源项，由质量守恒得：
[0100][0101]
其中，
[0102]
结合连续方程得到：
[0103][0104]
由此可得到：
[0105][0106]
差分格式为：
[0107][0108]
其中
ɑ
为网格单元边界被当前积冰研究区域边界所占据的比例，当网格单元边界
全部在当前积冰研究区域边界内部时，
ɑ
＝1。反之，当网格单元边界全部位于当前积冰研究区域边界外部时，
ɑ
＝0。
[0109]
将体积力源项、质量源项代入到流动方程组中，联立求解得到积冰研究区域的流体区域内的任一网格单元的u、v，得到积冰研究区域的流体区域内的任一网格单元的气流速度进而得到流体区域的速度场分布。
[0110]
本发明将积冰过程简化为两个子过程，即水滴收集过程和积冰生长过程。
[0111]
假设条件为：
[0112]
①
水滴近似为球形，且运动过程中形状不发生变化；
②
水滴碰撞前后温度、密度、粘性等性质不发生变化；
③
只考虑流场对水滴的影响，而不考虑水滴对于流场的扰动；
[0113]
本发明s3中输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程。
[0114]
针对水滴进行受力分析，仅考虑重力、浮力以及空气阻力的影响，由牛顿第二定律，采用拉格朗日方法得到水滴运动轨迹方程为：
[0115][0116]
上式中，m
d
为水滴质量，为水滴位移，ρ
a
、ρ
d
分别代表空气密度和水滴密度，v
d
为水滴体积，为当地重力加速度，a
d
代表水滴迎风面积，c
d
为阻力系数，为阻力系数，分别代表气流速度和水滴速度。
[0117]
然后根据计算得到的水滴移动轨迹信息并结合当前积冰研究区域边界进行碰撞检测，判断水滴是否接触积冰研究区域。将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程。接触当前积冰研究区域边界但未被收集的水滴标记为反弹水滴，利用水滴运动轨迹方程计算其反弹后的水滴运动轨迹，进一步进行碰撞检测。未接触当前积冰研究区域边界的水滴不予考虑。
[0118]
本发明s4中对水滴收集过程中得到的已收集水滴进行积冰生长。已经收集的水滴并非将全部积冰，根据质量守恒，对于当前积冰研究区域边界上的任一网格单元，其将凝结成冰的水滴质量等于流入该网格单元的水滴质量减去流出该网格单元的水滴质量。表达式为：
[0119][0120]
式中，表示落在某网格单元内的水滴质量总和，通过水滴运动轨迹可以确定。代表流入当前网格单元的水滴质量，代表由于蒸发或升华而损失的水滴质量，表示从当前网格单元流出的水滴质量，则单位时间内该网格单元内留下的水滴质量，即凝结成冰的质量。
[0121]
则第i个时间步内，当前积冰研究区域边界上的任一网格单元内的积冰生长高度为：
[0122][0123]
其中ρ
ice
表示冰的密度，a
s
代表网格单元面积。
[0124]
利用计算结果对积冰研究区域积冰后的边界进行更新。将积冰生长后新的积冰研究区域边界作为下一次迭代中的初始积冰研究区域边界输入，进行下一次的迭代。
[0125]
在本发明一实施例中，提供一种飞行器积冰冰形模拟装置，包括：
[0126]
第一模块，用于输入飞行器翼型，确定翼形上的积冰研究区域，采用正交网格进行网格划分；
[0127]
第二模块，用于根据当前积冰研究区域边界区分固体区域和流体区域；
[0128]
第三模块，用于通过求解流动方程组，获取当前积冰研究区域边界中流体区域的速度场分布；
[0129]
第四模块，用于输入水滴信息，获取水滴移动轨迹，利用水滴移动轨迹并结合当前积冰研究区域边界进行碰撞检测，将接触积冰研究区域且粘附在积冰研究区域的水滴作为已收集水滴，完成水滴收集过程；
[0130]
第五模块，用于对水滴收集过程中得到的已收集水滴进行积冰生长，积冰沿着当前积冰研究区域边界法向增长，计算积冰研究区域不同位置的积冰生长高度，对积冰研究区域积冰后的边界进行更新，得到积冰生长后新的积冰研究区域边界。
[0131]
上述飞行器积冰冰形模拟装置中的各组成模块的实现方法可以采用上述任一实施例中实现对应方法的方法实现，在此不再赘述。
[0132]
本发明一实施例中提供一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述任一实施例中所述飞行器积冰冰形模拟方法的步骤。
[0133]
本发明一实施例中提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一实施例中所述飞行器积冰冰形模拟方法的步骤。
[0134]
以上包含了本发明优选实施例的说明，这是为了详细说明本发明的技术特征，并不是想要将发明内容限制在实施例所描述的具体形式中，依据本发明内容主旨进行的其他修改和变型也受本专利保护。本发明内容的主旨是由权利要求书所界定，而非由实施例的具体描述所界定。