CN115659517B - 一种旋翼桨叶结冰准非定常数值模拟方法和系统 - Google Patents

一种旋翼桨叶结冰准非定常数值模拟方法和系统 Download PDF

Info

Publication number
CN115659517B
CN115659517B CN202211402273.XA CN202211402273A CN115659517B CN 115659517 B CN115659517 B CN 115659517B CN 202211402273 A CN202211402273 A CN 202211402273A CN 115659517 B CN115659517 B CN 115659517B
Authority
CN
China
Prior art keywords
icing
grid
time step
calculating
rotor blade
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211402273.XA
Other languages
English (en)
Other versions
CN115659517A (zh
Inventor
朱春玲
田成林
杨宇辰
赵宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202211402273.XA priority Critical patent/CN115659517B/zh
Publication of CN115659517A publication Critical patent/CN115659517A/zh
Application granted granted Critical
Publication of CN115659517B publication Critical patent/CN115659517B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种旋翼桨叶结冰准非定常数值模拟方法和系统，其中方法包括生成旋翼桨叶网格和背景网格；将总结冰时间均匀划分为多个结冰时间步；将旋翼的运动周期均匀划分为多个物理时间步；计算当前物理时间步的非定常流场；计算当前物理时间步的水滴场；计算每个结冰时间步的结冰计算周期总数；计算每个结冰时间；计算每个结冰时间步内结冰计算的次数；更新结冰后的旋翼桨叶网格；完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。本发明保留旋翼运动过程中流场、水滴场和结冰的非定常特性，适用于直升机在悬停、前飞状态下的旋翼结冰数值模拟，提高了旋翼桨叶冰形预测的准确度。

Description

一种旋翼桨叶结冰准非定常数值模拟方法和系统
技术领域
本发明属于结冰数值模拟技术领域，尤其涉及一种旋翼桨叶结冰准非定常数值模拟方法和系统。
背景技术
直升机在结冰条件下飞行时，旋翼表面会发生结冰现象。旋翼结冰会影响其气动性能，降低直升机的操纵性和稳定性，严重威胁直升机的飞行安全，因此开展旋翼结冰数值模拟研究对于保障直升机飞行安全至关重要。
直升机在悬停或者前飞过程中，流场和水滴场具有很强的非定常效应，旋翼结冰是一种非定常结冰现象。在旋翼结冰数值模拟研究中，应用较为广泛的是一种针对旋翼翼型的结冰计算方法：首先通过计算得到三维旋翼流场结果，然后沿桨叶展向提取二维旋翼翼型的流场结果，最终对旋翼翼型进行结冰热力学计算。通过这种方法可以简单地得到一些旋翼结冰特征，但是由于没有考虑旋翼三维效应以及离心力的影响。
另外，定常和准定常方法被提出并应用于旋翼结冰数值模拟上。例如，Rajmohan发展了一种旋翼前飞状态下的结冰数值模拟方法，搭建了耦合旋翼空气动力学与结冰热力学分析的直升机旋翼结冰计算平台，该计算平台基于OVERFLOW，TURNS和GT-Hybrid计算非定常流场，基于LEWICE3D代码进行水滴轨迹计算和结冰计算。该方法把结冰看作一个很慢的定常过程，把旋翼的三维流场信息按照桨叶方位角进行划分，并对每个特征方位角下的桨叶进行了单独的结冰计算。
陈希基于结构运动嵌套网格系统首先进行了旋翼非定常流场和水滴场数值模拟，然后采用考虑水膜流动的Messinger结冰热力学模型进行结冰计算。在旋翼流场和水滴场的基础上，通过线性插值来获得旋翼一个旋转周期内的结冰条件与方位角的对应关系，在一个旋转周期内通过对特征方位角上的桨叶进行结冰量计算并依次累加，最终获得一个旋转周期内的总结冰量，并反映在结冰外形上，该方法在得到一个旋转周期的结冰量后通过线性累加最终获得多个结冰周期的结冰量。
Myles Morelli采用开源软件SU2计算旋翼的非定常流场，采用基于拉格朗日法的内部代码PoliDrop计算每个非定常时间步的水滴运动轨迹和撞击位置，并在旋翼一个旋转周期后计算水收集系数，最后采用基于Myers结冰热力学模型的内部代码PoliMIce进行结冰计算。
上述定常或准定常结冰计算方法忽略了旋翼运动过程中的强非定常效应，无法准确预测冰形。
发明内容
本发明针对现有技术中的不足，提供一种旋翼桨叶结冰准非定常数值模拟方法和系统。
第一方面，本发明提供一种旋翼桨叶结冰的准非定常数值模拟方法，包括：
S1，根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格；
S2，将总结冰时间均匀划分为多个结冰时间步；
S3，将旋翼的运动周期均匀划分为多个物理时间步；
S4，获取当前物理时间步的旋翼桨叶网格；
S5，根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元；
S6，根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场；
S7，根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场；
S8，根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数；
S9，在旋翼的运动周期内计算每个结冰时间；
S10，计算每个结冰时间步内结冰计算的次数；
S11，在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形；
S12，采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格；
S13，在每个结冰时间步内重复步骤S3-S13，直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。
进一步地，所述根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场，包括：
根据以下公式计算当前物理时间步的非定常流场：
Figure 100002_DEST_PATH_IMAGE001
；
其中，
Figure 112308DEST_PATH_IMAGE002
为流场计算网格单元体守恒变量；
Figure 100002_DEST_PATH_IMAGE003
为流场计算网格单元面对流通量；
Figure 237259DEST_PATH_IMAGE004
为流场粘性通量；t为计算时间；Ω为计算网格单元的体积；S为计算网格单元其中一个 面的面积。
进一步地，所述根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场，包括：
根据以下公式计算当前物理时间步的水滴场：
Figure DEST_PATH_IMAGE005
；
其中，
Figure 409877DEST_PATH_IMAGE006
为水滴场计算网格单元体守恒变量；
Figure 52211DEST_PATH_IMAGE007
为水滴场计算网格单元面对流 通量；
Figure 656367DEST_PATH_IMAGE008
为水滴场计算网格单元的源项；t为计算时间；Ω为计算网格单元的体积；
Figure 100002_DEST_PATH_IMAGE009
为水 滴场计算网格单元面的法向量；S为计算网格单元其中一个面的面积。
进一步地，所述计算每个结冰时间步内结冰计算的次数，包括：
根据以下公式计算每个结冰时间步内结冰计算的次数：
k=j×m；
其中，k为每个结冰时间步内结冰计算的次数；j为每个结冰时间步所包含的结冰计算周期总数，j=T ice /T，T ice 为结冰时间步；T为旋翼的运动周期；m为小于物理时间步总数的常数。
第二方面，本发明提供一种旋翼桨叶结冰的准非定常数值模拟系统，包括：
网格生成模块，用于根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格；
第一划分模块，用于将总结冰时间均匀划分为多个结冰时间步；
第二划分模块，用于将旋翼的运动周期均匀划分为多个物理时间步；
第一获取模块，用于获取当前物理时间步的旋翼桨叶网格；
重叠网格装配模块，用于根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元；
第一计算模块，用于根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场；
第二计算模块，用于根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场；
第三计算模块，用于根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数；
第四计算模块，用于在旋翼的运动周期内计算每个结冰时间；
第五计算模块，用于计算每个结冰时间步内结冰计算的次数；
第二获取模块，用于在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形；
更新模块，用于采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格；
循环模块，用于在每个结冰时间步内依次重复第二划分模块至更新模块的操作，直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。
进一步地，所述第一计算模块，包括：
第一计算单元，用于根据以下公式计算当前物理时间步的非定常流场：
Figure 686640DEST_PATH_IMAGE010
；
其中，
Figure 100002_DEST_PATH_IMAGE011
为流场计算网格单元体守恒变量；
Figure 845089DEST_PATH_IMAGE012
为流场计算网格单元面对流通量；
Figure 100002_DEST_PATH_IMAGE013
为流场粘性通量；t为计算时间；Ω为计算网格单元的体积；S为计算网格单元其中一个面 的面积。
进一步地，所述第二计算模块，包括：
第二计算单元，用于根据以下公式计算当前物理时间步的水滴场：
Figure 320808DEST_PATH_IMAGE014
；
其中，
Figure 779471DEST_PATH_IMAGE006
为水滴场计算网格单元体守恒变量；
Figure 100002_DEST_PATH_IMAGE015
为水滴场计算网格单元面对流 通量；
Figure 980645DEST_PATH_IMAGE016
为水滴场计算网格单元的源项；t为计算时间；Ω为计算网格单元的体积；
Figure 32915DEST_PATH_IMAGE009
为水 滴场计算网格单元面的法向量；S为计算网格单元其中一个面的面积。
进一步地，所述第五计算模块，包括：
第三计算单元，用于根据以下公式计算每个结冰时间步内结冰计算的次数：
k=j×m；
其中，k为每个结冰时间步内结冰计算的次数；j为每个结冰时间步所包含的结冰计算周期总数，j=T ice /T，T ice 为结冰时间步；T为旋翼的运动周期；m为小于物理时间步总数的常数。
第三方面，本发明提供一种计算机设备，包括处理器和存储器；其中，处理器执行存储器中保存的计算机程序时实现第一方面所述的旋翼桨叶结冰的准非定常数值模拟方法的步骤。
第四方面，本发明提供一种计算机可读存储介质，用于存储计算机程序；计算机程序被处理器执行时实现第一方面所述的旋翼桨叶结冰的准非定常数值模拟方法的步骤。
本发明提供一种旋翼桨叶结冰准非定常数值模拟方法和系统，其中方法包括根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格；将总结冰时间均匀划分为多个结冰时间步；将旋翼的运动周期均匀划分为多个物理时间步；获取当前物理时间步的旋翼桨叶网格；根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元；根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场；根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场；根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数；在旋翼的运动周期内计算每个结冰时间；计算每个结冰时间步内结冰计算的次数；在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形；采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格；直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。本发明保留了旋翼运动过程中流场、水滴场和结冰的非定常特性，适用于直升机在悬停、前飞状态下的旋翼结冰数值模拟，提高了旋翼桨叶冰形预测的准确度。
附图说明
为了更清楚地说明本发明的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种旋翼桨叶结冰的准非定常数值模拟方法的流程示意图；
图2为本发明实施例提供的非结构重叠网格宿主单元搜索过程示意图；
图3为本发明实施例提供的三维空间线段和三角形位置关系示意图；
图4为本发明实施例提供的霜冰状态下控制体表面热流示意图；
图5为本发明实施例提供的明冰状态下控制体表面热流示意图；
图6为本发明实施例提供的三维曲面上的贴体坐标系示意图；
图7为本发明实施例提供的SC2110翼型俯仰振荡过程的结冰计算结果与文献参考值以及实验值的对比图；
图8为本发明实施例提供的旋翼翼型俯仰振荡过程中不同时刻的结冰计算冰形以及与实验值的对比图；
图9为本发明实施例提供的Schweizer 269旋翼液态水含量分布俯视图；
图10为本发明实施例提供的Schweizer 269旋翼单片桨叶的结冰分布图；
图11为本发明实施例提供的Schweizer 269旋翼单片桨叶在截面r=0.9R处的结冰冰形与实验结果以及Lewice计算结果的对比图；
图12为本发明实施例提供的UH-1H旋翼液态水含量分布俯视图；
图13为本发明实施例提供的计算得到的UH-1H旋翼三维结冰分布图(a)与Lewice计算结果(b)的对比图；
图14为本发明实施例提供的一种旋翼桨叶结冰的准非定常数值模拟系统的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在一实施例中，如图1所示，本发明实施例提供一种旋翼桨叶结冰的准非定常数值模拟方法，包括：
步骤S1，根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格。
利用网格生成软件，例如利用ICEM软件根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格，网格类型为非结构网格，网格文件类型为CGNS格式。
步骤S2，将总结冰时间均匀划分为多个结冰时间步。
将整个结冰周期T total 划分为N个结冰时间步，在每个结冰时间步T ice =T total /N内冰形保持不变，每经过T ice 后更新冰形。
步骤S3，将旋翼的运动周期均匀划分为多个物理时间步。
每个物理时间步Δt=T/n。
步骤S4，获取当前物理时间步的旋翼桨叶网格。
步骤S5，根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元。
每个桨叶网格与背景网格的非结构重叠网格装配过程如下：
基于相邻单元搜索算法和交替数字叉树（ADT）搜索算法相结合的方法搜寻网格单元的宿主单元，如图2所示。
（1）建立桨叶物面的二叉树数据结构，并初始化搜索的出发单元；
（2）将单元I形心与当前搜索单元J形心连线，基于相交算法判断线段IJ与当前搜索单元J各个面的相交情况：
1）如果无相交面，则判定当前单元J为单元I的宿主单元，结束搜索；
2）如果有相交面且该相交面不是边界面，则以该相交面的相邻单元为新的搜索单元，重新执行操作（2）；
3）如果有相交面且该相交面是除物面外的边界面，则判定单元I为活动单元，结束搜索；
4）如果有相交面且该相交面是物面，基于ADT算法判断该线段与整个物面的相交情况：如果相交面个数为奇数，则单元I的形心位于物面内，判定单元I为非活动单元，结束搜索；如果相交面个数为偶数，则以距离单元I最近的相交物面的相邻单元为新的搜索单元，重新执行操作（2）；
在所述的宿主单元搜索过程中，需要用到几何相交算法，非结构重叠网格系统中存在不同类型的网格单元，比如四面体、三棱柱、六面体等，这些网格单元的面由三角形或四边形组成，所以判断线段和单元面的相交情况，就是判断三维空间中一条线段和三角形或者四边形的相交情况：
在判断空间中线段和单元面的相交情况时，采用基于Plucker坐标系的相交判断算法；
三维空间中两个不同的点p = (px，py，pz)和q = (qx，qy，qz)可以定义一条有向直线，在Plucker坐标系中，这条有向直线L对应于6个系数：L=(L[0]，L[1]，L[2]，L[3]，L[4]，L[5])。定义L=plucker(p，q)，为了求这6个系数，首先构造一个2×4的矩阵：
Figure 876106DEST_PATH_IMAGE017
；
L中的每个系数等于这个矩阵的2×2子矩阵的行列式，即：
L[0]=pxqy-qxpy；
L[1]=pxqz-qxpz；
L[2]=px-qx；
L[3]=pyqz-qypz；
L[4]=pz-qz；
L[5]=qy-py；
如果a[6]和b[6]是基于Plucker坐标系定义的两条直线，通过Side Operator可以确定两条直线的相交关系：
Side（a，b）=a[0]×b[4]+a[1]×b[5]+a[2]×b[3]+a[3]×b[2]+a[4]×b[0]+a[5]×b[1]；
当且仅当Side（a，b）=0时，直线a和b平行或相交；
首先，基于Plucker坐标系和Side Operator判断线段所在直线和三角形的相交情况，由于单元I形心与当前搜索单元J的形心连线所在直线不可能与单元J的面共面，所以不考虑直线与三角形共面的情况，线段S的两个端点为p和q，三角形的三个顶点为v0、v1和v2，通过p和q建立直线L的Plucker坐标L=plucker(p, q)，通过v0、v1和v2分别建立三角形三条边的Plucker坐标：e1=plucker(v0，v1)，e2=plucker(v1，v2)和e3=plucker(v2，v0)，然后分别建立L和e1，e2，e3的Side Operator：
S1=side(L，e1)，S2=side(L，e2)，S3=side(L，e3)；
通过S1、S2和S3判定情况如下：
（1）如果S1、S2和S3均大于0或者均小于0，则直线与三角形相交（互相穿过），如图3中（a）所示；
（2）如果S1、S2和S3均不为0且符号不完全相同，则直线与三角形不相交，如图3中（b）所示；
（3）如果S1、S2和S3中任意两个为0，另一个不为0，则直线与三角形的一个顶点相交，如图3中（c）所示；
（4）如果S1、S2和S3中任意一个为0，另外两个符号相同，则直线与三角形的一条边相交，如图3中（d）所示；
如果直线与三角形位置关系的判定结果为（2），则线段与三角形也一定不相交，如果判定结果为（1）（3）（4），则需进一步地判断线段与三角形的相交情况：
首先构建Plucker坐标系下的直线L2，L3和L4，L2和L3是穿过线段端点和三角形公共顶点的直线，L4是公共顶点对面的边所在的直线，公共顶点的选择方法如下：
如果直线L和三角形相交（互相穿过）或者直线L和三角形的一条边相交，那么公共顶点可以从三角形的三个顶点任意选择；
如果直线L与三角形的一个顶点相交，那么公共顶点选择另外两个顶点中任意一个；
其次进行Side Operator运算，S4=side(L4，L3)，S5=side(L4，L2)；
通过S4和S5，对线段和三角形位置关系的判定情况如下：
（1）如果S4和S5异号，则线段和三角形不相交，如图3中（e）所示；
（2）如果S4和S5同号，则线段和三角形相交，如图3中（f）所示；
（3）如果S4和S5其中一个为0，则线段的其中一个端点在三角形上，为图3中（f）的一种特殊情况；
基于上述方法即可判断线段和三角形的相交情况，线段和四边形的相交判断算法与此类似，不再赘述。
在搜索单元I的宿主单元的过程中，初始出发单元的选择是影响搜索效率的关键，本发明提供了一种基于阵面推进技术的初始出发单元选择方法：
（1）确定任意网格单元的宿主单元，标记该单元并把该单元作为阵面链表的头节点；
（2）如果阵面链表非空，则取出链表的头节点I；
（3）循环单元I的邻居单元，如果邻居单元Jn尚未标记，则搜索Jn的宿主单元，搜索的初始出发单元设为单元I的宿主单元，搜索到Jn的宿主单元后，标记Jn并把其加入阵面链表中；
（4）删除阵面链表中的头节点，回到操作（2）；
通过以上方法可搜寻得到网格单元的宿主单元。
基于阵面推进技术计算桨叶网格单元形心的物面距，并将其作为单元物面距，对于不存在物面的背景网格预设参考物面距；
以物面距为分类准则参数，将网格单元分为活动单元、非活动单元和边界单元，其中，活动单元为参与流场计算的单元，非活动单元为不参与流场计算的单元，边界单元为负责网格间流场信息插值的单元；把单元I的物面距dI与其宿主单元J的物面距dJ进行比较，如果dI≤dJ，则单元I为活动单元，否则为非活动单元，基于所述准则判定全部单元的类型后，把与活动单元相邻的非活动单元重新标记为边界单元，记录边界单元及其对应的宿主单元用于流场和水滴场的重叠网格边界插值计算。
步骤S6，根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场。
采用双时间步法进行非定常流场求解，流场采用惯性坐标系下的任意拉格朗日-欧拉形式的非定常雷诺平均N-S方程，根据以下公式计算当前物理时间步的非定常流场：
Figure 330221DEST_PATH_IMAGE018
。
其中，
Figure 734920DEST_PATH_IMAGE019
为流场计算网格单元体守恒变量；
Figure 274485DEST_PATH_IMAGE020
为流场计算网格单元面对流通量；
Figure 921367DEST_PATH_IMAGE022
为流场粘性通量；t为计算时间；Ω为计算网格单元的体积；S为计算网格单元其中一个 面的面积。
基于格心格式的非结构网格有限体积法对流场控制方程进行离散求解：
采用HLLC格式计算对流通量，采用格林-高斯法计算流场变量梯度，采用左右单元值平均法计算粘性通量。采用非定常双时间步法，内迭代采用显式Runge-Kutta方法或隐式LU-SGS方法。
步骤S7，根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场。
水滴场采用基于欧拉法的水滴撞击特性求解方法，有如下假设：
（1）水滴为球形刚体，统一大小，不发生形变；
（2）水滴运动过程中，相互之间不发生碰撞、汇聚，撞击到壁面后不会飞溅；
（3）忽略水滴与周围空气间的传质、传热；
（4）水滴相动量方程中忽略粘性项及压力项，忽略湍流脉动的影响；
（5）作用在水滴上的力，仅考虑空气阻力与重力。
根据以下公式计算当前物理时间步的水滴场：
Figure 229989DEST_PATH_IMAGE023
。
其中，
Figure 304124DEST_PATH_IMAGE006
为水滴场计算网格单元体守恒变量；
Figure 330986DEST_PATH_IMAGE007
为水滴场计算网格单元面对流 通量；
Figure 781559DEST_PATH_IMAGE024
为水滴场计算网格单元的源项；t为计算时间；Ω为计算网格单元的体积；
Figure 944687DEST_PATH_IMAGE025
为水滴 场计算网格单元面的法向量；S为计算网格单元其中一个面的面积。
重复步骤S4-S7，计算2至3个旋翼的运动周期后，得到流场和水滴场的稳定周期解，选择m个点表征一个运动周期，记录下m个表征点的桨叶网格的物面流场、水滴场数据（坐标、法向量、空气速度、压力、剪切力、对流换热系数、水滴速度、水收集系数等）。
步骤S8，根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数。
步骤S9，在旋翼的运动周期内计算每个结冰时间。
步骤S10，计算每个结冰时间步内结冰计算的次数。
步骤S8-S10，计算每个结冰时间步T ice 所包含的结冰计算周期总数j=T ice /T，T为旋翼的运动周期，每个结冰计算周期中选择m个表征点进行结冰计算，计算每个表征点的结冰时间Δt ice =T/m，在每个结冰时间步T ice 内结冰计算的次数k=j×m；m为小于物理时间步总数的常数。
在每一次结冰计算中，读取所对应表征点的物面数据（坐标、法向量、空气速度、压力、剪切力、对流换热系数、水滴速度、水收集系数等），初始化冰层高度、水膜高度、冰水层表面温度等。
基于所述桨叶网格的物面数据进行当前表征点的结冰计算，结冰时间为Δt ice ，结冰计算中采用了考虑水膜流动的Myers结冰模型，考虑了离心力、空气压力梯度和剪切力对水膜流动的影响，Myers结冰模型包括水膜流动模型和结冰热力学模型，其中水膜流动模型通过计算连续性方程获得水膜高度，结冰热力学模型将冰层导热方程与能量方程联合求解，获得冰层高度。
在进行结冰数值模拟时，水滴撞击到温度远低于0℃的壁面时立刻冻结，结冰类型为霜冰，根据质量平衡方程确定结冰增长率，根据热量平衡方程确定冰层表面温度；水滴撞击到温度较高但仍低于0℃的壁面时只发生部分冻结，其余部分转化为水膜，在流动过程中逐渐冻结，结冰类型为明冰，冰层表面温度为0℃，根据能量平衡方程确定结冰增长率；
图4和图5展示了在霜冰和明冰状态下基于当前结冰模型旋翼桨叶表面的冰层、水膜控制体示意图，如图所示控制体内的热流被分为水滴撞击动能Q k ，气动加热热流Q a ，结冰冻结潜热Q l ，水滴显热Q d ，空气的对流换热热流Q c ，蒸发热流Q evap ，升华热流Q sub ，辐射热流Q r ，冰层向壁面导热热流Q cond ；
每次结冰计算时，结冰时间为Δt ice ，采用显式推进，时间步长取全局最小时间步dt，具体过程如下：
（1）更新冰层厚度和表面温度：
比较控制体表面温度T f 和T m =273.15K的大小：
如果T f ＜T m ，则当前为霜冰状态，根据质量平衡方程
Figure 699978DEST_PATH_IMAGE026
。
计算结冰增长率和冻结量，更新当前步冰层厚度
Figure 214135DEST_PATH_IMAGE027
：
Figure 202820DEST_PATH_IMAGE028
。
其中，
Figure 220455DEST_PATH_IMAGE029
为上一步冰层厚度，
Figure 901972DEST_PATH_IMAGE027
为当前步冰层厚度，m imp 表示撞击水量，m in 表示 上游溢流水量，m sub 表示升华质量，ρ i 表示冰的密度。
计算各项热流，并根据控制体内热量平衡方程
Q cond = Q k + Q a + Q l - Q c - Q d - Q sub - Q r 。
计算表面温度T f ：
Figure 637847DEST_PATH_IMAGE030
。
其中，T s 为壁面温度，k i 为冰层导热系数。
如果T f ≥T m ，则当前为明冰状态，计算各项热流，控制体内热量平衡方程为：
Q l = Q cond + Q d + Q c + Q evap + Q r - Q k - Q a 。
计算结冰增长速率：
Figure 430222DEST_PATH_IMAGE031
。
更新当前步冰层厚度
Figure 302363DEST_PATH_IMAGE027
：
Figure 656246DEST_PATH_IMAGE032
。
其中，L f 表示冰的融化潜热。
（2）更新水膜厚度：
图6为三维曲面上的贴体坐标系示意图，贴体坐标系用(s 1，s 2，η)表示。其中s 1，s 2与曲面相切，η是曲面在该点的法线，不同控制体的s 1，s 2和η的方向不同；
水膜在曲面上的流动受到空气压力梯度，质量力，空气剪切力等影响，在这里进行了如下假设：
1）水的粘度、密度在不发生相变时，随温度变化很小，在计算水膜流动时将其视为常数；
2）忽略机翼表面粗糙度对水膜流动的影响；
当前Myers模型基于润滑理论，简化了水膜在平板上的流动方程，得到了水膜在三维表面上沿向量s 1和s 2方向上单位长度的体积流量Q 1和Q 2，即水流速度u在该点水膜厚度上的积分：
Figure 879417DEST_PATH_IMAGE033
。
Figure 209904DEST_PATH_IMAGE034
。
其中，H w 表示水膜厚度，p表示空气压力，μ w 表示水的粘度，vf 1和vf 2分别表示质量力沿s 1和s 2方向上的分量，对于旋翼结冰问题，根据不同位置的惯性加速度计算出s 1和s 2方向上的离心力vf 1和vf 2，Sh 1和Sh 2分别表示气流剪切力沿s 1和s 2方向上的分量，由此得到水膜流动控制方程：
Figure 202131DEST_PATH_IMAGE035
。
式中，LWC表示来流液态水含量；V ∞ 表示空气的来流速度；β为水收集系数；m evap 为单位面积蒸发速率；ρ w 和ρ i 分别为水和冰的密度；H i 表示冰层厚度。
根据上式更新当前步水膜高度
Figure 225451DEST_PATH_IMAGE036
：
Figure 670339DEST_PATH_IMAGE037
。
其中，
Figure 804517DEST_PATH_IMAGE038
为上一步水膜厚度，
Figure 916829DEST_PATH_IMAGE039
为当前步水膜厚度；
（3）进入下一迭代时间步，重复（1）（2），直到达到结冰时间Δt ice 。
步骤S11，在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形。
每次结冰时间Δt ice 后记录冰层高度、水膜高度、冰层温度、水膜温度等数据，用来初始化下一次结冰计算的冰层高度、水膜高度、冰水层表面温度等。
重复进行结冰计算，直到完成k次结冰计算，累积结冰时间达到T ice 后记录并输出当前结冰时间步的冰形。
步骤S12，采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格。
采用动网格技术更新结冰后的桨叶网格的实现方法如下：
（1）提取并标记当前结冰时间步旋翼桨叶表面结冰后和结冰前的几何外形坐标数据，作为新、旧物面边界，标记桨叶网格远场边界；
（2）根据边界围成区域确定动网格生成区域，基于Delaunay准则对边界上的节点进行四面体划分，生成Delaunay背景图；
（3）基于阵面推进算法，寻找结冰前桨叶网格节点所属的背景图单元；
（4）根据结冰后的Delaunay背景图和结冰前的Delaunay背景图，检测各背景单元节点偏移情况；
（5）根据径向基函数，基于Delaunay背景图节点变形前后的位移变化，通过插值得到各网格节点位移s:
Figure DEST_PATH_IMAGE040
。
其中，M是Delaunay三角形的顶点数，三维情况下M=4；α i 是插值系数，φ是基函数，r i 是Delaunay三角形各顶点的坐标向量，r是当前节点的坐标向量。
系数α i 通过求解以下方程组获得：
dt=Mt,tα。
其中dt=[s(x1)，s(x2)，┄，s(xN)]T是Delaunay三角形节点的位移集合；Mt,t是一个M×M的矩阵，其元素为φ i,j =φ(||x i -x j ||)，即径向基函数，记为：
Figure 812848DEST_PATH_IMAGE041
。
其中，φ采用的一种径向基函数，记为：
φ（ξ）=(1-ξ)4(4ξ+1)。
对于Delaunay图中任意单元中节点v的位移
Figure DEST_PATH_IMAGE042
，通过α可写 为：
Δx v =Aαx；
Δy v =Aαy；
Δz v =Aαz；
A=[φ v,1,φ v,2,┄,φ v,N]。
对于网格节点v，其原始坐标为
Figure 604086DEST_PATH_IMAGE043
，则该点的新坐标
Figure DEST_PATH_IMAGE044
为：
Figure 73114DEST_PATH_IMAGE045
。
基于所述过程可以得到结冰后新的桨叶网格。
步骤S13，在每个结冰时间步内重复步骤S3-S13，直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。
下面给出几个算例作为本发明所公开方法的具体实施例。
实施例一
依据结冰风洞对Sikorsky SC2110翼型的俯仰振荡过程的结冰实验，采用本发明提出的准非定常结冰计算方法对SC2110翼型的俯仰振荡结冰过程进行了数值模拟研究。旋翼在直升机前飞过程中会进行周期性变距运动，在二维问题上反映为旋翼翼型的俯仰振荡运动，本实施例通过研究旋翼翼型俯仰振荡结冰过程来验证本发明方法的可行性和准确性。SC2110翼型弦长为0.381m，翼型以正弦运动形式进行俯仰振荡过程。本实施例结冰工况：结冰温度为259.15K，空气静压为101325Pa，来流速度为77.17m/s，初始迎角为5°，迎角脉动为±6°，振荡频率为2.8Hz，水滴平均容积直径为22μm，液态水含量为0.5g/m3，结冰时间为600s。总结冰时间T total =600s，设置总结冰步数N=5，每个结冰时间步长T ice =120s，根据振荡频率求得振荡周期T=0.357s，将每个运动周期T划分为180个非定常物理时间步，求得每个物理时间步长Δt=0.00198s，在一个振荡周期中选择m=4个表征点，分别为初始迎角位置以及最大、最小迎角位置，每个表征点的结冰时间Δt ice =0.08925s，每个结冰时间步长T ice 所包含的运动周期数j=T ice /T=336，在每个结冰时间步长内需要进行k=j×m=1344次结冰计算，每一次结冰计算的时间推进步长dt=10-4s。根据上述参数，采用准非定常结冰数值模拟方法进行计算，图7给出了本发明的冰形计算结果与文献参考值以及实验值的对比结果，可以看出，本发明计算结果的结冰极限、结冰量、冰角位置与文献参考值吻合良好，冰角位置比实验值略小。本实施例计算结果证明了本发明所提出方法的可行性和准确性。
实施例二
根据某俯仰振荡旋翼翼型的结冰实验结果，验证明冰情况下本发明方法的可靠性。该旋翼翼型弦长为0.5m，翼型以正弦运动形式进行俯仰振荡过程。本实施例结冰工况：结冰温度为263.05K，空气静压为101325Pa，来流速度为50m/s，初始迎角为-3°35′，迎角脉动为±3°，振荡频率为1Hz，水滴平均容积直径为20μm，液态水含量为1g/m3，结冰时间为900s。总结冰时间T total =900s，设置总结冰步数N=5，每个结冰时间步长T ice =180s，根据振荡频率求得振荡周期T=1s，将每个运动周期T划分为180个非定常物理时间步，求得每个物理时间步长Δt=0.00556s，在一个振荡周期中选择m=4个表征点，分别为初始迎角位置以及最大、最小迎角位置，每个表征点的结冰时间Δt ice =0.25s，每个结冰时间步长T ice 所包含的运动周期数j=T ice /T=180，在每个结冰时间步长内需要进行k=j×m=720次结冰计算，每一次结冰计算的时间推进步长dt=10-3s。根据上述参数，采用准非定常结冰数值模拟方法进行计算，图8给出了不同时刻冰形计算结果以及与实验值的对比结果，900s为最终计算冰形，可以看出本发明计算结果的结冰极限、冰角位置与实验值吻合良好，结冰量略小于实验结果。
实施例三
以AERTS试验台使用的小展弦比的Schweizer 269旋翼叶片为研究对象开展结冰数值模拟研究，来验证本发明所提出方法的准确性。Schweizer 269旋翼叶片截面为NACA0015翼型，截面恒定弦长为0.172m，叶片半径为1.18m，根切长度为20%半径长。本实施例结冰工况：结冰温度为262.45K，空气静压为101325Pa，旋翼转速为500rpm，水滴平均容积直径为25μm，液态水含量为2g/m3，结冰时间为180s。总结冰时间T total =180s，设置总结冰步数N=6，每个结冰时间步长T ice =30s，根据旋翼转速求得旋翼每转一圈的周期T=0.12s，将每个运动周期T划分为120个非定常物理时间步，即旋翼每3°转动一次，求得每个物理时间步长Δt=0.001s，在旋翼桨叶的一个运动周期中选择m=24个表征点，每个表征点的结冰时间Δt ice =0.005s，每个结冰时间步长T ice 所包含的运动周期数j=T ice /T=250，在每个结冰时间步长内需要进行k=j×m=6000次结冰计算，每一次结冰计算的时间推进步长dt=10-5s。根据上述参数，采用准非定常旋翼结冰数值模拟方法进行计算，图9给出了水滴场液态水含量LWC分布俯视图，在该工况下，受桨叶旋转影响，沿桨根向桨尖速度逐渐增加，水滴遮蔽区后移，桨叶运动区域LWC值沿桨根向桨尖方向逐渐减小。图10展示了小展弦比Schweizer 269旋翼单片桨叶的结冰分布图。本实施例由于旋翼展弦比较小，翼尖旋转速度偏小，因此旋翼旋转产生的气动加热效应和离心力等因素对最终结冰外型影响较小，旋翼结冰厚度沿展向近似呈线性分布，与试验结果展现的趋势一致。图10展示了截面r=0.9R处的结冰冰形与试验结果和Lewice结冰计算软件结果的对比图，可以看出，相比Lewice的计算结果，本发明计算结果的结冰厚度与结冰范围都与实验值更加吻合。
实施例四
以Bell UH-1H直升机旋翼为研究对象进行结冰数值模拟研究，进一步验证在大展弦比情况下本发明的鲁棒性。该旋翼截面为NACA0012翼型，恒定弦长为0.5334m，叶片半径为7.3152m，沿展向有10.9°的线性负扭转。本实施例结冰工况：结冰温度为254.15K，空气静压为101325Pa，旋翼转速为33.9rad/s，水滴平均容积直径为30μm，液态水含量为0.7g/m3，结冰时间为180s。总结冰时间T total =180s，设置总结冰步数N=6，每个结冰时间步长T ice =30s，根据旋翼转速求得旋翼每转一圈的周期T=0.185s，将每个运动周期T划分为120个非定常物理时间步，即旋翼每3°转动一次，求得每个物理时间步长Δt=0.00154s，在旋翼桨叶的一个运动周期中选择m=24个表征点，每个表征点的结冰时间Δt ice =0.00771s，每个结冰时间步长T ice 所包含的运动周期数j=T ice /T=162，在每个结冰时间步长内需要进行k=j×m=3888次结冰计算，每一次结冰计算的时间推进步长dt=10-5s。根据上述参数对UH-1H旋翼进行结冰数值模拟。图12给出了水滴场液态水含量LWC分布俯视图，可以看出当前工况下，由于桨尖和桨根马赫数相差较大，旋翼运动区域内LWC值变化更加明显，同样的变化反映在桨尖后侧遮蔽区范围上。图13给出了基于本发明计算方法得到的UH-1H旋翼三维结冰分布图(a)与Lewice结果(b)的对比，两者分布趋势一致。
基于同一发明构思，本发明实施例还提供了一种旋翼桨叶结冰的准非定常数值模拟系统，由于该系统解决问题的原理与前述旋翼桨叶结冰的准非定常数值模拟方法相似，因此该系统的实施可以参见旋翼桨叶结冰的准非定常数值模拟方法的实施，重复之处不再赘述。
在另一实施例中，本发明一个实施例提供的旋翼桨叶结冰的准非定常数值模拟系统，如图14所示，包括：
网格生成模块10，用于根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格。
第一划分模块20，用于将总结冰时间均匀划分为多个结冰时间步。
第二划分模块30，用于将旋翼的运动周期均匀划分为多个物理时间步。
第一获取模块40，用于获取当前物理时间步的旋翼桨叶网格。
重叠网格装配模块50，用于根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元。
第一计算模块60，用于根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场。
第二计算模块70，用于根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场；
第三计算模块80，用于根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数；
第四计算模块90，用于在旋翼的运动周期内计算每个结冰时间；
第五计算模块100，用于计算每个结冰时间步内结冰计算的次数；
第二获取模块110，用于在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形；
更新模块120，用于采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格；
循环模块130，用于在每个结冰时间步内依次重复第二划分模块至更新模块的操作，直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。
可选地，所述第一计算模块，包括：
第一计算单元，用于根据以下公式计算当前物理时间步的非定常流场：
Figure DEST_PATH_IMAGE046
。
其中，
Figure 134873DEST_PATH_IMAGE047
为流场计算网格单元体守恒变量；
Figure DEST_PATH_IMAGE048
为流场计算网格单元面对流通量；
Figure 703257DEST_PATH_IMAGE049
为流场粘性通量；t为计算时间；Ω为计算网格单元的体积；S为计算网格单元其中一个面的 面积。
可选地，所述第二计算模块，包括：
第二计算单元，用于根据以下公式计算当前物理时间步的水滴场：
Figure DEST_PATH_IMAGE050
。
其中，
Figure 247371DEST_PATH_IMAGE051
为水滴场计算网格单元体守恒变量；
Figure DEST_PATH_IMAGE052
为水滴场计算网格单元面对流通 量；
Figure 192194DEST_PATH_IMAGE053
为水滴场计算网格单元的源项；t为计算时间；Ω为计算网格单元的体积；
Figure 13519DEST_PATH_IMAGE025
为水滴场 计算网格单元面的法向量；S为计算网格单元其中一个面的面积。
可选地，所述第五计算模块，包括：
第三计算单元，用于根据以下公式计算每个结冰时间步内结冰计算的次数：
k=j×m。
其中，k为每个结冰时间步内结冰计算的次数；j为每个结冰时间步所包含的结冰计算周期总数，j=T ice /T，T ice 为结冰时间步；T为旋翼的运动周期；m为小于物理时间步总数的常数。
关于上述各个模块更加具体的工作过程可以参考前述实施例公开的相应内容，在此不再进行赘述。
在另一实施例中，本发明实施例还提供了一种计算机设备，包括处理器和存储器；其中，处理器执行存储器中保存的计算机程序时实现前述实施例公开的旋翼桨叶结冰的准非定常数值模拟方法。
关于上述方法更加具体的过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
在另一实施例中，本发明实施例还提供一种计算机可读存储介质，用于存储计算机程序；计算机程序被处理器执行时实现前述公开的旋翼桨叶结冰的准非定常数值模拟方法。
关于上述方法更加具体的过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似部分互相参见即可。对于实施例公开的系统、设备、存储介质而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本领域的技术人员可以清楚地了解到本发明实施例中的技术可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本发明实施例中的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）执行本发明各个实施例或者实施例的某些部分所述的方法。
以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。

Claims (10)

1.一种旋翼桨叶结冰的准非定常数值模拟方法，其特征在于，包括：
S1，根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格；
S2，将总结冰时间均匀划分为多个结冰时间步；
S3，将旋翼的运动周期均匀划分为多个物理时间步；
S4，获取当前物理时间步的旋翼桨叶网格；
S5，根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元；
S6，根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场；
S7，根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场；
S8，根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数；
S9，在旋翼的运动周期内计算每个结冰时间；
S10，计算每个结冰时间步内结冰计算的次数；
S11，在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形；
S12，采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格；
S13，在每个结冰时间步内重复步骤S3-S13，直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。
2.根据权利要求1所述的旋翼桨叶结冰的准非定常数值模拟方法，其特征在于，所述根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场，包括：
根据以下公式计算当前物理时间步的非定常流场：
Figure DEST_PATH_IMAGE001
；
其中，
Figure 929455DEST_PATH_IMAGE002
为流场计算网格单元体守恒变量；
Figure DEST_PATH_IMAGE003
为流场计算网格单元面对流通量；
Figure 908912DEST_PATH_IMAGE004
为流 场粘性通量；t为计算时间；Ω为计算网格单元的体积；S为计算网格单元其中一个面的面 积。
3.根据权利要求1所述的旋翼桨叶结冰的准非定常数值模拟方法，其特征在于，所述根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场，包括：
根据以下公式计算当前物理时间步的水滴场：
Figure 157491DEST_PATH_IMAGE006
；
其中，
Figure 896908DEST_PATH_IMAGE008
为水滴场计算网格单元体守恒变量；
Figure DEST_PATH_IMAGE009
为水滴场计算网格单元面对流通量；
Figure DEST_PATH_IMAGE011
为水滴场计算网格单元的源项；t为计算时间；Ω为计算网格单元的体积；
Figure DEST_PATH_IMAGE013
为水滴场计 算网格单元面的法向量；S为计算网格单元其中一个面的面积。
4.根据权利要求1所述的旋翼桨叶结冰的准非定常数值模拟方法，其特征在于，所述计算每个结冰时间步内结冰计算的次数，包括：
根据以下公式计算每个结冰时间步内结冰计算的次数：
k=j×m；
其中，k为每个结冰时间步内结冰计算的次数；j为每个结冰时间步所包含的结冰计算周期总数，j=T ice /T，T ice 为结冰时间步；T为旋翼的运动周期；m为小于物理时间步总数的常数。
5.一种旋翼桨叶结冰的准非定常数值模拟系统，其特征在于，包括：
网格生成模块，用于根据旋翼桨叶的几何模型，生成旋翼桨叶网格和背景网格；
第一划分模块，用于将总结冰时间均匀划分为多个结冰时间步；
第二划分模块，用于将旋翼的运动周期均匀划分为多个物理时间步；
第一获取模块，用于获取当前物理时间步的旋翼桨叶网格；
重叠网格装配模块，用于根据当前物理时间步的旋翼桨叶网格和背景网格进行非结构重叠网格装配，确定旋翼桨叶网格和背景网格的插值边界单元和插值边界单元对应的宿主单元；
第一计算模块，用于根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的非定常流场；
第二计算模块，用于根据当前物理时间步的非结构重叠网格装配结果计算当前物理时间步的水滴场；
第三计算模块，用于根据旋翼的运动周期计算每个结冰时间步的结冰计算周期总数；
第四计算模块，用于在旋翼的运动周期内计算每个结冰时间；
第五计算模块，用于计算每个结冰时间步内结冰计算的次数；
第二获取模块，用于在每个结冰时间步内依次获取对应结冰时间的物面流场和水滴场数据，直至完成每个结冰时间步内结冰计算的次数，得到当前结冰时间步的冰形；
更新模块，用于采用动网格技术更新结冰后的旋翼桨叶网格，将背景网格和更新结冰后的旋翼桨叶网格作为下一个结冰时间步的非定常流场和水滴场求解的计算网格；
循环模块，用于在每个结冰时间步内依次重复第二划分模块至更新模块的操作，直至完成所有数量次结冰时间步，得到总结冰时间内旋翼桨叶结冰后的冰形结果。
6.根据权利要求5所述的旋翼桨叶结冰的准非定常数值模拟系统，其特征在于，所述第一计算模块，包括：
第一计算单元，用于根据以下公式计算当前物理时间步的非定常流场：
Figure 85182DEST_PATH_IMAGE001
；
其中，
Figure 376486DEST_PATH_IMAGE002
为流场计算网格单元体守恒变量；
Figure DEST_PATH_IMAGE015
为流场计算网格单元面对流通量；
Figure DEST_PATH_IMAGE017
为流 场粘性通量；t为计算时间；Ω为计算网格单元的体积；S为计算网格单元其中一个面的面 积。
7.根据权利要求5所述的旋翼桨叶结冰的准非定常数值模拟系统，其特征在于，所述第二计算模块，包括：
第二计算单元，用于根据以下公式计算当前物理时间步的水滴场：
Figure 315623DEST_PATH_IMAGE018
；
其中，
Figure DEST_PATH_IMAGE019
为水滴场计算网格单元体守恒变量；
Figure 576840DEST_PATH_IMAGE020
为水滴场计算网格单元面对流通量；
Figure DEST_PATH_IMAGE021
为水滴场计算网格单元的源项；t为计算时间；Ω为计算网格单元的体积；
Figure 557303DEST_PATH_IMAGE022
为水滴场计 算网格单元面的法向量；S为计算网格单元其中一个面的面积。
8.根据权利要求5所述的旋翼桨叶结冰的准非定常数值模拟系统，其特征在于，所述第五计算模块，包括：
第三计算单元，用于根据以下公式计算每个结冰时间步内结冰计算的次数：
k=j×m；
其中，k为每个结冰时间步内结冰计算的次数；j为每个结冰时间步所包含的结冰计算周期总数，j=T ice /T，T ice 为结冰时间步；T为旋翼的运动周期；m为小于物理时间步总数的常数。
9.一种计算机设备，其特征在于，包括处理器和存储器；其中，处理器执行存储器中保存的计算机程序时实现权利要求1-4任一项所述的旋翼桨叶结冰的准非定常数值模拟方法的步骤。
10.一种计算机可读存储介质，其特征在于，用于存储计算机程序；计算机程序被处理器执行时实现权利要求1-4任一项所述的旋翼桨叶结冰的准非定常数值模拟方法的步骤。
CN202211402273.XA 2022-11-10 2022-11-10 一种旋翼桨叶结冰准非定常数值模拟方法和系统 Active CN115659517B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211402273.XA CN115659517B (zh) 2022-11-10 2022-11-10 一种旋翼桨叶结冰准非定常数值模拟方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211402273.XA CN115659517B (zh) 2022-11-10 2022-11-10 一种旋翼桨叶结冰准非定常数值模拟方法和系统

Publications (2)

Publication Number Publication Date
CN115659517A CN115659517A (zh) 2023-01-31
CN115659517B true CN115659517B (zh) 2023-02-28

Family

ID=85015926

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211402273.XA Active CN115659517B (zh) 2022-11-10 2022-11-10 一种旋翼桨叶结冰准非定常数值模拟方法和系统

Country Status (1)

Country Link
CN (1) CN115659517B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116562192B (zh) * 2023-07-06 2023-09-12 中国空气动力研究与发展中心计算空气动力研究所 一种飞机结冰冰形预测方法、装置、设备及存储介质
CN116738576B (zh) * 2023-07-06 2024-01-16 中国空气动力研究与发展中心计算空气动力研究所 一种旋翼结冰冰形预测方法、装置、设备及存储介质
CN118220493B (zh) * 2024-05-24 2024-08-13 天津云圣智能科技有限责任公司 一种可移动平台旋翼的结冰监测方法、程序产品及设备

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682145A (zh) * 2011-11-30 2012-09-19 天津空中代码工程应用软件开发有限公司 飞行结冰的数值模拟方法
WO2013078629A1 (zh) * 2011-11-30 2013-06-06 天津空中代码工程应用软件开发有限公司 飞行结冰的数值模拟方法
CN108460217A (zh) * 2018-03-13 2018-08-28 西北工业大学 一种非稳态三维结冰数值模拟方法
CN109376403A (zh) * 2018-09-29 2019-02-22 南京航空航天大学 一种基于笛卡尔自适应重构技术的二维结冰数值模拟方法
CN111396269A (zh) * 2020-06-08 2020-07-10 中国空气动力研究与发展中心低速空气动力研究所 一种多时间步非定常结冰计算方法、系统及存储介质
CN113792387A (zh) * 2021-09-30 2021-12-14 中国人民解放军国防科技大学 飞行器积冰冰形模拟方法、装置、计算机设备及存储介质
CN114462330A (zh) * 2022-01-19 2022-05-10 清华大学 飞机结冰冰形预测方法、装置、计算机设备和存储介质
CN115292806A (zh) * 2022-07-05 2022-11-04 北京航空航天大学 考虑周期性边界的三维热气防冰系统表面温度计算方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102522026B (zh) * 2011-11-29 2013-09-18 天津空中代码工程应用软件开发有限公司 飞行结冰模拟器

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682145A (zh) * 2011-11-30 2012-09-19 天津空中代码工程应用软件开发有限公司 飞行结冰的数值模拟方法
WO2013078629A1 (zh) * 2011-11-30 2013-06-06 天津空中代码工程应用软件开发有限公司 飞行结冰的数值模拟方法
CN108460217A (zh) * 2018-03-13 2018-08-28 西北工业大学 一种非稳态三维结冰数值模拟方法
CN109376403A (zh) * 2018-09-29 2019-02-22 南京航空航天大学 一种基于笛卡尔自适应重构技术的二维结冰数值模拟方法
CN111396269A (zh) * 2020-06-08 2020-07-10 中国空气动力研究与发展中心低速空气动力研究所 一种多时间步非定常结冰计算方法、系统及存储介质
CN113792387A (zh) * 2021-09-30 2021-12-14 中国人民解放军国防科技大学 飞行器积冰冰形模拟方法、装置、计算机设备及存储介质
CN114462330A (zh) * 2022-01-19 2022-05-10 清华大学 飞机结冰冰形预测方法、装置、计算机设备和存储介质
CN115292806A (zh) * 2022-07-05 2022-11-04 北京航空航天大学 考虑周期性边界的三维热气防冰系统表面温度计算方法

Also Published As

Publication number Publication date
CN115659517A (zh) 2023-01-31

Similar Documents

Publication Publication Date Title
CN115659517B (zh) 一种旋翼桨叶结冰准非定常数值模拟方法和系统
CN108460217B (zh) 一种非稳态三维结冰数值模拟方法
CN109376403B (zh) 一种基于笛卡尔自适应重构技术的二维结冰数值模拟方法
CN115563906B (zh) 一种基于非定常欧拉两相流的多步长结冰计算方法和系统
Aftosmis et al. Adaptation and surface modeling for Cartesian mesh methods
Pena et al. A single step ice accretion model using Level-Set method
Wang et al. Multi-body separation simulation with an improved general mesh deformation method
CN111396269B (zh) 一种多时间步非定常结冰计算方法、系统及存储介质
CN111310381B (zh) 一种三维水滴收集系数计算方法
CN112818573B (zh) 一种用于非结构网格的获取边界层非当地变量信息的方法
Prospathopoulos et al. Modelling wind turbine wakes in complex terrain
CN109751204A (zh) 一种风力机结冰数值模拟方法
WO2021087011A1 (en) System and method for simulating turbulence
Hedde et al. Improvement of the ONERA 3D icing code, comparison with 3D experimental shapes
Macarthur Numerical simulation of airfoil ice accretion
Wang et al. A new wind turbine icing computational model based on Free Wake Lifting Line Model and Finite Area Method
Donizetti et al. A three-dimensional level-set front tracking technique for automatic multi-step simulations of in-flight ice accretion
US20230185994A1 (en) Computational analysis of physical systems
BATINA Three-dimensional flux-split Euler schemes involving unstructured dynamic meshes
CN114692328A (zh) 一种风力发电机叶片结冰仿真方法及装置
Thompson et al. Discrete surface evolution and mesh deformation for aircraft icing applications
Murman et al. A vortex wake capturing method for potential flow calculations
Yoon et al. Effect of wing deformation by camber angle on aerodynamic performance of flapping micro air vehicles
CN111199093B (zh) 再入飞行器端头烧蚀的耦合方法、存储介质及终端
Li et al. A numerical study on rime ice accretion characteristics for wind turbine blades

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant