本发明涉及飞机飞行安全设计技术领域，特别是一种非稳态三维结冰数值模拟方法。
背景技术：
结冰严重威胁着飞机的飞行安全。对结冰问题的研究可以采用实验和计算这样两种方式。结冰实验所用的冰风洞造价及运行费用高昂，并且我国目前投入使用的冰风洞都是小型冰风洞，受其尺寸限制，某些工况是无法实验的，因此采用实验的方法研究结冰问题存在着诸多困难和障碍。而随着数值计算方法以及计算机技术的发展，数值模拟越来越成为结冰研究的一个重要手段。数值模拟能够以较低的成本，较快的速度对所研究问题进行定性或定量的分析。国外对结冰的数值模拟发展较早，并且已经具有一些较为成熟的计算软件。国内近年来对结冰的数值模拟也逐渐增多，二维结冰计算研究开展较全面，三维结冰计算也逐渐展开，三维结冰计算中存在三个显著的难点：滞止单元的确定、溢流水的流向和结冰表面网格遍历。另外目前的三维结冰计算是按准定常方法进行，往往将结冰过程作为一个时间步，这与真实结冰过程有一定差别，导致了明冰、混合冰的计算结果往往与实验误差较大。
技术实现要素：
因此，针对现有技术的上述问题，本发明为克服现有技术的不足，提出一种非稳态三维结冰的计算方法，能够更加真实地描述结冰过程，得到结冰冰形。
具体的，所述方法具体包括：
步骤1确定计算时间；
步骤2开始外循环，搜索当前计算单元组，首次计算时为滞止单元组；
步骤3开始内循环，在计算单元组内确定当前计算单元；
步骤4设定当前计算单元平衡温度；
步骤5计算当前计算单元的结冰量、溢流水量、平衡温度、冻结系数；
步骤6对当前计算单元进行结冰热力学计算，判断是否达到热力平衡；如果达到热力平衡，则进入步骤7，否则返回步骤4；
步骤7搜索溢流水将要流向的单元编号，计算流入水量及水温；
步骤8判断是否遍历当前计算单元组中所有单元，如果判断为是，结束内循环，进入步骤9，否则返回步骤3，进行下一次内循环；
步骤9判断是否遍历结冰表面所有网格单元，如果判断为是，结束外循环，进入步骤10，否则返回步骤2，进行下一次外循环并以步骤7搜索到的单元编号作为当前计算单元组。
步骤10判断是否达到设定的时间，如果判断为是，完成计算，否则返回步骤1开启下一时间步计算。
进一步的，所述方法中结冰表面网格的遍历是由内循环和外循环构成，外循环将计算单元组从前向后推进，内循环是完成计算单元组内网格的遍历。
进一步的，所述方法中滞止单元的确定方法为，根据单元的法向量与撞击速度矢量的夹角，如果小于某个数值，则该单元为滞止单元。
进一步的，所述方法中步骤7溢流水将要流向的单元的确定方法为，以当前单元的中心点与相邻单元的中心点形成矢量，点乘空气的速度矢量，获得两个矢量之间的夹角，如果夹角小于90度，则判定溢流水会流向该相邻单元。
进一步的，所述方法中步骤7计算流入水量的方法为根据两个矢量之间的夹角的大小比例分配，夹角越小，流入水量越大。
本发明针对三维结冰计算中存在的三个难点：滞止单元的确定、溢流水流向、结冰表面网格遍历，提出了如下解决方案：结冰表面网格的遍历采用两个循环过程，内循环和外循环。内循环负责当前网格单元组内的网格遍历，外循环负责从前缘向后的不断推进。分别计算了三维naca0012翼型的混合冰、明冰算例并与实验结果进行了对比，结果表明本发明提出的非稳态三维结冰计算方法是有效的。
附图说明：
图1为本发明三维结冰计算方法流程图；
图2为微元体的质量平衡图；
图3为微元体的能量平衡图；
图4为当前单元的相邻网格图；
图5为空气流速与网格中心点矢量夹角图；
图6为确定当前计算单元的示意图；
图7为温度-4.4℃冰形(工况1)；
图8为温度-28.3℃冰形(工况2)；
图9为温度-2.2℃冰形(工况3)；
图10为温度-3.89℃冰形(工况4截面a)；
图11为温度-3.89℃冰形(工况4截面b)
图12为温度-3.89℃冰形(工况4截面c)。
具体实施方式
下面对本发明的具体实施方式进行说明：
本发明提出的非稳态三维结冰计算方法由一个内循环和一个外循环组成。首先在结冰表面上根据水滴的撞击方向和计算单元法方向之间的夹角选取滞止单元群，由这些滞止单元(起始计算单元)出发，采用质量守恒和能量守恒对每个网格单元进行结冰热力学计算，获得当前单元的结冰量、溢流水量。根据当前网格单元中空气的流动方向确定溢流水的流向及流量，并标记下一次循环需要计算的网格单元及各网格单元流入的水量，如图1中的内循环，内循环完成时当前计算单元组中所有单元已被遍历，然后进入下一次当前单元组的计算，如图1中的外循环，如此不断地向后推进，直到遍历结冰表面上所有的网格单元，至此，完成一个时间步长的所有网格单元的结冰计算，而后进入下一个时间步长的循环。当达到设定时间时，计算完成。
气液两相流计算
结冰气象条件下的流场计算是空气、水滴共存的两相流计算。本发明采用欧拉—拉格朗日方法进行空气—液态水滴两相流计算。在欧拉—拉格朗日两相流计算方法中，空气相为连续相，直接求解时均n-s方程，湍流模型采用标准k-ε湍流模型；水滴相为离散相，通过拉格朗日法计算空气流场中水滴的运动轨迹方程。两相流计算前，对水滴运动做如下假设：a.气流不因含有过冷水滴而影响气流流过型面的流场；b.水滴在运动过程中，既不凝聚也不破碎，水滴的尺寸保持不变；c.温度、粘度、密度等介质参数在水滴运动过程中保持不变。当水滴穿过空气运动时，水滴的轨道通过当地空气作用于水滴上的各种平衡作用力来进行计算。
结冰热力学模型
水滴在壁面上的撞击、合并、迁移等过程十分复杂，为了方便计算，进行如下简化：
1.撞击到壁面上的水滴温度与壁面周围气体温度相等，由于过冷水滴直径很小，所以认为水滴与气流间换热迅速，保持温度相等；
2.撞击到壁面上的水滴全部停留在壁面上，即忽略撞击时溅离壁面的液体质量；
3.气流流过壁面时不携带起液体；
4.液体在壁面沿流线方向铺展，直到结冰或蒸发完毕；
5.结冰过程处于准稳态，任意时刻进入壁面的水滴总质量与结冰、蒸发以及流出的液体质量相等；
6.壁面和冰的换热量不计，即壁面为绝热壁面；
采用messinge的控制容积法对结冰表面收集的水建立质量守恒和能量守恒方程。
微元体的质量平衡如图2所示，对应的质量守恒方程可由下式确定：
为描述结冰量引入冻结系数，即表征每个单元的结冰量与进入该单元所有液态水质量之比，可以表达成：
微元体的能量平衡如图3所示，对应的能量守恒方程可由下式确定：
其中
联立方程(1)、(2)、(3)，并且假设表面温度，通过迭代求解可得微元体的结冰量
溢流水的分配
混合冰、明冰计算中，壁面某一网格单元上由于水没有完全结冰，未冻结的水在空气的作用下进入到相邻的网格中，这就是溢流问题。对于二维计算，水流向的下一个单元是非常明确的。对于三维计算，主要难点在于当前网格单元上未冻结的水到底流向哪个网格单元，如何确定该网格单元。这给结冰计算带来很大难度，也是本发明研究的核心问题之一。
水流动的动力主要是靠近壁面流体域单元中空气流动产生的拖曳力，因此本发明认为水的流动方向与壁面附近单元的空气流动方向是一致的。从某一网格单元(面网格)流出的水量分配到临近网格单元(面网格)的水量主要是根据几何分析得到。在三维非结构网格中，面上的网格为三角形，如图4，与之有边相邻的单元有3个(边界上是两个)，这里假设溢流水只流向与当前网格单元有边相邻的单元。水溢流方向与空气流动产生的拖曳力有关，因此水的溢流方向由空气的速度方向决定。如图5，当前单元的中心点o和与之相邻的单元的中心点o1形成矢量
滞止单元的选取
结冰计算首先要找到滞止单元，然后从滞止单元开始遍历结冰表面所有网格。滞止单元是起始计算单元，因此不存在上游溢流水流入的问题，即可认为流入水量为零。若是二维计算，滞止单元只有一个，较容易确定，但是三维计算中滞止单元有很多，如图6，单元1、3、5、7、9都是滞止单元，在程序中如何寻找、确定这些单元是一个较困难的问题，如果通过坐标寻找就会限制该方法的通用性。
本发明滞止单元的确定基于以下思想：首先根据单元的法向量与撞击速度的矢量二者之间夹角小于某个角度，确定一部分网格单元为滞止单元，如图6，1-20可能都被当作滞止单元，而这些网格并不都是真正的滞止单元，有些单元有溢流水流入，因此我们称这些单元为当前单元。首先计算这些当前单元上撞击水的结冰量，确定流出量和将要流入的单元的单元编号，溢流水流入的单元即为下一次循环的当前单元。
以上过程可以看出，在初始计算时，不需要精确地知道滞止单元都有哪些，只需要大概地选定一些网格单元作为起始计算单元即可，但是选定的这些网格单元必须全部包含真正的滞止单元。
结冰表面网格的遍历
t＝0时刻时，首先确定当前计算单元组，初次计算的当前单元组为滞止单元组，之后根据能量守恒和质量守恒原理，对当前单元组中每个网格单元进行结冰热力学计算，此时进行的是内循环过程，通过内循环获得当前单元组中每个单元的结冰量、溢流水量。根据以上溢流水的分配方法确定溢流水的流向及流量，并标记水流向的网格单元编号，也即下一次内循环需要计算的网格单元，下一次内循环要计算的网格单元组成了下一次计算的当前单元组，当前计算单元组中所有单元已被遍历后，本次内循环结束。然后进入外循环，搜索下一次当前单元组，之后再次进入内循环，如此不断地向后推进，直到遍历结冰表面上所有的网格单元。至此，一个时间段的所有网格单元的结冰计算完成，而后进入t＝t+δt的循环，当达到设定时间时，计算完成。
对比实验验证
为了与现有技术的模拟结果进行对比，本实施例截取计算获得的naca0012平直翼冰形、30°后掠角ms-317翼型冰形沿流线方向的二维截面，截面位于翼展中点。对于截面为glc-305的后掠翼，截取了垂直于前缘的截面a和截面b，以及平行于来流的截面c三个截面的冰形图与现有技术进行了对比，工况如表1所示。
表1
图7是温度为-4.4℃时naca0012翼型计算冰形与现有技术的对比，从图上可以看出计算冰形和实验结果基本符合，结冰上下极限、驻点附近最大结冰厚度和实验结果较一致，冰层生长趋势是向左上方突起，与实验结果吻合。
对于温度较低的霜冰，可认为溢流水为零，本实施例计算了温度为-28.3℃时naca0012翼型冰形，并与实验进行了对比，如图8，结冰上、下极限以及最大冰层厚度与实验吻合较好。
对于30°后掠角ms-317翼型，本实施例计算了-2.2℃的冰形，与现有技术对比如图9，结冰上极限比实验略大，结冰下极限和冰层厚度与实验一致，冰层的生长趋势与实验一致。
截面为glc-305的后掠翼中截面a、截面b、截面c的冰形与现有技术对比如图10至图12。图10中，结冰上、下极限比实验值略大，结冰厚度与实验基本一致，冰层的生长趋势与实验一致。图11中，结冰上、下极限以及最大结冰厚度与实验基本一致，但左上方突起的冰角模拟欠佳。图12中，计算冰形与实验非常相似，吻合较好。
通过以上对比看出本发明模拟方法的计算结果与实验结果基本一致，计算结果能够反映冰层的生长趋势和结冰最大厚度，说明本发明模拟方法三维结冰计算方法是有效的。
以上是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。