空中结冰极易影响飞机的飞行和操作性能，对飞行安全构成严重威胁。在对1994年Roselawn空难的调查中，人们发现过冷大水滴（supercooled large droplet，SLD）环境的危险性[1] 。2014—2015年，美国联邦航空局（FAA）和欧洲航空安全局（EASA）相继发布了SLD环境结冰的适航条款。我国C919大型客机将成为全球首个面对该条款审查的机型。SLD结冰的数值和试验验证技术至今仍不成熟，严重影响了我国民机安全设计与国际适航取证工作。试飞数据显示，SLD环境相较于小水滴环境（CCAR-25部附录C）的不同是含有超过100μm的大粒径水滴。FAA根据大量结冰环境试飞研究制定了SLD的环境标准14 CFR-25部附录O [2] ，包括冻毛毛雨和冻雨的4种SLD环境。欧洲EASA也采用了该环境标准。其中，水滴粒径分布均具有“双峰分布”特征，即水滴分别集中在大、小粒径两个区域中。但到目前为止，国际上还未能实现完整的SLD结冰环境试验模拟能力。冰风洞试验是飞机结冰模拟分析的主要手段之一，也是结冰数值模拟技术的最终校准对象。经过长期研究，美国国家航空航天局（NASA）Glenn中心冰风洞[3] 、加拿大国家研究委员会（NRCC）结冰风洞[5] 、意大利航天研究中心（CIRA）结冰风洞[6] 通过降低喷口气压及多种喷头组合的方式均实现了附录O 中的冻毛毛雨结冰环境，但粒径分布准确模拟还较为困难。在中国，空气动力发展与研究中心结冰风洞具备水滴平均粒径（MVD）10～300μm的模拟能力[7] ，但尚未看到其SLD试验能力的报道。据报道，航空工业气动院结冰风洞[8] 已经实现部分冻毛毛雨结冰条件。由于冰风洞中大水滴沉降[9] 和冷却慢[10-11]的限制，现有国内外冰风洞均只能实现冻毛毛雨环境（最大粒径500μm）的模拟。在SLD结冰数值计算方面，研究主要关注大水滴动力学效应的建模和计算，包括变形/破碎[12] 、飞溅[13] ，以及相应的数值算法[14] ，以准确模拟SLD的水滴收集率。Zhang等[15] 提出SLD撞击对非稳态传热的影响模型。研究表明，SLD碰撞回缩和结冰存在耦合[16-17]，使得SLD结冰和温度存在非线性关系[18] 。相关模型尚未应用到SLD数值模拟中。尽管SLD结冰模型和数值方法有了许多进展，但另一方面，SLD复杂的水滴粒径分布在数值模拟中考虑的仍然较少。目前，主流结冰软件通常以中值体积直径（median volume diameter，MVD）来代表水滴粒径分布，未具备复杂粒径分布的模拟能力。早期，SLD冰风洞研究已发现粒径分布对冰形有影响[19] ，这也成为其提高云雾粒径分布模拟精度的依据。朱程香等[20] 通过多尺度分布水滴环境的撞击特性和冰形研究认为，粒径分布形式对冰形有一定影响。显然，考虑SLD粒径分布对提升模拟精度十分有必要。对于水滴多尺度分布的模拟，现有研究主要关注粒径分布函数，如Shah等[19]研究了用Γ分布函数拟合粒径分布 ；Kollar等[21] 和国内符澄等[22] 均采用Rosin-Rammler分布函数对粒径的双峰分布进行拟合，但未看到SLD结冰计算中应用。精确模拟SLD粒径分布的难点在于数值计算代价过大。在结冰数值计算中，为更好地模拟SLD动力学效应，大多使用拉格朗日法[23] 计算大水滴轨迹和收集率，因其可以细致描述水滴流场信息。但其缺点在于需要计算大量粒子以获得较好的空间分辨率。近年来，欧拉两相流法[24] 因其计算效率优势被越来越多地用于SLD的运动模拟[25-26]。考虑到SLD环境的粒径分布特征，其中既包含稠密的小水滴群，也包含稀疏、随机的大水滴。单纯用一种方式难以准确模拟其分布特性。可以结合以上方法，配合粒径分布函数抽样的方式实现更真实的SLD分布特性。本文提出了一种可模拟复杂云雾粒径分布的欧拉-拉格朗日抽样混合算法，并基于开源CFD软件包OpenFOAM开发了SLD结冰模拟算法，实现了SLD“双峰分布”的高效准确数值模拟。通过相关文献和气动院的结冰风洞试验对比验证了算法，进行了粒径分布影响冰形规律的分析。1　SLD结冰数值计算方法1.1　数值计算框架本文结冰数值计算方法基于开源CFD工具包OpenFOAM平台开发，耦合了网格自适应生成、空气-水滴流场求解、表面水膜流动传热和结冰计算功能。计算将每个结冰时间步内的翼型绕流近似为定常流动，通过SIMPLE算法[15]计算。计算网格为自动生成的混合笛卡儿-贴体网格[27] ，如图1所示。每个结冰时间步后均根据当前冰形重新生成网格并计算空气流场。湍流模型采用SST模型[28] ，整个计算流程可MPI并行。10.12050/are20220308.F001图1结冰计算的笛卡儿-贴体网格Fig.1Cartesian body-fitted grids applied in icing calculation1.2　SLD粒径分布的抽样混合与收集算法根据适航标准14CFR-25部附录O，SLD的滴谱具有明显“双峰分布”特性[29] ，可以50μm粒径为界分为大/小两个水滴群[7] 。本文分别对两个水滴群进行抽样和轨迹计算。步骤如下。（1）将水滴分布离散为多个不同粒径的组分，将其水含量（LWC）以最小二乘法进行Rosin-Rammler分布函数拟合。其累计体积频度分布形式为CDF(d)=1-e- (d/d¯)n (1)式中：d¯为颗粒的体积中位数粒径；n为分布参数，用来描述颗粒群分布的集中度。这里使用两个Rosin-Rammler分布函数叠加表现SLD双峰分布[22] 。不同粒径水滴组分的体积频度函数为PDF'(d)=ωe- (d/d1¯)n1n1dn1-1d1n1+(1- ω)e- (d/d2¯)n2n2dn2-1d2n2 （2）式中：ω为颗粒群的质量分数；d¯1和d¯2分别为大小颗粒群的平均粒径；n1和n2分别为大小颗粒群的分布参数。这5个参数将通过将规定的SLD粒径分布代入式（2）中，以最小二乘法拟合得到。首先根据输入的小水滴的离散分布计算出相应的CDF'SD，SLD的粒径分布需要满足CDF'SD≤CDF'LD≤1，采用随机数的方式计算出SLD的每个抽样粒径对应的CDF'LD为CDF'LD=1-(1-CDF'SD)*random(0,1) （3）式中：random(0，1)为[0，1]之间的一个随机数。再通过式（2）利用二分法可求出抽样粒径的值。（2）对50μm以下的水滴组分，直接根据拟合的离散粒径di和LWC通过欧拉-欧拉方法求解流场中水滴的连续和动量方程[19]，得出水滴容积αi和水滴速度Udi，并得到在翼面上的边界值αi0和Udi0。（3）对50μm以上的水滴组分，对拟合的分布函数进行粒径和LWC的随机抽样，通过拉格朗日法计算水滴的运动轨迹，得到水滴收集量Mcollected。（4）将最终水滴收集结果进行叠加Mcollected=Mcollected+∑iαi0×(Udi0• S)S （4）式中：αi0×(Udi0•S)S为欧拉法算出的翼面的质量源项，S为翼面面元矢量。针对14CFR-25部附录O环境4条分布曲线的离散拟合如图2所示，离散间距为20μm，拟合参数见表1。其中，分别针对冻毛毛雨（freezing drizzle, FZDZ）、冻雨（freezing rain, FZRA）的SLD不同形成条件拟合的结果进行了对比罗列。可以看出拟合结果与SLD结冰条件滴谱吻合良好。可以看出抽样的粒径分布与附录O曲线吻合较好，有效模拟了“双峰分布”的规律。同时可以看出大水滴整体分布较稀疏，粒径域宽。图2附录O条件下4种SLD环境的粒径分布拟合结果（FZDZ代表冻毛毛雨，FZRA代表冻雨）Fig.2Fitted droplet size distribution of four SLD conditions in Appendix O(“FZDZ” represents thefreezing drizzle, “FZRA” represents the freezing rain)10.12050/are20220308.F2a110.12050/are20220308.F2a2(b)体积频度曲线（蓝色曲线代表附录O给出的水滴粒径分布关系，红色柱状图为本文对曲线进行的粒径离散拟合结果） 10.12050/are20220308.T001表 1粒径分布的拟合参数Table 1Fitting parameters for droplet size distributionsSLD形成条件颗粒群ωMVD/μmn冻毛毛雨小粒径水滴0.82212.6MVD40μm大粒径水滴0.181091.3冻毛毛雨小粒径水滴0.32242.1MVD40μm大粒径水滴0.682181.8冻雨小粒径水滴0.70172.0MVD40μm大粒径水滴0.306832.0冻雨小粒径水滴0.24221.5MVD40μm大粒径水滴0.767922.11.3　水滴飞溅模型本文通过OpenFOAM中内嵌的水滴飞溅模型[28] 来实现大水滴撞击飞溅的模拟。模型参数采用了LEWICE模型参数[29] 。其中，飞溅临界撞击系数、质量损失率和飞溅水滴大小为Kctr=Ky(sinθ)1.25 （5）f=0.7(1-sinθ){1-exp[-0.0092(Kctr- 200)]} （6）dsd0=8.72e- 0.028K(0.05≤dsd0≤1) （7）式中：Ky=A- 38(Oh- 25Wen)45为撞击系数，Oh=μd/ρwσd为Ohnesorge数；A=1.5(LWC⋅ρw- 1)13为水滴入射频率；θ为水滴入射角；Wen=ρwud,n2d/σ为第n个水滴的撞击韦伯数，ρw，ud,n2，d和σ分别为水密度、水滴垂直撞击速度、直径和表面张力。当KctrKy时，收集的水滴质量需要乘以系数f，否则不变。1.4　流动、传热和结冰模型本文利用OpenFOAM平台内嵌水膜模型[27] 实现流动与传热模型，结冰模型为Myers模型[30] 。结冰时的质量和能量守恒方程ρi∂B∂t+ρw∂h∂t=Mcollected (8)ρilf∂B∂t=ki∂T∂z-kw∂θ∂z (9)式中：B和h分别为冰高和水膜厚度；ρi和ρw分别为冰和水的密度；Mcollected为该网格单位时间的收集水量；lf为单位质量水的结冰潜热；ki和kw分别为冰和水的导热系数；T和θ分别为冰层和水膜内部的温度；z为冰生长高度。表面未冻结水膜动量方程为∂ρwhU∂t+∇ρwhUU=- h∇(pimp+pspl)+                                                      Simp+Sspl (10)式中：U为水膜速度矢量；pimp和pspl为水滴撞击/飞溅的压力源项；Simp和Sspl为水滴撞击/飞溅动量源项。在求解水膜方程时，由于水膜内速度、温度等变化均为线性变化，故以一层虚拟网格对其进行计算。空气流场与水膜有共同的边界条件。首先得到空气流场中水的质量、动量和能量于水膜边界处的值，传入水膜得到方程源项，再通过求解水膜方程传出水膜厚度、温度与速度。冰层推移速率为ρilf∂B∂t=kiTf-Tsb+kwQcf+Qd+Qe-Qaf-Qkkw+Qcf+Qd+QeTf-Tah+                (Qout-Qin) (11)式中：Tf和Ts分别为冰的熔点和壁面温度；Qaf和Qcf分别为空气对水层的气动加热和对流换热热流；Qk为水滴动能；Ql为该网格单位的结冰潜热；Qd为水滴显热；Qe为蒸发热流；Qout和Qin为流出/流入水的热量；kw为水的导热系数；Ta为空气温度。将式（11）代回质量守恒式（8）积分可得水膜厚度改变量。当结冰量超过落下水滴和流入水质量，结冰进入霜冰状态，生成的冰质量等于进入网格的水质量。2　结冰计算对比校验2.1　校验算例和条件为校验计算方法的可靠性，本文选择了参考文献[31]和[1]中的冰风洞试验结果作为对比算例，其中包含了小水滴云雾的霜冰和明冰（见表2），以及SLD结冰（见表3）。结冰外形均为弦长0.5334m的NACA0012翼型。由于两篇文献并未给出云雾环境的具体粒径分布情况，故这里使用其给出的MVD作为单一水滴粒径计算。10.12050/are20220308.T002表2小水滴结冰对比算例工况[31]Table 2Icing conditions of cases with small droplets[31]编号结冰时间/s速度/(m/s)迎角/(°)温度/℃含水量/(g/m3)粒径/μm136067.054-19.4120236067.054-4.412010.12050/are20220308.T003表3SLD结冰对比算例工况[1]Table 3Icing conditions of cases with SLDs[1]编号结冰时间/s速度/(m/s)迎角/（°）温度/℃含水量/(g/m3)粒径/μm3804510-19.60.91704336774-19.21.04160本文采用的结冰计算方法包括流场和结冰过程均为三维，但对比的冰风洞试验结果为近似二维流动条件，冰外形结果也为二维曲线。故我们在计算时将对比算例的二维翼型在展向上拉伸一个弦长的长度，用二维工况计算得到2.5维冰形。由于大水滴的空间分布稀疏性，冰形在展向分布上会有一定随机性。为了减少冰形的随机性，在校验时在冰形沿展向于不同截面上截取多个冰外形曲线，将其坐标点进行平均处理后再和试验冰形进行对比。2.2　冰形对比情况2.2.1　小水滴云雾结冰冰形对比情况针对小粒径水滴霜冰和明冰条件的计算与对比可以验证本文算法在基本结冰算法功能和溢流传热方面的计算可靠性，其计算工况见表2，结果如图3和图4所示。计算结果与试验结果吻合良好，表明小水滴条件下的收集、溢流传热等计算可靠。10.12050/are20220308.F003图3算例1的冰形对比情况（霜冰）Fig.3Comparison of Ice shapes for case 1 (rime ice)10.12050/are20220308.F004图4算例2的冰形对比情况（明冰）Fig.4Comparison of Ice shapes for case 2 (glaze ice)2.2.2　SLD结冰冰形对比情况在小水滴结冰计算校验基础上，进一步计算和对比校验SLD结冰冰形可以有效验证本算法关于SLD动力学模型的计算准确性。SLD校验算例选择了参考文献[1] 中条件差异较大的两个冰风洞试验结果，以更好验证算法在不同条件下的计算可靠性，其工况见表3。图5（a）和图6（a）给出了平均冰形与试验冰形的对比结果。算例3的计算结果与试验结果吻合较好，算例4的计算结果中结冰厚度比试验冰形稍大，基本轮廓相似。图5算例3的冰形对比情况Fig.5Comparison of ice shape for case 310.12050/are20220308.F5a110.12050/are20220308.F5a2 图6算例4的冰形对比情况Fig.6Comparison of ice shape for case 410.12050/are20220308.F6a110.12050/are20220308.F6a2 SLD的冰形在展向上有较明显的随机性，并且与MVD成正相关，如图5（b）和图6（b）所示。由于飞机结冰数值计算中水滴收集、传热和结冰等过程均只考虑来流的平均值，理论上不应产生随机偏差。这说明大水滴分布会导致冰形随机性。但现有数值模拟方法还无法准确计算该粗糙特征[33] 。3　SLD冰风洞试验对比分析3.1　试验条件为进一步校验本算法并且分析SLD粒径组分的影响，本节基于航空工业气动院国内首次实现的SLD冰风洞试验进行对比校验和分析。这里选择了其调试的三种粒径分布曲线和6个试验条件进行计算对比。结冰外形为0.6m弦长的NACA0012翼型。如图7所示，三种云雾条件的水滴平均粒径均在20μm左右，均为MVD40μm的冻毛毛雨条件，但有着不同的LWC和大水滴占比。其中，云雾1的分布最接近附录O中冻毛毛雨分布（MVD40μm），LWC则接近云雾3，为0.32g/m3；云雾2的LWC较大（为0.49g/m3），分布曲线形状接近云雾1；云雾3的大水滴占比最大，粒径分布接近于MVD40μm的冻毛毛雨分布，LWC为0.37g/m3。算例5~7重点对比三种云雾分布的结冰形状，以分析SLD条件下LWC和粒径分布对冰形的影响。算例8~10则基于云雾1分析了霜冰到明冰条件的冰形变化。冰形计算对比条件见表4。10.12050/are20220308.F007图7SLD冰风洞试验的云雾粒径分布Fig.7Cloud and fog particle size the diameter distribution in the of SLD ice-cloud in icing wind tunnel experiment10.12050/are20220308.T004表4气动院SLD结冰风洞试验工况Table 4Conditions of SLD icing wind tunnel experimentsby AVIC ARI编号结冰时间/s速度/(m/s)迎角/(°)温度/℃云雾平均粒径/μm5815950-15.01196539950-15.02207710950-15.03238815950-6.61199815950-10.011910815950-18.31193.2　冰形对比情况这里应用第1.2节介绍的粒径分布模拟方法拟合了图7所示的云雾粒径分布函数，并进行冰形计算。由图8可见，6个工况下的计算结果能与试验结果较好符合，只有图8（a）的上冰角与试验相差较多。图8（d）~图8（f）显示的不同温度冰形表明算法可以较好地模拟不同温度的SLD结冰。10.12050/are20220308.F008图8计算与试验的SLD冰形对比情况Fig.8Comparison of the comparison between calculated and experimental SLD icing wind tunnel results对比图3、图4和图8可看到，SLD环境结冰的显著特征是冰形后部存在大量随机突起，而小水滴条件下冰形相对平滑。如图8（a）~图8（c）所示，更大的LWC和大水滴组分都可产生较大的冰角和随机凸起。目前，算法虽然能一定程度上模拟冰形表面随机性，如图6（b）所示，但与实际冰形相比还有很大差距。这是SLD结冰模拟还需要进一步研究的方向。3.3　冰形对SLD粒径组分的敏度讨论由于云雾1最接近14CFR-25部附录O中MVD40μm的冻毛毛雨环境，从图8（d）~图8（f）可以明显看到该环境的结冰特征：更接近小粒径云雾结冰，没有很明显的后部随机凸起。值得注意的是，云雾1和3有着十分接近的LWC和MVD，但两者的冰形差异极大，如图8（a）和图8（c）所示。两种云雾条件最大的差异在于环境中大/小水滴组分比例（约20%），如图7所示。由此可知，SLD云雾中的水滴组分对冰形有显著的影响，这将造成飞机带冰气动性能评估的较大偏差[34] 。而在SLD结冰风洞试验中，水滴粒径分布的精确模拟在国际上也未能全面实现[3] 。SLD粒径分布偏差对冰形的影响规律值得探讨。4　SLD粒径分布影响冰形规律分析4.1　SLD分布环境模拟精度对冰形的影响为量化SLD粒径分布模拟精度对冰形的影响，本节基于14CFR-25部附录O条件对比了不同粒径分布模拟方式得出的冰形差异。对比条件包括：（1）按单一粒径计算（以云雾整体MVD和LWC模拟）；（2）双粒径组合（大、小双峰分别取MVD和LWC后叠加）；（3）按图2的方式拟合粒径分布并抽样。通过对比冰形可评估SLD粒径分布对冰形的影响。本处计算条件以工况一为基础，并为了突出冰形的差异将结冰时间延长到600s，具体见表5。10.12050/are20220308.T005表 5云雾粒径分布影响冰形的计算工况Table 5Conditions of icing evaluation on the effect ofdroplet size distribution编号结冰时间/s速度/(m/s)迎角/(°)温度/℃含水量/(g/m3)1160067.054-101图9显示了不同粒径模拟方式计算得到的冰形。当MVD40μm时，以云雾MVD单一粒径计算的冰形较另外两种模拟方式的小很多。双粒径组合与完整分布模拟的冰形结果有一定差异，但整体形状接近。而当MVD40μm时，单一粒径模拟与另两种模拟方式的冰形差异较小，特别是图9（d）显示的冻雨结冰情况。这说明在MVD较小时，较少的大水滴组分可对冰形造成极大改变。推测在SLD结冰过程中，大、小水滴组分的质量并非线性叠加到冰形上，而是有明显的耦合作用。故精确模拟SLD的粒径分布是准确计算冰形的基本条件。10.12050/are20220308.F009图9三种云雾分布模拟方法的冰形对比Fig.9Ice shape comparison between methods using droplet distribution and single/double MVD(s)不过，从图9（b）和图9（d）中也可看出，在大水滴为主的条件下，水滴粒径分布形式对冰形的影响将减弱。此时降低SLD粒径分布模拟精度不会导致太大的冰形模拟精度损失。4.2　粒径分布组合对收集率的影响为理解水滴粒径分布对冰形的影响规律，本节计算并对比了几种典型粒径及其组合的收集率。为了更直观地观察水滴直径对收集的影响，这里的计算未使用破碎、飞溅等水滴动力学模型。图10中显示了20μm、100μm、200μm的收集率及其组合的情况。其计算条件与表5一致。图10（a）显示20μm和100μm水滴的收集率差异很大，而100μm和200μm水滴收集率十分接近。这是因为超过50μm的大水滴惯性明显，不易在翼型扰流流场中改变轨迹。在此流场条件下，100μm和200μm粒径的水滴应有相近的飞行轨迹。因此，当水滴粒径增加到100μm以上时，粒径的改变对收集率贡献很小。10.12050/are20220308.F010图10不同单一粒径水滴与组合后收集率Fig.10Collect efficiencies of droplets with three diameters and their composed clouds从图10（b）中可见，当20μm粒径的云雾中加入20%的100μm水滴后收集率曲线发生很大改变，包括驻点处收集率峰值和两侧收集范围。收集率的变化幅度取决于大水滴组分的比例。而在100μm粒径的云雾中增加200μm粒径水滴组分时，收集率结果并未有明显改变，如图10（c）所示。这一结果与图10（a）显示的结果相符。水滴收集是冰形的直接质量来源。特别在明冰条件下，驻点附近收集率的增加相当于水含量增加，导致溢流增强而进一步改变冰角形状。因此，当MVD较小时大水滴组分的比例对冰形有极大影响。而当MVD大时，大水滴组分的影响将不明显。因此在SLD环境的数值和试验模拟中，当MVD40μm时云雾粒径分布的精确度需要重点保证；而在MVD40μm时粒径分布模拟精度可适当放宽。4.3　SLD分布环境结冰的计算效率讨论目前的主流结冰软件（如FENSAP-ICE）只能计算单一粒径水滴条件结冰，在面对SLD环境的复杂粒径分布时只能将其人为划分为多个离散组分，分别计算后再将水滴收集率叠加。考虑到SLD中50μm以上大水滴粒径的宽广分布，必须要划分5个以上组分才可满足模拟精度要求。考虑到需要人工介入计算过程，所需计算时间将远大于5倍单粒径条件计算时间。本文方法可以在同样收集率计算精度下一次完成任意粒径分布的计算。相比传统方法的计算效率提升不言而喻。5　结论本文针对SLD双峰分布环境提出了水滴粒径分布的抽样混合结冰计算方法，同国内首次实现的MVD40μm冻毛毛雨条件结冰风洞试验结果进行了对比，并探讨了SLD粒径分布对冰形的影响规律。主要结论如下：（1）水滴粒径分布抽样混合算法可高效准确地模拟SLD双峰分布云雾结冰。冰形计算结果在多种条件下均与冰风洞试验结果符合较好，表明本算法可有效预测MVD40μm的冻毛毛雨结冰冰形。（2）在冻毛毛雨条件下，水滴粒径分布对冰形有显著影响。在MVD不变的情况下，20%的大/小水滴组分差异可导致冰形张角与高度50%以上改变。数值计算显示，在MVD较小时，SLD云雾结冰冰形对水滴组分的变化十分敏感。但随着云雾的MVD增加，粒径组分对冰形的影响逐渐减弱。（3）导致冰形对SLD组分敏感性的原因是，小水滴为主的条件下大水滴组分的加入极大地改变了水滴收集范围和质量。而在水滴总体粒径大时分布的改变对水滴收集影响较小。（4）在冰风洞试验和计算中，模拟MVD40μm的SLD环境需要精确模拟云雾粒径分布。而在模拟MVD40μm的SLD环境时，可适当放宽分布的模拟精度。