飞机结冰是飞行安全的重大隐患之一，结冰造成的飞行事故约占总事故的12%，无人机由于结冰造成的事故约占总事故的 25%[1]。飞机结冰适航取证依赖于冰形预测。相比于成本高、周期长的冰风洞试验，结冰数值模拟是一种更方便的预测手段。结冰数值模拟主要分三个步骤：空气流场计算、水滴撞击特性计算和结冰热力学分析[7]。其中，最复杂的是第三步。从20世纪50—60年代至今，欧美各国针对结冰热力学分析做了大量工作。1953年，Messinger[14]提出了第一个用于航空的结冰热力学模型。2001年，Myers等[15-17]在Messinger模型的基础上进一步发展，提出了Myers模型。在结冰热力学分析中，对流换热是有显著影响的源项。对流换热现象十分复杂，前人对此做了大量研究。Frossling[18]计算了绕二维和轴对称的任意形状物体的层流边界层内的温度分布；Kayalar [20]从理论和试验两方面研究了湍流度对圆柱体换热的影响；Karman [22]将边界层细分为三个区域，并推导了传热系数与表面摩擦力系数之间近似关系式；Dipprey等[23]研究了粗糙表面的传热。以这些研究为基础，形成了目前常用的两种对流换热系数计算方法：边界层积分法和计算流体力学（CFD）直接计算对流换热系数的方法（简称为CFD直接法）[8]。边界层积分法产生于计算资源有限、难以求解N-S方程的年代，是一种半经验的方法。边界层积分法只需要物面粗糙度分布和边界层外速度分布即可求出对流换热系数。此方法在简单的二维表面有较高的准确度，并得到了广泛的应用和验证。随着计算能力的提升，求解RANS方程不再困难，出现了用于计算对流换热的CFD直接法[7]。此方法通过求解RANS方程得到壁面附近温度分布，用差分方法求壁面法向的温度梯度后，再用傅里叶导热定律求出对流换热系数[26]。CFD直接法可用于各种复杂形状和三维冰形计算，但由于梯度的精确计算易受网格密度、质量和湍流模型的影响，计算结果分散度较大[32-34]。目前，市面上已经出现了多套成熟的结冰模拟软件，多采用上述两种对流换热计算方法。如美国国家航空航天局（NASA）开发的LEWICE及LEWICE3D软件[2-3]，使用边界层积分法；加拿大NTI公司开发的FENSAPE-ICE软件[4-5]，使用CFD直接法；法国ONERA开发的Icing Code[6]，使用边界层积分法。AERO-ICE是清华大学自主开发的结冰数值模拟软件，基于RANS方程求解计算带冰流场，基于Messinger模型和Myers模型计算冰形，通过迭代实现二者耦合。AERO-ICE实现了完整的结冰流程模拟，所得的冰形计算结果得到了初步的验证。该软件原本利用CFD直接法计算对流换热系数。本文利用RANS方程求解计算的外流结果与边界层积分法相耦合，实现了对流换热系数的另一种计算方式。本文将对上述两种对流换热计算方法进行介绍，以一个对流换热算例和R203、R204两个结冰算例对比验证两种方法的精度等特性。1　对流换热系数的求解方法1.1　对流换热系数用热传导的傅里叶定律类比对流换热现象，可定义对流换热系数h。对流换热系数是当流体与固体表面之间的温度差为单位温度时，单位时间、单位面积上流体与固体间的对流换热量h=q˙wTw-Tt, ∞ （1）式中：q˙w为单位时间、单位面积的对流换热量；T为流体温度，Tw为壁面温度，Tt,∞为来流总温，如图1所示。10.12050/are20220310.F001图1对流换热示意图Fig.1Schematic diagram of convective heat transfer对流换热系数的地位与热导率相似，但本质不同。h不是流体本身的性质，而与当地流速、湍流度、温度、表面粗糙度等参数有关[9]。1.2　CFD直接法CFD直接法求解对流换热系数利用了壁面处对流换热量q˙w等于导热量q˙cond[26]的性质，用傅里叶导热定律求q˙wq˙w=q˙cond=k(∂T∂y)w （2）再代入式(1)可得h=k∂T∂ywTw-Tt,∞ （3）式中：k为空气热导率。本文使用专门针对前缘非流线体大分离改进的SPF k-v2-ω湍流模型求解RANS方程。对于积冰高度较高、流场中分离严重的算例，该模型能极大地提高计算准确性[24]。在该模型中，壁面粗糙度的影响通过设置ω的边界条件体现[24]。ω=0.09×τwμ⋅SR(ks+) （4）式中：τw为壁面切应力；μ为动力黏度；ks+为壁面粗糙度雷诺数。ks+=τwρksv （5）式中：ρ为流体密度；ks为等效沙粒粗糙高度；v为运动黏度。Wilcox基于试验结果推导了SR表达式[28]SR=(50/ks+)2ks+25100/ks+ks+≥25 （6）CFD求解RANS方程得到温度T的分布，用差分方法求壁面梯度后即可求得对流换热系数。CFD直接法在实践中还存在一些问题：（1）对同一算例，在其他设置基本相同的情况下，不同的湍流模型计算得到的对流换热系数h不同。Abdollahzadeh等[34]用不同湍流模型模拟了热平板流动中的对流换热。结果显示，对流换热系数计算结果确受湍流模型影响，这与本节结论一致。（2）对流换热系数计算结果有较强的网格敏感性[41-42]。对不同的流动条件，需要的网格也不同[40，43]，不同的网格甚至会使计算结果相差一个量级[39]。阎超等[39]指出，第一层网格高度对热流影响明显。研究者普遍用网格雷诺数来反映第一层网格高度Recell=ρ∞u∞Δyμ （7）式中：ρ∞为来流密度；u∞为来流速度；Δy为第一层网格高度。Zhang等[45]提出，对于圆柱的对流换热问题，当使用SST模型时，网格雷诺数应小于10。阎超等[39]的研究则显示，对于不同的算例存在不同的最佳网格雷诺数。钝双锥算例的最佳网格雷诺数约为200，钝椎算例的最佳网格雷诺数约为950。同时，Loomans[40]和阎超等[39]都强调，过大或过小的网格雷诺数都会使热流计算结果偏离真实值。对k-ε模型，过于精细的网格会高估换热量；而对零方程模型，更粗糙的网格会导致更大的换热量[40]。Olynick等[43]研究发现，对流换热计算对网格的敏感性还与差分格式有关，三阶格式的网格敏感性显著低于一阶格式。1.3　边界层积分法边界层积分法计算分为三个区域进行：驻点、层流区和湍流区。Smith等[11]近似给出了驻点对流换热系数的表达式h=kc[1.5629.89·Re∞·(d(ue/u0)d(s/c))s=sstag]0.5 （8）式中：c为翼型弦长；Re∞为来流雷诺数；ue为边界层外速度，可通过势流理论或求解N-S方程得到；s是距离驻点的弧长。转捩位置通过粗糙度雷诺数Rek判断[12]Rek=ukksν≥600 （9）式中：uk为y=ks处的空气流速；ν为运动黏度。依据Pohlhausen近似[13]，边界层内速度可表示为uue=2yδ-2(yδ)3+(yδ)4+16δ2νduedsyδ(1-yδ)3 （10）式中：δ为边界层厚度。代入y=ks可求得ukukue=2ksδ-2(ksδ)3+(ksδ)4+16δ2νduedsksδ(1-ksδ)3 （11）在层流区内，Smith等[13]发现对流换热系数h可由式（12）预测h=ka- 12ν(ue- b∫0xueb- 1dx)- 12 （12）对大部分工况下的空气，有a- 12=0.296，b=2.88，即h=0.296kν(ue- 2.88∫0xue1.88dx)- 12 （13）式（13）称为Smith & Spalding关系。湍流区的对流换热系数表达式为h=cf2ueρcpPrt+cf2Stk （14）cf2=0.168[ln(864δ2/ks)]2 （15）式中：cf为表面摩擦系数，与边界层动量厚度δ2和等效沙粒粗糙高度ks有关[9]。Prt为湍流普朗特数，通常取0.9；Stk为粗糙度斯坦顿数，由剪切速度uτ=uecf /2和等效沙粒粗糙高度ks决定[12]。Stk=1.16(uτksν)- 0.2 （16）边界层积分法不用求温度梯度，并且比较清楚地显示了物面粗糙度对对流换热系数的影响。但此计算过程需要沿物面积分，且默认空气依次流经驻点—层流区—转捩点—湍流区，因此难以用于计算复杂外形和三维机翼的对流换热。同时，边界层积分法中的转捩判据式(8)是基于经验，不一定适用于所有工况。但对于常见的较为规则的外形，该方法在计算效率上具有明显优势。2　AERO-ICE简介AERO-ICE是实验室近年自主开发的三维结冰模拟软件，有4个核心模块：网格自动生成、空气流场RANS计算、水滴场欧拉方法计算和结冰热力学分析。其中，结冰热力学分析模块同时具有Messinger和Myers两个热力学模型；将Messinger模型和Myers模型结合求解冰-水交界面能量方程。AERO-ICE的RANS计算采用专门针对前缘非流线体大分离计算改进的SPF k- v2- ω湍流模型。软件具有三维多块网格的计算能力，具有多重网格加速策略，并可以进行大规模并行计算。图2展示了AERO-ICE计算翼型结冰的流程。空气流场计算后，从流场数据中读入所需信息，使用选定的方法计算对流换热系数。10.12050/are20220310.F002图2AERO-ICE软件模拟翼型结冰流程图Fig.2Flowchart of AERO-ICE simulating airfoil icing常见的结冰软件（如LEWICE[2-3]）常采用势流模型计算边界层外速度ue。对于AERO-ICE，使用边界层积分法时，ue由RANS计算输出的壁面压力pw求得。在边界层内，沿物面法向压力不变[29]，边界层外压力pe等于壁面压力pe=pw （17）根据等熵方程可得边界层外马赫数Mae与边界层外压力pe的关系[30]Mae=2κ-1ptpeκ-1κ-1 （18）式中：κ为绝热指数，对于空气有κ=1.4；pt为来流总压。由边界层外马赫数可求得边界层外温度TeTe=Tt,∞1+κ-12Mae2- 1 （19）最后，边界层外速度为ue=MaeκRgTe （20）式中：Rg是气体常数，对于空气有Rg=287.053J/(kg•K)。ue=0处即为驻点。3　结果与分析3.1　对流换热算例本节通过计算NACA0012的对流换热算例，验证AERO-ICE中对流换热计算的准确性。算例工况见表1。10.12050/are20220310.T001表1对流换热算例工况[36]Table 1Working condition of convection heat transfer case[36]翼型c/mRe(Ts-T∞)/KNACA00120.5332.5×10630计算使用经过收敛性验证的网格（见3.2.2节），为方便与试验值对比，计算结果表示为基于弦长的Frossling数Frc=hc/kρuec/μ （21）图3为Frossling数分布，包含试验结果、AERO-ICE采用两种对流换热计算方法的结果和Samad等[35]的计算结果。AERO-ICE的边界层积分法结果（AERO-ICE IBL）与试验值吻合，这说明AERO-ICE的边界层积分法是可信的。在同一网格下，AERO-ICE的CFD直接法结果（AERO-ICE CFD）与试验值偏差较大，大于Samad等[35]的计算结果。10.12050/are20220310.F003图3Frossling数分布Fig.3Distribution of Frossling number3.2　结冰算例3.2.1　算例工况本节展示AERO-ICE的冰形计算结果和对流换热系数计算结果，包括R203和R204两个算例。算例的结冰工况见表2。10.12050/are20220310.T002表2验证算例的结冰工况[37]Table 2Working conditions of verified cases[37]算例R203R204翼型GLC305GLC305时间/min422.5MVD /μm2020LWC/（g/m3）0.540.54速度 /(m/s)9090温度 /K268.05268.053.2.2　网格收敛性验证为验证流场计算的网格收敛性，使用切向、法向节点数分别为185×69（粗网格）、361×137（中网格）、457×173（密网格）的三套网格模拟R203算例。为消除不同冰形的影响，只比较第一个时间步，即干净翼型上的计算结果。压力分布结果如图4所示，中网格和密网格的压力分布非常接近。CFD直接法求得的对流换热系数结果如图5所示，中网格与密网格的结果差别明显小于粗网格和中网格。但在驻点附近，中网格与密网格的结果差异仍大于边界层积分法，这会对冰形预测产生不利影响。边界层积分法求得的对流换热系数结果如图6所示，从中网格到密网格，对流换热系数改变量很小，已经达到收敛。本文的对流换热算，和结冰算例计算均采用节点数为361×137的中网格（见图7）。10.12050/are20220310.F004图4压力分布网格收敛性验证结果Fig.4Results of grid convergence verification of pressure coefficient10.12050/are20220310.F005图5CFD直接法的对流换热系数网格收敛性验证结果Fig.5Results of grid convergence verification of convectionheat transfer coefficient by CFD direct method10.12050/are20220310.F006图6边界层积分法的对流换热系数网格收敛性验证结果Fig.6Results of grid convergence verification of convectionheat transfer coefficient by Boundary-Layers Integral method10.12050/are20220310.F007图7后续计算所用网格（节点数为361×137）Fig.7Grid for subsequent calculation(The number of nodes is 361×137)3.2.3　结果对比分析图8和图10是冰形计算结果，图中包含试验测得冰形、LEWICE计算结果和AERO-ICE采用两种对流换热计算方法的结果，图 8还包括Cao等[38]的计算结果，冰角的定义采用2012年美国汽车工程师协会（SAE）制定的标准[31]。图9和图11是对流换热系数h分布，为消除冰形不同对h的影响，只展示第一步结果的对比，即在干净翼型GLC305上的对流换热系数。图中包含LEWICE计算结果和AERO-ICE采用两种对流换热计算方法的结果。10.12050/are20220310.F008图8R203算例冰形Fig.8Ice shape in case R20310.12050/are20220310.F009图9R203算例第一步结冰中的对流换热系数Fig.9Convective heat transfer coefficient in the first step oficing in case R20310.12050/are20220310.F010图10R204算例冰形Fig.10Ice shape in case R20410.12050/are20220310.F011图11R204算例第一步结冰中的对流换热系数Fig.11Convective heat transfer coefficient in the first step oficing in case R204对两个算例的冰形计算，AERO-ICE边界层积分法结果与LEWICE结果接近，均与试验冰形符合较好；R203算例中，AERO-ICE边界层积分法结果与Cao等[38]的结果基本吻合。这证明了本文实现的边界层积分法是准确可信的。在R204算例中，计算LEWICE结果与试验冰形偏差较大，而AERO-ICE边界层积分法结果则与试验冰形非常接近。分析表明，算例R204积冰高度较高，流场中分离严重（见图12）。此时LEWICE用势流理论计算难以获得准确流场，而AERO-ICE的流场采用SPF k-v2-ω模型求解RANS方程，显著提高了对分离流动计算的准确性。图13展示了AERO-ICE和LEWICE求出的压力分布，结果存在较大的偏差。10.12050/are20220310.F012图12使用SPF k- v2- ω 模型计算的R204算例流场Fig.12Flow field in case R204 calculated using SPF k- v2- ωmodel10.12050/are20220310.F013图13R204算例的压力分布对比Fig.13Comparison of pressure coefficient in case R204对比图8和图10中的AERO-ICE分别采用边界层积分法和CFD直接法的结果，可以发现，CFD直接法求得的冰角偏小。分析两个算例的冻结系数（见图14）和溢流水通量（见图15）分布可以发现，CFD直接法算得的对流换热系数h偏小，导致冻结系数f偏小，更多水向后溢流，使冰角偏小。R203和R204算例的工况完全相同，只是积冰时间不同。当采用相同的计算时间间隔时，第一步的冻结系数和溢流水通量完全相同，因此两个算例在图14和图15中完全重合。10.12050/are20220310.F014图14R203和R204算例第一步结冰中的冻结系数Fig.14Freezing fraction in the first step of icing in case R203and R20410.12050/are20220310.F015图15R20和R204算例第一步结冰中的溢流水通量Fig.15Overflow flux in the first step of icing in case R203 and R204根据Kim等[25]的试验结果，积冰靠近翼型前缘时（x/c3.5%），冰角越小，对升力系数的影响越大。即上述两个算例中CFD直接法的计算结果高估了结冰对升力系数的影响，会使容冰设计更加保守。在结冰过程中，随着积冰逐渐破坏了翼型的流线体外形，CFD直接法出现了振荡现象。图16是R204算例中，AERO-ICE分别用两种方法得到的第1080s的对流换热系数结果。从图中可以看到， CFD直接法结果有明显的振荡，而边界层积分法结果未发生震荡。Hanson等[46-47]的计算结果也存在相同的现象。10.12050/are20220310.F016图16R204算例结冰过程中1080s的对流换热系数Fig.16Convective heat transfer coefficient in 1080s of icing in case R204Henry等[48]测量了带冰翼型的对流换热系数，在他测量的4个冰形中，无论是明冰还是霜冰，都未发现类似的振荡现象，这表明这种振荡可能不是真实存在的。4　结论数值模拟是飞机结冰适航取证的重要手段，对流换热计算是结冰数值模拟的关键部分。本文基于自行发展的结冰模拟软件AERO-ICE，验证和对比了对流换热计算的CFD直接法和边界层积分法。具体结论如下：（1）在基于RANS求解流场的AERO-ICE软件中，利用边界层积分法实现对流换热系数计算。与CFD直接法计算对流换热系数相比，此方法对湍流模型和网格质量依赖性较低，并可较为简单地引入壁面粗糙度的影响。（2）边界层积分法能得到比较准确地对流换热系数分布，从而使得到的冰形与试验结果高度吻合。（3）对流换热系数对冰形影响强烈。对流换热系数偏小会导致算出的冰角偏小，高估结冰对升力系数的不利影响。