CN114940266A - 一种复杂防冰腔所能维持的蒙皮表面温度预测方法及系统 - Google Patents
一种复杂防冰腔所能维持的蒙皮表面温度预测方法及系统 Download PDFInfo
- Publication number
- CN114940266A CN114940266A CN202210302199.8A CN202210302199A CN114940266A CN 114940266 A CN114940266 A CN 114940266A CN 202210302199 A CN202210302199 A CN 202210302199A CN 114940266 A CN114940266 A CN 114940266A
- Authority
- CN
- China
- Prior art keywords
- surface temperature
- assumed
- skin
- heat
- temperature
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 84
- 238000004364 calculation method Methods 0.000 claims abstract description 127
- 238000012546 transfer Methods 0.000 claims abstract description 104
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 51
- 230000008569 process Effects 0.000 claims description 30
- 238000012937 correction Methods 0.000 claims description 21
- 238000004590 computer program Methods 0.000 claims description 10
- 239000007788 liquid Substances 0.000 claims description 9
- 230000004907 flux Effects 0.000 claims description 5
- 238000012360 testing method Methods 0.000 description 40
- 238000009826 distribution Methods 0.000 description 35
- 239000007789 gas Substances 0.000 description 33
- 239000010410 layer Substances 0.000 description 23
- 238000013461 design Methods 0.000 description 13
- 238000010586 diagram Methods 0.000 description 13
- 238000004088 simulation Methods 0.000 description 11
- 230000000694 effects Effects 0.000 description 10
- 239000000306 component Substances 0.000 description 9
- 230000002829 reductive effect Effects 0.000 description 9
- 239000007787 solid Substances 0.000 description 9
- 239000000243 solution Substances 0.000 description 9
- 238000010168 coupling process Methods 0.000 description 6
- 238000011156 evaluation Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 5
- 230000008878 coupling Effects 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 239000000463 material Substances 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000005507 spraying Methods 0.000 description 4
- 238000009529 body temperature measurement Methods 0.000 description 3
- 238000012938 design process Methods 0.000 description 3
- 238000010438 heat treatment Methods 0.000 description 3
- 238000009434 installation Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 230000036961 partial effect Effects 0.000 description 3
- 230000002265 prevention Effects 0.000 description 3
- 230000003068 static effect Effects 0.000 description 3
- 239000000126 substance Substances 0.000 description 3
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000004134 energy conservation Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 229910021389 graphene Inorganic materials 0.000 description 2
- 238000002347 injection Methods 0.000 description 2
- 239000007924 injection Substances 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 238000003801 milling Methods 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000001681 protective effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000007921 spray Substances 0.000 description 2
- 230000006641 stabilisation Effects 0.000 description 2
- 238000011105 stabilization Methods 0.000 description 2
- 239000010935 stainless steel Substances 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 229910000838 Al alloy Inorganic materials 0.000 description 1
- 229910001069 Ti alloy Inorganic materials 0.000 description 1
- RTAQQCXQSZGOHL-UHFFFAOYSA-N Titanium Chemical compound [Ti] RTAQQCXQSZGOHL-UHFFFAOYSA-N 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000000712 assembly Effects 0.000 description 1
- 238000000429 assembly Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 239000008358 core component Substances 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005206 flow analysis Methods 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 239000003595 mist Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000004080 punching Methods 0.000 description 1
- 238000004064 recycling Methods 0.000 description 1
- 239000002356 single layer Substances 0.000 description 1
- 229910001220 stainless steel Inorganic materials 0.000 description 1
- 229910001256 stainless steel alloy Inorganic materials 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000002352 surface water Substances 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
- 239000010936 titanium Substances 0.000 description 1
- 230000005514 two-phase flow Effects 0.000 description 1
- 239000002912 waste gas Substances 0.000 description 1
- 238000003466 welding Methods 0.000 description 1
Images
Classifications
-
- B—PERFORMING OPERATIONS; TRANSPORTING
- B64—AIRCRAFT; AVIATION; COSMONAUTICS
- B64D—EQUIPMENT FOR FITTING IN OR TO AIRCRAFT; FLIGHT SUITS; PARACHUTES; ARRANGEMENT OR MOUNTING OF POWER PLANTS OR PROPULSION TRANSMISSIONS IN AIRCRAFT
- B64D15/00—De-icing or preventing icing on exterior surfaces of aircraft
- B64D15/20—Means for detecting icing or initiating de-icing
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Aviation & Aerospace Engineering (AREA)
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
- Investigating Or Analyzing Materials Using Thermal Means (AREA)
Abstract
本发明涉及一种复杂防冰腔所能维持蒙皮表面温度的预测方法，包括步骤：S1对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和内表面温度；S2基于假定外表面温度、内表面温度计算蒙皮导热量；S3基于假定外表面温度计算外表面的第一总热流；S4基于假定内表面温度计算内表面的第一对流换热系数，并基于第一对流换热系数计算内表面的第一对流换热量；S5分别判断第一总热流、第一对流换热量是否等于导热量；当判断出第一总热流、第一对流换热量均等于导热量时，执行S6输出假定内、外表面温度，否则，执行S7对假定内、外表面温度进行修正，并基于修正后的蒙皮内、外表面温度返回到步骤S1。本发明的预测方法具有较高计算效率和计算精度。
Description
技术领域
本发明涉及飞机除冰领域，具体涉及一种复杂防冰腔所能维持的蒙皮表面温度预测方法。
背景技术
在飞机除防冰系统设计中，如何用最低的加热功率维持机翼表面一定的温度是一个关键性问题，传统的基于理论模型的工程估算方法越来越不使用于愈加复杂的机翼内部结构，例如，复合材料、石墨烯材料等新型材料的应用为采用理论建模估算飞机除防冰热载荷，或者估算蒙皮表面温度带来了极大的困难。
热气防冰是飞机结冰防护中最常采用的方法之一，热气防冰的其中一个主要防冰模块为防冰腔，在防冰腔的结构设计过程中，需要开展防冰腔的综合性能评估，主要用于检验防冰腔的结构设计是否合理以及是否达到设计要求。为了能够在设计过程中提前评估防冰腔的结构是否合理或者是否达到设计要求，可以通过预测防冰腔在引气量一定时能够维持的蒙皮表面温度来评估防冰腔的防冰性能，一般来说，当引气量一定时，防冰腔所能维持的蒙皮表面温度越高，防冰腔的防冰性能越好。
随着飞机除防冰领域的发展，对防冰性能的要求也进一步提升，防冰腔结构的设计也随之变得更加复杂。然而，随着防冰腔的结构越来越复杂(如防冰腔内设有笛形管，笛形管上可以设有不同间距、夹角或大小的射流孔)，现有技术中用于蒙皮表面温度的计算或仿真方法的适应性逐渐变差。通常为了追求计算精度，往往会在计算过程中考虑到多个影响参数，但这些参数的计算难度较大或者计算过程较为复杂(如计算过程中的迭代步骤较为繁琐)，例如，转捩处的临界雷诺数、转捩处以焓厚度表示的雷诺数、以弧长为特征长度的雷诺数等，由此使得计算流程较为复杂，计算效率也较低，并且当遇到复杂防冰腔时，由于考虑的参数过多而难以建立相适应的计算体系。
因此，目前亟需一种简单且精度高的复杂防冰腔所能维持的蒙皮表面温度计算方法。
发明内容
为了部分地解决或部分缓解上述技术问题，本发明基于内外流联合仿真可以快速预测防冰腔结构在一定量的引气量条件下，防冰腔所能维持的蒙皮表面温度，或者，基于某一特定的复杂防冰腔结构，计算选定最佳的防冰腔引参数(引气速率或引气量)，具有工程实用性。
本发明基于结冰数值模拟法预测防冰表面温度，主要分为外部热问题模拟和内部热问题模拟。外部热问题分析包括外部空气流场的计算分析、外部水滴运动的计算分析、蒙皮外表面热流分析及防冰热载荷计算等，首先求解空气流场，然后对过冷水滴运动方程进行求解，得出物面水滴收集系数，进而根据质量守恒方程和热平衡方程对物面的热力学过程进行求解，最后求解出蒙皮的表面温度。然后对防冰腔的防冰性能进行评估，例如，根据物面防冰温度要求，将表面温度需求作为边界条件，耦合求解所需的物面传递输入热量，即为防冰系统(或防冰腔)在当地需要提供的能量，以此可以确定如石墨烯多层复杂结构的防冰系统的功率分布，并结合防冰系统所要提供的能量和功率分布对防冰腔的防冰性能进行评估。
本发明第一方面提供了一种复杂防冰腔所能维持的蒙皮表面温度预测方法，包括步骤：
S1：对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和假定内表面温度；
S2：基于所述假定外表面温度和所述假定内表面温度计算得到蒙皮的导热量；
S3：基于所述假定外表面温度计算外表面的第一总热流；
S4：基于所述假定内表面温度计算内表面的第一对流换热系数，并基于所述第一对流换热系数计算内表面的第一对流换热量，其中，所述蒙皮的非驻点区域的所述第一对流换热系数的计算公式为：
S5：判断所述第一总热流是否等于所述导热量以及所述第一对流换热量是否等于所述导热量；当判断出所述第一总热流等于所述导热量，且所述第一对流换热量等于所述导热量时，执行S6：输出此时的所述假定外表面温度和假定内表面温度，否则，执行S7：对所述假定外表面温度和假定内表面温度进行修正，基于修正后假定外表面温度和假定内表面温度返回到步骤S1。
在一些实施例中，所述修正系数的计算公式为：
在一些实施例中，所述边界层雷诺数的计算公式为：
在一些实施例中，所述S3还包括步骤：
基于所述假定外表面温度计算外部流场，基于所述外部流场确定蒙皮外表面的第二对流换热量；
基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流。
在一些实施例中，所述基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流的步骤，包括步骤：
基于所述第二对流换热量，采用欧拉法计算水滴撞击特性，得到水滴撞击特性结果；
基于所述水滴撞击特性结果，采用基于改进Messinger模型部件表面热载荷的计算方法计算蒙皮外表面的所述各项热流；
基于所述各项热流计算所述第一总热流。
本发明还提供了一种用于预测复杂防冰腔所能维持的蒙皮表面温度的预测系统，包括：
温度假定模块，被配置为用于对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和假定内表面温度；
导热量计算模块，被配置为基于所述假定外表面温度和所述假定内表面温度计算所述蒙皮的导热量；
外表面总热流计算模块，被配置为基于所述假定外表面温度计算外表面的第一总热流；
内表面对流换热量计算模块，被配置为基于所述假定内表面温度计算内表面的第一对流换热系数，并基于所述第一对流换热系数计算内表面的第一对流换热量，其中，所述蒙皮的非驻点区域的所述第一对流换热系数的计算公式为：
第一判断模块，被配置为用于判断所述第一总热流是否等于所述导热量，以及所述第一对流换热量是否等于所述导热量；
温度输出模块，用于当所述第一判断模块判断出所述第一总热流等于所述导热量，以及所述第一对流换热量等于所述导热量时，输出此时的所述假定外表面温度和所述假定内表面温度；
温度修正模块，用于当所述第一判断模块判断出所述第一总热流不等于所述导热量，或者所述第一对流换热量不等于所述导热量时，对所述假定外表面温度和所述假定内表面温度进行修正，并将修正后的所述假定外表面温度和所述假定内表面温度发送到所述温度假定模块。
在一些实施例中，所述修正系数的计算公式为：
其中χ表示修正系数，θ为蒙皮的表面节点的两侧线段夹角。
在一些实施例中，所述边界层雷诺数的计算公式为：
在一些实施例中，所述外表面总热流计算模块包括：
外部流场计算单元，被配置为基于所述假定外表面温度计算外部流场，并基于所述外部流场计算蒙皮外表面的第二对流换热量；
外部热流计算单元，被配置为基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流。
在一些实施例中，所述外部热流计算单元包括：
第一计算子单元，被配置为基于所述第二对流换热量，采用欧拉法计算水滴撞击特性，得到水滴撞击特性结果；
第二计算子单元，被配置为基于所述水滴撞击特性结果，采用基于改进Messinger模型部件表面热载荷的计算方法计算蒙皮外表面的所述各项热流；
第三计算子单元，被配置为基于所述各项热流计算所述第一总热流。
本发明还提供了一种用于预测复杂防冰腔能够维持的蒙皮表面温度的计算机程序产品，其中存储有计算机程序，所述计算机程序被处理器执行时，控制所述计算机程序产品所在设备执行如上述实施例中任一所述方法的步骤。
有益技术效果：
相较于传统的多参数的计算方式，本发明提供了一种优化的表面温度预测方法，在计算非驻点区域的对流换热系数时，本发明在计算过程中不对非驻点区域进行区分，而是将非驻点区域直接作为湍流区域进行计算，并考虑到了蒙皮的形状结构上对蒙皮表面流场的影响，从而对预测方法进行了简化，并对预测方法进行修正，显著地提高了预测效率，并且经过实验数据验证，本发明的预测结果和试验结果相关性和一致性良好，准确性高(本发明的计算误差在工程应用上可以忽略不计，或者说对工程应用的实际影响很小)。
因此，基于本发明的预测方法或系统能够对蒙皮表面温度进行高效且高精度的预测，尤其是能够对复杂防冰腔进行可靠的防冰性能评估，由此可以在防冰腔的结构设计初期对防冰腔的结构进行防冰性能评估，减小对风洞试验的依赖性，缩短防冰腔结构的设计周期，并减少了防冰腔设计过程中的设计成本。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
图1是飞机热气防冰系统的工作原理图；
图2是防冰腔的典型结构示意图；
图3a是防冰腔展向切面的六种典型型式的结构示意图；
图3b是笛形管结构示意图；
图4a为本发明一示例性实施例中的方法流程示意图；
图4b为本发明一示例性实施例中的系统的装置结构示意图；
图5a是某型飞机的翼段几何外形示意图；
图5b为本发明一示例性实施例中选用的计算网格示意图；
图6a为三维翼段的表面水收集系数分布云图；
图6b为图6a所示的三维翼段展向50％剖面水收集系数曲线图；
图7a为三维翼段表面防冰热载荷分布云图；
图7b为图7a所示三维翼段的展向50％剖面的防冰热载荷分布曲线示意图；
图8为本发明一具体实施例中的蒙皮表面温度计算流程示意图；
图9a为防冰腔模型示意图；
图9b为图9a所示的防冰腔模型的网格划分示意图；
图10a为防冰腔表面压力分布示意图；
图10b为中间喷口截面Ma数分布示意图；
图10c为防冰腔内流线示意图；
图10d为防护区上、下表面水膜厚度分布示意图；
图10e为防护区上、下表面积冰增长率分布示意图；
图10f为防护区的蒙皮表面温度分布示意图；
图11a为结冰试验模型示意图；
图11b为结冰风洞试验段模型安装示意图；
图12为截面1#上温度传感器位置分布示意图；
图13a为防冰试验的第一角度的结果图；
图13b为防冰试验的第二角度的结果图；
图14a为中截面前缘点防冰过程温度随时间变化图；
图14b为防冰稳定后三截面温度分布图；
图15a为表面温度分布计算与试验结果对比图；
图15b为防冰效果测量示意图；
图16为流场中某一空间区域的区域划分示意图。
其中，P1为发动机的压气机，P2为流量限制器，P3为单向阀，P4为机翼防冰阀，P5为机翼前缘防冰腔；1为蒙皮，2为墙，3为波纹板，4为梁，5为挡板，6为笛形管(分配管)，8为射流孔，11a为上连接件，11b为截面1#，11c为截面2#，11d为截面3#，11e为下连接件，11f为前缘部段，11g为后缘部段；
02为温度假定模块，04为导热量计算模块，06为外表面总热流计算模块，08为内表面对流换热量计算模块，10为第一判断模块，12为温度输出模块，14为温度修正模块。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本文中，使用用于表示元件的诸如“模块”、“部件”或“单元”的后缀仅为了有利于本发明的说明，其本身没有特定的意义。因此，“模块”、“部件”或“单元”可以混合地使用。
本文中，术语“上”、“下”、“内”、“外”“前端”、“后端”、“两端”、“一端”、“另一端”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。
本文中，除非另有明确的规定和限定，术语“安装”、“设置有”、“连接”、“相连”等，应做广义理解，例如“相连”，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是无线连接，也可以是无线通信连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
本文中，术语“非驻点区域”指的是在蒙皮表面的计算区域(即计算域)中，除了驻点处的其他区域。
为了使本发明的目的、技术方案及优点更加清楚明白，首先对飞机的热气防冰工作原理进行说明。
飞机热气防冰系统工作原理如图1所示，当机翼防冰阀P4打开时，热空气由发动机的压气机P1引出，经过流量限制器P2及单向阀P3，进入机翼的笛形管，再通过管壁上的孔喷入机翼前缘防冰腔P5，防冰腔内热气沿通道的流动过程将热量传递到蒙皮，使防冰表面温度高于冰点，保证表面不产生结冰。单向阀P3的作用是当发动机损坏时，避免高压气流向发动机倒流；流量限制器P2的作用是限制进入防冰系统的热气流量，避免过高的空气量流入防冰腔导致发动机进气不足从而影响发动机的性能。
其中，防冰腔是热气防冰系统中的重要组成部分，是既有传热又有传质的热交换器。对于翼型来讲，前缘防冰腔的结构型式对防冰的效果影响很大，理想的防冰腔设计是希望热空气与蒙皮间的热交换足够大，即希望热气与防冰通道间的对流换热系数要大，防冰腔向外的传热面积要尽可能大。
防冰腔的结构构型很多，现有的飞机防冰腔主要形式如图2所示，图中方案(a)为最简单的防冰腔，热气在空腔A中流动，由于A腔的横断面积较大，造成热气的流速很小，使得热气与防冰表面间的换热系数也很小，即热空气与外蒙皮1之间的换热情况很差，此外，通过墙2的热损失比较大。这种型式的防冰腔热效率低，一般只能用于尺寸很小的部件(如发动机叶片等)，对于尺寸较大的机(尾)翼很少采用。
方案(b)在A腔内放一笛形管，用于热气的分配，这样可以减少防冰热气的能量损失，换热情况比方案(a)好，但换热效果仍有较大的提升空间。采用方案(b)型式的防冰腔，防冰热气通过波纹通道把热量传给外蒙皮，由于防冰通道的流通面积小，流速大，放热系数也比较大，这种防冰腔的缺点是沿展向传热面积比较小，因为防冰通道外无热气通过(图2中的方案(a))。
为此，提出了方案(c)，它是在方案(b)的基础上加了挡板5，使其展向切面如图2中的方案(b)所示，其传热面积得到有效增加，但结构也变复杂了。
方案(d)是在方案(a)的防冰腔中加入了一根笛形管6，热气通过笛形管射流，可增强与前缘蒙皮的换热，通过这种简易的形式可以有效的提高换热效率，但仍有很大的提升空间。
方案(e)称为“微引射”式，相较于方案(d)增加了一波纹板3，这种方案热气从笛形管的射流孔喷出，热气引射防冰后的空气，使得部分防冰后的“废气”循环使用，这样可节省防冰所需的热气，提高热能利用效率。
方案(f)是利用分配管喷出的热气使A腔增压，利用A腔与B腔(B腔为墙2、梁4以及波纹板3围合成的区域)的压差使气体流动，这种型式的热量利用率也较高，但结构也比较复杂。
参见图3a，为了实现更好的换热防冰效果，防冰腔应该具有尽可能大的传热面积。由图3a可知，a)为单层波纹壁通道，铆接处无热气流通，其传热面积比较小，防冰效果会受到较大影响；b)和c)分别为双层波纹壁通道和铆接处冲窝的双壁通道，这两种形式解决了铆接处无热气流道的问题，增大了传热的面积，有利于更多热量传导至蒙皮表面；随着加工工艺的不断发展，出现了内蒙皮铣槽铆接、外蒙皮铣槽铆接、外蒙皮铣槽点焊等形式的热气流通通道，这样可以大大提高热利用率，同时减轻防冰腔的重量。
笛形管是热气防冰系统中另一核心部件，其工作原理主要是利用管内外的高压比使得气流在射流孔处达到音速，从而控制热量和流量的分配，以满足防冰系统所需的防冰热载荷。如图3b所示，笛形管形如一根长笛，其表面分布有各种排布形式的射流孔8。笛形管通常横置于防冰腔内前缘，为了保证下游射流压力，笛形管8沿着展向通常采用渐变管径设计。
笛形管参数的选择对热气防冰效率影响较大，其主要包括：笛形管材料、笛形管厚度、笛形管内径，射流孔间距以及射流夹角等。通常在工程上，笛形管材料一般选择金属钛、不锈钢或铝合金等；笛形管厚度一般在1-2mm；射流孔直径一般在0.25mm-3mm之间；考虑到热气利用效率以及蒙皮表面加热的均匀性，射流孔间距与孔径的比值一般在3.1-12.5之间。
防冰腔以及笛形管的具体参数要根据防冰系统防冰热载荷需要以及实际加工情况而定，对于不同的模型，其考虑不尽相同，如模型展长较大就会导致笛形管较长，笛形管各参数需要与之匹配，尽量使热效率较高。
为了确定防冰腔以及笛形管的具体参数(设计参数)需要对防冰腔的防冰性能进行综合评估，然而随着防冰腔结构的复杂化，对防冰腔所能维持的蒙皮表面温度的预测难度逐渐增大，现有计算方法的预测误差也随之增大，难以应用于工程实际。
本文中，“复杂防冰腔”指的是内部结构多重、复杂、异型、非常规的防冰腔结构，这一类型的防冰腔的内部结构往往数量较多(如图2中的方案(d)中在原方案(a)基础上增设了笛形管等结构，又如方案(e)在原方案(a)基础上增设了笛形管、波纹板等结构)，并且与此类防冰腔内部结构搭配的外形结构在某些实施例中也相对较为复杂(例如，图2中的方案(b)、方案(d)、方案(e)等)。
实施例一
参见图4a，本发明第一方面提供了一种复杂防冰腔所能维持的蒙皮表面温度预测方法，包括步骤：
S1：对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和假定内表面温度；
优选地，在一些实施例中，外表面温度和内表面温度的假定数值在一定范围内，例如，如初次对内、外表面温度进行假定时，假定外表面温度的取值范围在大约2°-大约5°之间，假定内表面温度的取值范围在大约100°-大约150°之间，当然，内外表面温度的假定数值的选取范围可以基于实际情况，由工作人员基于历史工作经验或除冰需求进行自适应的设定，例如，在一些实施例中，对防冰需求较高时，外表面温度的取值范围可以在大约2°-15°之间。
S2：基于假定外表面温度和假定内表面温度计算得到蒙皮的导热量；
在一些实施例中，导热量计算过程中选用的公式包括：
这里K表示蒙皮的导热系数，Cps表示蒙皮的定压比热，ρs表示蒙皮的密度，ΔT表示相邻两个迭代步之间(即当前迭代步和上一个迭代步之间)的温度差，ΔH表示相邻两个迭代步之间的焓。
求解上述方程时，从初始温度和初始焓开始迭代，新的温度采用下式计算：
Tnew＝Told+ΔT (2a)
其中，Tnew表示当前迭代步骤得到的温度，Told表示上一个迭代步骤的温度，ΔT表示相邻两个迭代步(即当前迭代步和上一个迭代步)之间的温度差。
对应的焓的计算公式为：
Hnew＝Hold+ΔH＝Hold+CpsΔT (2b)
其中，Hnew表示当前迭代步骤得到的温度对应的焓值，Hold表示上一个迭代步骤得到的温度对应的焓值，ΔT表示相邻两个迭代步之间的温度差，ΔH表示相邻两个迭代步之间的焓。
当达到稳态，即计算结果收敛时，ΔT趋于极小值。
S3：基于假定外表面温度计算外表面的第一总热流；
S4：基于假定内表面温度计算内表面的第一对流换热系数，并基于第一对流换热系数计算内表面的第一对流换热量，其中，蒙皮内表面计算区域的非驻点区域的第一对流换热系数的计算公式为：
其中，hs,c表示第一对流换热系数，s表示沿弦向位置，ReΔ表示边界层雷诺数，ρa表示空气密度，ue表示空气附面层外边界处的气流速度，cp,a表示空气定压比热容，χ表示修正系数(用于修正蒙皮的表面曲率对计算结果的影响)；
蒙皮的计算区域的驻点区域的对流换热系数计算公式为：
其中，hstag表示驻点处的对流换热系数，Nustag表示驻点处的努赛尔数，Cm为常数，可以取为1.56，Re∞为空气来流雷诺数，ue为附面层外边界处的气流速度，u∞为来流空气速度，s为沿弦向位置，c为特征弦长，ka为空气的热传热率。
在一些实施例中，ue可以采用以下方法计算：
在求解对流换热系数时，需要求解附面层外气流场参数，定义压力系数表达式为：
其中，
表示压力系数，p表示空气压力，p∞表示来流的环境压力，u∞为来流的速度，ρa表示空气密度。
对压力系数进行可压缩的修正，修正后的可压缩压力系数可以表示为：
其中，
表示可压缩压力系数，Ma∞表示总马赫数。
局部压力利用可压缩压力系数进行修正，得到：
其中，κa为空气导热系数，而来流空气的总马赫数Ma∞、总温Tt和总压pt都可以由空气的熵方程获得：
附面层外边界处的马赫数Mae、温度Te和速度ue可以表示为：
其中，Ra为空气气体常数，κa为空气导热系数，T∞为来流温度。
S5：判断第一总热流是否等于导热量以及第一对流换热量是否等于导热量；当判断出第一总热流等于导热量，且第一对流换热量等于导热量时，执行S6：输出此时的假定外表面温度和假定内表面温度，否则，执行S7：对假定外表面温度和假定内表面温度进行修正，基于修正后假定外表面温度和假定内表面温度返回到步骤S1。往复循环，继续基于修正后的假定外表面温度和假定内表面温度进行后续计算，直至计算出的第一对流换热量等于导热量，且第一总热流等于导热量时，输出此时的蒙皮的假定外表面温度和假定内表面温度，也即是计算预测得到的蒙皮外表面温度和内表面温度。
在一些实施例中，S4包括步骤：基于假定内表面温度计算蒙皮的内流场，并进一步地基于计算所得的内流场来计算第一对流换热系数。其中，蒙皮的内流场与蒙皮的几何结构、防冰腔内的几何结构以及防冰腔的工作参数(如射流孔的喷射速率或喷射量，或所喷射的热流的温度等)密切相关。因此，在现有技术中，为了提高预测方案对复杂防冰腔结构的适应性(即针对复杂防冰腔结构构建相适应的表面温度预测方案，以提高预测精度)，通常会选择对内流场的计算过程、计算方式进行优化。而与现有技术不同的是，本发明是基于预测方案中的另一个关键参数计算过程(即第一对流换热系数计算过程)对预测方案进行优化，也即是从内流场计算的结果上对内流场计算过程中可能产生的误差进行修正。
本实施例中，为了提高预测方法的预测效率，本发明不对非驻点区域进行进一步的区分，而将非驻点区域直接作为湍流区域进行计算，从而避免或减少了部分复杂的计算迭代过程；另外，本发明对预测方案中第一对流换热系数的计算进行了修正，考虑到了机翼表面的曲面外形对空气流场的影响(相反地，现有技术中则通常会选用平板模型来进行计算)，从而保证了预测结果的准确性。
具体地，在一些实施例中，选取蒙皮内表面任意一个节点，节点两侧线段的夹角为θ，当表面为凸表面时，θ＜180；当表面为凹表面时，θ＞180，相应地，考虑几何外形的对流换热系数的修正系数的计算公式为：
其中χ表示修正系数，θ为表面节点的两侧线段夹角。
本实施例中对预测方案中第一对流换热系数的计算进行了修正，考虑到了蒙皮的曲面结构对蒙皮内流场产生的影响，进而保证了第一对流换热系数以及最终计算结果的精度，并且计算数据量较少。换句话说，即上述技术方案在未牺牲计算效率的同时还保证最终预测结果的精度。
进一步地，在一些实施例中，基于上述的对流换热系数计算方法，由于仅计算湍流区域的对流换热系数，边界层雷诺数的计算公式也进一步地简化为：
其中，μa表示空气的动力粘度，ρa表示空气密度，ΔT1表示空气附面层外边界气流温度与液态水表面温度之差的绝对值。
本实施例中，基于上述的第一对流换热系数计算方案，本发明还提供了新的边界层雷诺数计算方案，进一步地简化了表面温度预测方案，从而减少了预测过程中的数据计算量，以提高计算效率。
进一步地，在一些实施例中，S3还包括步骤：
基于假定外表面温度计算外部流场，基于外部流场确定蒙皮外表面的第二对流换热量；
基于第二对流换热量计算得到外表面的各项热流，并基于各项热流计算第一总热流。
例如，在一些实施例中，外部流场通过质量守恒方程、动量守恒方程和能量守恒方程进行求解，具体地，将上述三个守恒方程联立，可以得到三维无量纲积分形式的N-S方程如下：
其中，
为守恒变量,
为对流通量,
为粘性通量,
为控制体内的源项,Ω—为控制体体积,
为面上的法向矢量。其中，控制体指的流场中的某一确定的空间区域，这个空间区域的边界称为控制面，如图16所示，图16所示的区域被曲线划分为了多个区域，每个区域就是一个控制体。
进一步地，在一些实施例中，基于第二对流换热量计算得到外表面的各项热流，并基于各项热流计算第一总热流的步骤，包括：
基于第二对流换热量，采用欧拉法计算水滴撞击特性，得到水滴撞击特性结果；
基于水滴撞击特性结果，采用基于改进Messinger模型部件表面热载荷的计算方法，计算蒙皮外表面的各项热流；
基于各项热流计算第一总热流。
在一些实施例中，水滴撞击特性计算方法包括：拉格朗日方法和欧拉方法。
优选地，本实施例中选用欧拉方法计算水滴撞击特性。
为了进一步地说明本方法所采用的技术方案以及能够实现的技术效果，下文选取了某型飞机全尺寸翼段作为计算模型，对计算具体步骤加以说明，并对计算结果进行验证：
具体地，选取了某型飞机全尺寸翼段作为计算模型，模型为等直翼段，其弦长和展长均为1.4m，翼段几何外形如图5a所示。采用ICEM软件对翼段进行结构网格划分，网格总数为57.8万，网格正交性大于0.2，网格角度大于18°，图5b中展示了计算域整体和局部网格。
根据该型飞机的飞行包线和遭遇的云雾包线，选取了如表1所示的计算状态，利用本文建立的三维翼面防冰热载荷计算方法进行了翼段表面防冰热载荷的计算，其中防冰表面温度取15℃(即假定的蒙皮外表面温度取15°)。
表1为防冰热载荷计算状态
水滴收集系数是防冰热载荷的重要影响因素，参见图6a，图6a中的机翼表面中机翼前缘的驻点区域(即虚线区域)的水收集系数最大，且机翼的水收集系数从机翼前缘区域朝向两个箭头所示方向逐渐减小。从图6a和图6b中可以看出，水滴主要撞击在机翼前缘附近，最大水收集系数约为0.65。
图7a、图7b分别给出了三维翼段表面防冰热载荷分布云图和展向50％剖面的防冰热载荷分布曲线。参见图7a和图7b，图7b中的a点对应为图7a中的驻点(如图7a中的虚线所示区域)的热载荷值，从图中可以看出，热载荷主要分布在机翼前缘部分，最大热载荷达到了34000W/m2；热载荷曲线呈现中间小，随后向两边增大，最后为0的趋势，这是由于驻点位置静温较高，出现了较小的防冰热载荷需求量，而沿着弦向水收集量逐渐减少其热载荷也不断减小最后直至为0。对表面热载荷进行累积计算，得到计算状态总热载荷为8103.7W。
飞机防冰腔热效率通常在60％左右，在防冰热气引气量估算中本实施例取60％，热效率η可以表示为：
其中Tin为热气防冰系统引气温度，Tout为防冰系统出口温度，Ts为环境温度。由防冰腔热效率60％可以得到下式：
Tin-Tout＝(Tin-Ts)×η (9)
根据热载荷的定义，蒙皮表面总热载荷应等于热气流量的能量损失，可以得到如下计算公式：
Qr＝Cp×Qm×(Tin-Tout) (10)
其中Qr为总热载荷，Cp为定压比热容，Qm为引气流量。
由于该飞机发动机引气温度大约在246℃，结合表1，由式(8)-(10)可以计算出热气引气流量大约为0.052kg/s。
热气防冰系统工作的物理过程是一个外部流动传热(外部流动传热分析包括外部流场的计算、外部水滴运动场的计算、蒙皮外表面热流分析等)、内部流动传热以及固体导热的复杂耦合过程。在外部流动传热过程中，机翼蒙皮外表面是暴露在结冰气象条件环境中，空气中的过冷水滴直接撞击到机翼表面，将发生复杂的质量和能量传递，使得蒙皮外表面温度发生变化；在内部流动传热过程中，机翼蒙皮内表面与防冰腔内的热空气发生对流换热，使得蒙皮内表面温度发生变化；在固体导热过程中，机翼蒙皮本身会发生热传导，影响了蒙皮内外表面温度分布。因此，在热气防冰系统性能仿真实际上是要同时考虑外部热问题、内部热问题以及固体导热问题，进而获取机翼蒙皮内外表面的温度分布。
内部流动传热问题是笛形管射流热气在防冰腔内的流动并与蒙皮内表面强迫对流换热的问题。采用有限体积法通过求解质量守恒方程、动量守恒方程和能量守恒方程，可获得内部流场、温度分布以及蒙皮内表面温度分布情况。
防冰腔内部笛形管射流通常为超音速射流，对于可压缩流体，密度随压力和速度的变化而变化，其质量守恒方程为：
式中：
为热空气流速，ρa为热空气密度。
动量守恒方程为：
式中σij为剪切应力；
式中，p为静压，δjk、δik、δij分别表示沿ij,jk和ik方向的剪切应力，υi、υj、υk
分别表示热空气沿i,j,k方向的速度，τij表示沿ij方向的应力。
对于粘性流体，其计算公式为：
μ为空气的动力粘度，T∞和μ∞表示参考变量，分别为273.15K和17.9×10-6Pa·s。
基于能量守恒原理，系统的总焓保持不变，其数学表达式为：
式中，E为总能，Va为热空气流速，H为总焓，T为温度，υi为沿i方向的热空气速度，τij为沿ij方向的应力，κ为热传导率，其计算公式为：
值得注意的是，对于三维问题，控制方程组共有8个方程，但却有9个未知数(ρa,p,T,H,vx,vy,vz,μ,k)，为封闭方程组，添加理想气体状态方程：
pa＝ρaRTa (17)
其中，pa为压强，R为摩尔气体常数，Ta为温度。
蒙皮表面温度是外部气液两相流场，内部热空气流场和蒙皮自身导热综合作用的结果。其中，固体非稳态导热计算方法如下：
固体的热传导可采用如下的偏微分方程表示：
这里K表示导热系数，Cps表示固体定压比热，ρs表示固体密度，ΔT表示相邻两个迭代步之间的温度差。
求解控制方程时，从初始温度和初始焓开始迭代，新的温度(当前迭代步骤的温度)采用公式(2)计算，对应的焓的计算公式参见公式(2b)。
达到稳态，即计算结果收敛时，ΔT趋于极小值。
在上述内外流场以及固体导热的计算方法的基础上，将三者进行耦合，计算蒙皮表面温度，其流程如图8所示。
具体计算步骤为：
(1)假定蒙皮外表面温度为T0，此温度作为输入参数，由设计人员给定，假定的蒙皮外表面温度越接近最终收敛结果，则蒙皮表面温度耦合计算过程进行的循环步数越少，收敛速度越快；
(2)计算的部件外部空气流场结果初始化外流场；
(3)基于步骤(2)空气流场的计算结果，采用欧拉法计算水滴撞击特性；
(4)基于步骤(3)计算得到的水滴撞击特性结果，采用基于改进Messinger模型部件表面防冰热载荷的计算方法，计算蒙皮外表面的各项热流，从而确定蒙皮外表面需要吸收的热量(即第一总热流)；
(5)假定蒙皮内表面温度，此温度也是由设计人员给定，计算的热气防冰腔内部流场结果初始化内流场；
(6)根据蒙皮的内外表面温度(即假定内表面温度和假定外表面温度)，计算蒙皮的导热量；
(7)比较步骤(6)中计算出的导热量与步骤(5)中计算出的对流换热量(即第一对流换热量)，若不相等，则修正对流换热量及蒙皮内表面温度(即修正蒙皮的假定内表面温度)；
(8)比较步骤(6)计算出的导热量与步骤(4)计算出的总换热量(即第一总热流)，若不相等，则修正总换热量及蒙皮外表面温度(即修正蒙皮的假定外表面温度)；
(9)运用新的蒙皮外表面温度T1，重新计算外流场，获得外表面对流换热系数(第二对流换热系数)及对流换热量(第二对流换热量)，计算水滴撞击特性和蒙皮外表面各项热流，确定通过蒙皮外表面的总热流(第一总热流)；
(10)运用新的蒙皮内表面温度(即修正后的假定外表面温度和假定内表面温度)，计算内流场，获得内表面对流换热系数(第一对流换热系数)及对流换热量(第一对流换热量)；
(11)根据新的蒙皮内外表面温度，基于上述公式(18)计算蒙皮导热量；
(12)比较步骤(11)中计算出的导热量与步骤(10)中计算出的对流换热量(第一对流换热量)，若不相等，则修正对流换热量及蒙皮内表面温度；
(13)比较步骤(11)计算出的导热量与步骤(9)计算出的总换热量，若不相等，则修正总换热量及蒙皮外表面温度；
(14)重复循环步骤(9)～(13)，直至蒙皮外部换热量，蒙皮导热量及蒙皮内表面换热量达到平衡，此时蒙皮的表面温度分布即为所要求的温度分布。
进行耦合计算时，外部流场计算模型、内部流场计算模型以及固壁计算模型的几何位置的对应关系非常重要，具体方案为：固壁外表面与外部流场计算模型蒙皮内表面重合；固壁内表面与内部流场计算模型外边界重合。这主要是因为内外及固壁耦合计算时，数据的传输是通过节点数据交换获得，因此三者在几何位置上必须对应。
在一些实施例中，步骤(9)中的第二对流换热系数也可以选用公式(3a)、(5)以及(6a)进行计算。
下文进一步地对本实施例中的仿真计算及结构进行说明：
计算模型是在上述的全尺寸翼段数模基础上，采用图2中的方案(d)的结构形式，设计了热气防冰系统，如图9a所示，其中笛形管管径为17mm，射流孔直径为1.6mm，射流孔夹角为30°且等距交错排布。计算网格采用C形结构网格，并对前缘防护区域、防冰腔前缘、笛形管等区域进行了进行了加密，网格总数量为800万左右，如图9b所示。
外部流传传热计算工况以及蒙皮材料的相关参数如表2、表3所示。基于上述步骤中计算得到的引气量，以及该飞机发动机引气的压力，确定引气参数条件如表4所示。
表2
表3
表4
图10a为笛形管防冰腔静压分布图，图10a所示的笛形管防冰腔中的相邻喷孔间压差较大且呈周期分布(图10a中的喷孔对应区域(即图10a中的白色圆点区域)的表面压力明显地高于其他区域的表面压力)，这也表明翼段前缘换热是不均匀的。图10b为翼段中间位置喷孔截面上的Ma数分布，图中表明了腔内的流动主要由喷口处的高速射流和腔内其它位置的低速绕流组成，热气由喷口高速喷出到达物面后沿防冰腔内表面流动，经过上下端的扁平通道后进入大气环境(图10b中Ma数沿虚线所在区域分别沿两个箭头所示方向逐渐减小)。图10c中流线分布充分展现了防冰腔内部的流动的三维效应，说明了防冰腔内流场并不适合采用二维的方式进行仿真计算。
耦合计算获得稳态的残留水膜分布(图10d中的(a)、(b)分别展示了防护区上、下表面水膜厚度分布，且图10d中的蒙皮最前缘的区域的溢流液态水膜最薄，且上下表面的水膜厚度在横向上接近周期性变化)、可能的溢流冰(图10e中的(a)、(b)分别展示了防护区上、下表面积冰增长率分布)、以及蒙皮表面温度(图10f，图10f中的蒙皮前缘对应于喷孔的区域表面温度最高，蒙皮前缘上下表面的表面温度沿逐渐远离喷孔处逐渐降低)，由图可知，在防护去上下表面存在溢流液态水膜，水膜流出防护区后会行程溢流结冰。
为了验证计算结果的准确性，进一步地，通过风动试验评估热气防冰性能，并与上述计算得到的蒙皮表面温度(具体地，计算蒙皮的外表面温度)对比。
具体地，试验在中国空气动力研究与发展中心3米×2米结冰风洞主试验开展。
试验模型与本实施例中选用的计算仿真模型相同，采用不锈钢加工制成。如图11a所示，试验模型包括前缘部段11f、后缘部段11g。试验模型通过上连接件11a和下连接件11e分别与主试验段上、下转盘固连，分别采用定位销+螺钉的连接方式固定连接。模型迎角变化由上下转盘同步转动实现。模型迎角控制精度为±0.05°。模型在风洞试验段内的安装情况如图11b所示。模型前缘部段内设有防冰腔和笛形管。防冰腔同样采用图2中的方案(d)的形式，蒙皮厚度约为1.7mm。笛形管长为1.4m，管外径为17mm，管壁厚为2mm；射流孔为圆形，孔直径为1.6mm，射流孔夹角为30°且等距交错排布。
试验共选取了三个横截面作为温度测量面，由下到上分别位于z＝220mm、z＝660mm和z＝1100mm处(坐标系如图11a所示)。温度测量采用Omega公司TT-K-36/SLE型热电偶。表5给出了各测温截面坐标以及传感器详细编号，表中“11b”、“11c”、“11d”截面分别表示由下到上的截面1#、截面2#、截面3#等三个测温截面(如图11a所示)，在每个测温截面的上下表面按相同弦长的位置分别对称布置了12个传感器，同时在前缘驻点布置了一个传感器，每个测温截面共25个热电偶传感器，三个测温截面传感器编号方式相同。以截面1#为例，测温传感器位置分布及编号如图12所示，模型上表面处，以约为0.17c处E0点开始，沿前缘点方向依次为E1、E2逐渐到E11，每个传感器之间沿弦向距离约为0.014c，此处c为弦长；模型下表面编号方式相同。
表5
试验条件与上述步骤中的计算状态相同，具体试验工况见表3。防冰热气参数同样参考上述的计算条件，引气温度为246℃，引气流量为0.052kg/s，引气压力为0.29MPa。
表6
为了确保试验安全进行、试验数据可靠，试验严格按照相应的标准流程进行，具体步骤主要包括：
(1)模型检查和安装
试验准备前将试验件置于平台上检测，确保模型状态良好，配件齐全。翼型前后缘表面中心对称线清晰。模型进洞安装后，确保模型、连接附件、转盘之间固定连接可靠，安装到位。检测模型姿态角，初始姿态角误差应小于0.05°。
(2)防冰试验流程
a.风洞稳风速运行，同时启动运行高度模拟系统，调节控制压力和风速；
b.当风速和压力稳定后，再启动制冷系统，温度开始下降；
c.启动防冰供气系统，加热机翼前缘；
d.当气流温度达到试验条件，且模型表面温度稳定后，按预先标定的云雾模拟水压、气压组合开启喷雾系统；
e.监视防冰过程记录喷雾时间，待达到试验条件要求的防冰时间后，同时关闭喷雾系统和防冰供气系统、制冷系统、动力系统及高度模拟系统；
f.工作人员穿防寒服进入试验段对防冰结果进行检测，包括照相、测量等；
g.对模型表面进行除冰，除冰结束后可进行下一次试验。
(3)结果记录与测量
试验结束后，利用相机、卷尺对防冰效果进行拍照测量；利用相机对冰形正面、侧面、细节进行取照，记录结冰区域及外形特征；模型防冰过程由风洞摄像监控系统录制。
进一步地，对试验结果进行说明与分析。
图13a和图13b为防冰试验结果图，由图可以看出，整个翼段前缘防冰区域均未发现结冰现象，说明系统达到了防冰的要求。通过试验结果看出，在翼段后缘出现了部分溢流冰，这是因为前缘收集的水在防冰表面并未完全蒸发，流出防冰区域后逐渐冻结进而形成了流线状的结冰；此外，上连接件11a、下连接件11e前缘出现了明显的积冰，这是由于在上、下连接件前缘未设有防冰系统所致。
图14a为测温中截面前缘点防冰过程度随时间变化图，由图可知，随着喷雾开始，测温点的度迅速下降，最后稳定在5℃左右。图14b为防冰稳定后三截面温度分布图，由图可知，三个截面测温点温度均高于0℃，且截面由下到上对应位置的温度逐渐降低，这是由于供气气流沿笛形管流动也是由下向上，且笛形管未采用变径设计，气流沿程温度、压力均有所下降所致；此外，防冰表面前缘驻点位置温度较低，这可能是由于前缘水收集量较大和笛形管射流孔夹角过大所致。
表7给出了3个测温截面温度的统计结果。1#、2#、3#截面测量的最高温度值分别为14.46℃、12.08℃、13.42℃，平均值分别为9.8 8℃、5.74℃、2.53℃。从统计数据可℃。从统计数据可以看出，沿着模型由下向上，表面温度平均值都呈减小趋势。
表7
| 截面编号 | 最大值℃ | 最小值℃ | 平均值℃ | 标准差 |
| 1# | 14.46 | 3.27 | 9.88 | 3.14 |
| 2# | 12.08 | 2.06 | 5.74 | 2.69 |
| 3# | 13.42 | -1.43 | 2.53 | 4.10 |
选取1#截面温度分布的预测值与试验值进行了对比，如图15a所示，由图可知，模型上表面温度较高，预测值和试验值的最大温度值接近，温度变化趋势相似。图15a中温度分布的预测值与试验结所得的试验值虽然存在一定的偏差(误差)，但该偏差处于本领域的合理偏差范围(误差范围)，或者说，该误差在实际的工程应用上影响很小，可以忽略不计。如图15a所示，基于本发明预测方法得到的预测值大部分都与试验值相吻合(与试验值趋势一致)，该预测实例中所得到的约80％的预测值都与实验值相吻合，仅有少部分落在误差范围外(从附图15a可见，共计19个预测值，其中15个预测值与试验值相同或基本吻合)，符合实际工程应用过程中对预测结果精度的要求。
进一步地，上述偏差的产生主要原因有以下两点：1)试验中温度传感器采用开孔安装的方式安装在机翼表面，虽然对机翼表面做了平滑处理，但仍会小幅影响机翼前缘部分的流场，进而影响前缘表面的温度分布的测量；2)三维水滴撞击、水膜流动的实际过程非常复杂，开展计算需要对数值模型进行合理假设和简化，这也会与试验产生一定的偏差。此外，在防冰试验中观测到了明显的溢流冰，如图15b所示，这与计算获得的现象相吻合(如图10e所示，图10e的(a)、(b)中浅色区域的的溢流冰明显地多于深色区域的溢流冰，蒙皮前缘(或靠近蒙皮前缘)区域没有溢流冰生成，或者溢流冰生成量很小，而蒙皮表面除了靠近蒙皮前缘的区域均有一定量的溢流冰生成)，这进一步相互验证了本发明的预测方法的正确性。
本实施例中以某型飞机全尺寸翼段模型为研究对象，系统开展了热气防冰系统性能评估研究，完成了典型状态下该全尺寸翼段的防冰热载荷和热气引气参数计算。然后，基于上述实施例中的预测方法对该翼段模型进行了表面温度预测和防冰性能评估。
如前所述，蒙皮表面温度的预测和试验的最大温度值接近，温度变化趋势相似，表明本发明的预测方法得到的计算结果和试验结果相关性和一致性良好，验证了本发明建立的防冰系统性能评估方法的正确性。
如前所述，采用基于本发明建立的热载荷计算方法得到的热气引气参数(引气温度为246℃，引气流量为0.052kg/s，引气压力为0.29MPa)能够保证防冰区域无积冰，有效实现了防冰功能；
实施例二
参见图4b，本发明第二方面在于，基于上述实施例中的预测方法，还提供了一种用于预测复杂防冰腔所能维持的蒙皮表面温度的预测系统，包括：
温度假定模块，被配置为用于对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和假定内表面温度；
导热量计算模块，被配置为基于所述假定外表面温度和所述假定内表面温度计算所述蒙皮的导热量；
外表面总热流计算模块，被配置为基于所述假定外表面温度计算外表面的第一总热流；
内表面对流换热量计算模块，被配置为基于所述假定内表面温度计算内表面的第一对流换热系数，并基于所述第一对流换热系数计算内表面的第一对流换热量，其中，所述蒙皮的非驻点区域的所述第一对流换热系数的计算公式为：
第一判断模块，被配置为用于判断所述第一总热流是否等于所述导热量，以及所述第一对流换热量是否等于所述导热量；
温度输出模块，用于当所述第一判断模块判断出所述第一总热流等于所述导热量，以及所述第一对流换热量等于所述导热量时，输出此时的所述假定外表面温度和所述假定内表面温度；
温度修正模块，用于当所述第一判断模块判断出所述第一总热流不等于所述导热量，或者所述第一对流换热量不等于所述导热量时，对所述假定外表面温度和所述假定内表面温度进行修正，并将修正后的所述假定外表面温度和所述假定内表面温度发送到所述温度假定模块。
在一些实施例中，所述修正系数的计算公式为：
在一些实施例中，所述边界层雷诺数的计算公式为：
其中，μa表示空气的动力粘度，ρa表示空气密度，ΔT 1表示空气附面层外边界气流温度与液态水表面温度之差的绝对值。
在一些实施例中，所述外表面总热流计算模块包括：
外部流场计算单元，被配置为基于所述假定外表面温度计算外部流场，并基于所述外部流场计算蒙皮外表面的第二对流换热量；
外部热流计算单元，被配置为基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流。
在一些实施例中，所述外部热流计算单元包括：
第一计算子单元，被配置为基于所述第二对流换热量，采用欧拉法计算水滴撞击特性，得到水滴撞击特性结果；
第二计算子单元，被配置为基于所述水滴撞击特性结果，采用基于改进Messinger模型部件表面热载荷的计算方法，计算蒙皮外表面的所述各项热流；
第三计算子单元，被配置为基于所述各项热流计算所述第一总热流。
实施例三
本发明第三方面在于，基于上述实施例中的预测方法，还提供了一种用于预测复杂防冰腔能够维持的蒙皮表面温度的计算机程序产品，其中存储有计算机程序，计算机程序被处理器执行时，控制计算机程序产品所在设备执行上述实施例中任一方法的步骤。
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者装置不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者装置所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者装置中还存在另外的相同要素。
以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。
Claims (10)
1.一种复杂防冰腔所能维持的蒙皮表面温度预测方法，其特征在于，包括步骤：
S1：对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和假定内表面温度；
S2：基于所述假定外表面温度和所述假定内表面温度计算得到蒙皮的导热量；
S3：基于所述假定外表面温度计算外表面的第一总热流；
S4：基于所述假定内表面温度计算内表面的第一对流换热系数，并基于所述第一对流换热系数计算内表面的第一对流换热量，其中，所述蒙皮的非驻点区域的所述第一对流换热系数的计算公式为：
S5：判断所述第一总热流是否等于所述导热量以及所述第一对流换热量是否等于所述导热量；当判断出所述第一总热流等于所述导热量，且所述第一对流换热量等于所述导热量时，执行S6：输出此时的所述假定外表面温度和假定内表面温度，否则，执行S7：对所述假定外表面温度和假定内表面温度进行修正，基于修正后假定外表面温度和假定内表面温度返回到步骤S1。
2.根据权利要求1所述的方法，其特征在于，所述修正系数的计算公式为：
3.根据权利要求1所述的方法，其特征在于，所述边界层雷诺数的计算公式为：
4.根据权利要求1-3任一所述的方法，其特征在于，所述S3还包括步骤：
基于所述假定外表面温度计算外部流场，基于所述外部流场确定蒙皮外表面的第二对流换热量；
基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流。
5.根据权利要求4所述的预测方法，其特征在于，所述基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流的步骤，包括步骤：
基于所述第二对流换热量，采用欧拉法计算水滴撞击特性，得到水滴撞击特性结果；
基于所述水滴撞击特性结果，采用基于改进Messinger模型部件表面热载荷的计算方法计算蒙皮外表面的所述各项热流；
基于所述各项热流计算所述第一总热流。
6.一种用于预测复杂防冰腔所能维持的蒙皮表面温度的预测系统，其特征在于，包括：
温度假定模块，被配置为用于对蒙皮的外表面温度和内表面温度进行假定，得到假定外表面温度和假定内表面温度；
导热量计算模块，被配置为基于所述假定外表面温度和所述假定内表面温度计算所述蒙皮的导热量；
外表面总热流计算模块，被配置为基于所述假定外表面温度计算外表面的第一总热流；
内表面对流换热量计算模块，被配置为基于所述假定内表面温度计算内表面的第一对流换热系数，并基于所述第一对流换热系数计算内表面的第一对流换热量，其中，所述蒙皮的非驻点区域的所述第一对流换热系数的计算公式为：
第一判断模块，被配置为用于判断所述第一总热流是否等于所述导热量，以及所述第一对流换热量是否等于所述导热量；
温度输出模块，用于当所述第一判断模块判断出所述第一总热流等于所述导热量，以及所述第一对流换热量等于所述导热量时，输出此时的所述假定外表面温度和所述假定内表面温度；
温度修正模块，用于当所述第一判断模块判断出所述第一总热流不等于所述导热量，或者所述第一对流换热量不等于所述导热量时，对所述假定外表面温度和所述假定内表面温度进行修正，并将修正后的所述假定外表面温度和所述假定内表面温度发送到所述温度假定模块。
7.根据权利要求6所述的预测系统，其特征在于，所述修正系数的计算公式为：
所述边界层雷诺数的计算公式为：
8.根据权利要求6-7任一所述的预测系统，其特征在于，所述外表面总热流计算模块包括：
外部流场计算单元，被配置为基于所述假定外表面温度计算外部流场，并基于所述外部流场计算蒙皮外表面的第二对流换热量；
外部热流计算单元，被配置为基于所述第二对流换热量计算得到所述外表面的各项热流，并基于所述各项热流计算第一总热流。
9.根据权利要求7所述的预测系统，其特征在于，所述外部热流计算单元包括：
第一计算子单元，被配置为基于所述第二对流换热量，采用欧拉法计算水滴撞击特性，得到水滴撞击特性结果；
第二计算子单元，被配置为基于所述水滴撞击特性结果，采用基于改进Messinger模型部件表面热载荷的计算方法计算蒙皮外表面的所述各项热流；
第三计算子单元，被配置为基于所述各项热流计算所述第一总热流。
10.一种用于预测复杂防冰腔能够维持的蒙皮表面温度的计算机程序产品，其中存储有计算机程序，其特征在于，所述计算机程序被处理器执行时，控制所述计算机程序产品所在设备执行如权利要求1至5中任一所述方法的步骤。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2021116725491 | 2021-12-31 | ||
| CN202111672549 | 2021-12-31 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114940266A true CN114940266A (zh) | 2022-08-26 |
| CN114940266B CN114940266B (zh) | 2023-04-25 |
Family
ID=82905779
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210302199.8A Active CN114940266B (zh) | 2021-12-31 | 2022-03-25 | 一种复杂防冰腔所能维持的蒙皮表面温度预测方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114940266B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115524131A (zh) * | 2022-09-13 | 2022-12-27 | 中国航发沈阳发动机研究所 | 一种基于非结冰条件的整机防冰系统验证方法 |
| CN115817822A (zh) * | 2023-02-09 | 2023-03-21 | 中国空气动力研究与发展中心低速空气动力研究所 | 电加热防冰系统的热载荷分布设计方法 |
| CN116611266A (zh) * | 2023-07-19 | 2023-08-18 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种非稳态电加热防冰系统模拟方法 |
| CN116644258A (zh) * | 2023-07-19 | 2023-08-25 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种防冰热载荷计算方法及系统 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140257771A1 (en) * | 2011-11-30 | 2014-09-11 | Ming Lu | Numerical simulation method for aircrasft flight-icing |
| CN106650081A (zh) * | 2016-12-19 | 2017-05-10 | 北京航空航天大学 | 冲击‑热气膜复合式防冰结构设计方法 |
| US20180229850A1 (en) * | 2017-02-15 | 2018-08-16 | Pratt & Whitney Canada Corp. | Anti-icing system for gas turbine engine |
| CN109883660A (zh) * | 2017-12-01 | 2019-06-14 | 中国飞机强度研究所 | 一种热模拟试验控制方法 |
| CN112197932A (zh) * | 2020-12-07 | 2021-01-08 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种飞机防除冰试验模型表面温度修正方法及测量方法 |
| CN112989727A (zh) * | 2021-05-10 | 2021-06-18 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种防冰系统的壁面温度模拟方法 |
-
2022
- 2022-03-25 CN CN202210302199.8A patent/CN114940266B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140257771A1 (en) * | 2011-11-30 | 2014-09-11 | Ming Lu | Numerical simulation method for aircrasft flight-icing |
| CN106650081A (zh) * | 2016-12-19 | 2017-05-10 | 北京航空航天大学 | 冲击‑热气膜复合式防冰结构设计方法 |
| US20180229850A1 (en) * | 2017-02-15 | 2018-08-16 | Pratt & Whitney Canada Corp. | Anti-icing system for gas turbine engine |
| CN109883660A (zh) * | 2017-12-01 | 2019-06-14 | 中国飞机强度研究所 | 一种热模拟试验控制方法 |
| CN112197932A (zh) * | 2020-12-07 | 2021-01-08 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种飞机防除冰试验模型表面温度修正方法及测量方法 |
| CN112989727A (zh) * | 2021-05-10 | 2021-06-18 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种防冰系统的壁面温度模拟方法 |
Non-Patent Citations (1)
| Title |
|---|
| 杨秋明等: "基于流-固耦合传热的热气防冰系统干空气飞行蒙皮温度场计算研究" * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115524131A (zh) * | 2022-09-13 | 2022-12-27 | 中国航发沈阳发动机研究所 | 一种基于非结冰条件的整机防冰系统验证方法 |
| CN115524131B (zh) * | 2022-09-13 | 2024-03-19 | 中国航发沈阳发动机研究所 | 一种基于非结冰条件的整机防冰系统验证方法 |
| CN115817822A (zh) * | 2023-02-09 | 2023-03-21 | 中国空气动力研究与发展中心低速空气动力研究所 | 电加热防冰系统的热载荷分布设计方法 |
| CN115817822B (zh) * | 2023-02-09 | 2023-04-18 | 中国空气动力研究与发展中心低速空气动力研究所 | 电加热防冰系统的热载荷分布设计方法 |
| CN116611266A (zh) * | 2023-07-19 | 2023-08-18 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种非稳态电加热防冰系统模拟方法 |
| CN116644258A (zh) * | 2023-07-19 | 2023-08-25 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种防冰热载荷计算方法及系统 |
| CN116611266B (zh) * | 2023-07-19 | 2023-09-29 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种非稳态电加热防冰系统模拟方法 |
| CN116644258B (zh) * | 2023-07-19 | 2023-12-22 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种防冰热载荷计算方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114940266B (zh) | 2023-04-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114940266A (zh) | 一种复杂防冰腔所能维持的蒙皮表面温度预测方法及系统 | |
| Cao et al. | Numerical simulation of three-dimensional ice accretion on an aircraft wing | |
| Lynch et al. | The effect of the combustor-turbine slot and midpassage gap on vane endwall heat transfer | |
| CN109145388B (zh) | 航空发动机部件的热分析方法 | |
| Hasanzadeh et al. | Quasi-steady convergence of multistep Navier–Stokes icing simulations | |
| Radenac | Validation of a 3D ice accretion tool on swept wings of the SUNSET2 program | |
| Son et al. | Quantitative analysis of a two-dimensional ice accretion on airfoils | |
| Chen et al. | Hot-air anti-icing heat transfer and surface temperature modeling | |
| CN115292806A (zh) | 考虑周期性边界的三维热气防冰系统表面温度计算方法 | |
| CN114398844B (zh) | 基于连续水膜流动的稳态防冰仿真方法 | |
| McClain et al. | A model for ice accretion roughness evolution and spatial variations | |
| Wang et al. | A new wind turbine icing computational model based on Free Wake Lifting Line Model and Finite Area Method | |
| Saeed | Numerical simulation of surface heat transfer from an array of hot-air jets | |
| Liu et al. | Three-dimensional integrated thermodynamic simulation for wing anti-icing system | |
| Mei et al. | Implicit numerical simulation of transonic flow through turbine cascades on unstructured grids | |
| Levin et al. | Optimization of a low Reynolds number airfoil with flexible membrane | |
| Anderson et al. | Evaluation and validation of the messinger freezing fraction | |
| Agricola et al. | Turbine vane leading edge impingement cooling with a sweeping jet | |
| Özgen et al. | In-flight ice formation simulation on finite wings and air intakes | |
| Wong et al. | Computational investigation of a bleed air ice protection system | |
| El Ayoubi et al. | Aerothermal optimization and experimental verification for discrete turbine airfoil film cooling | |
| Holgate et al. | An experimental-numerical method for transient infrared measurement of film cooling effectiveness and heat transfer coefficient in a single test | |
| Ragab et al. | Heat transfer analysis of the surface of a nozzle guide vane in a transonic annular cascade | |
| Tufts et al. | Stability Equation Based Transition Prediction | |
| Ingram et al. | a superposition technique for multiple-row film-cooling for calculation of 2-d effectiveness and heat transfer coefficients |
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 |