李伟斌，魏东，杜雁霞, *，易贤，杨肖峰

1.中国空气动力研究与发展中心 空气动力学国家重点实验室, 绵阳 621000 2.中国空气动力研究与发展中心 计算空气动力研究所, 绵阳 621000

过冷水滴撞击低温基底在满足一定条件后会发生冻结，随着夹杂过冷水气流的不断撞击，基底表面形成越来越厚的结冰。与传统结冰较为不同的是，这种结冰具有典型的动态过程，在微观上，表现为水滴的不断冻结累积，并且相互间形成孔隙。由于结冰条件的不同，结冰微观结构中的孔隙特性会存在不同程度差异，直接影响结冰的热/力学特性。

在结冰冰形计算[1-2]中，结冰的密度是关乎冰形结果的重要参数；在结冰探测[3-4]中，波速与结冰的微观结构密切相关；在防除冰模拟与试验[5-6]中，黏附力、导热系数等受结冰的致密程度影响。然而这些本质由结冰微观孔隙结构决定的参数大多采用经验设置方式给定，缺少定量化的手段，这极大限制了冰形精确预测、结冰高精度探测、防除冰系统的高效设计等。

目前，结冰微观的定量研究多集中于晶核的形成、晶体生长、表面接触形核等数值模拟[7-12]，虽然有利用孔隙率和分形维数等对结冰微观结构进行定量刻画[13-14]，但是缺少类似于多孔材料模型[15-17]的系统研究，无法实现对不同结冰条件下动态结冰孔隙的数学化表达。

本文针对动态结冰微观孔隙结构定量研究的不足，以不同结冰条件下的结冰显微图像为对象，采用前期所提的变分分割模型[18-19]，对图像中的孔隙进行分割提取，并分析了所使用分割方法的适用性；以孔隙提取结果为基础，系统分析并提出了孔隙形态、孔径分布等相关结构特征的数学表述方式。通过研究，以期为动态结冰微观结构的数学模型建立提供理论依据和数据支撑。

1 结冰试验手段及方法

1.1 结冰试验

本文动态结冰试验在中国空气动力研究与发展中心(CARDC)小型结冰风洞中进行，其试验段尺寸为0.3 m×0.2 m×0.65 m(宽×高×长)，如图1所示。所使用的喷雾用水经过有效过滤，含杂质较少。试验所用结冰模型为超临界翼型剖面的缩比模型，弦长100 mm，展长200 mm，材料为7075铝合金，导热系数为130 W/(m·K)，如图2所示。

图1 0.3 m×0.2 m×0.65 m结冰风洞Fig.1 0.3 m×0.2 m×0.65 m icing wind tunnel

图2 试验模型Fig.2 Test model

1.2 结冰微观图像采集

结冰显微图像由电子显微镜连接电脑采集得到，显微镜型号为Olympus CX31，如图3所示。由于基底表面和结冰表面的不同，动态结冰基底部分和其余部分具有明显不同，本文主要讨论动态结冰过程的微观孔隙结构，为排除基底对结冰的影响，试验结冰选取远离基底的部分。结冰切片是低温环境下从翼型结冰中切割而成，且表面处理光滑平整。文中试验所使用的结冰显微图像均为40倍放大率下采集得到。

图3 结冰显微图像采集系统Fig.3 Microscopy image acquisition system of icing

2 结冰微观孔隙结构提取方法

结冰微观结构的不同会直接影响结冰的热/力学特性，开展结冰微观孔隙结构提取是定量认识其变化规律的前提。

首先，采用前期所提变分分割方法[18]对结冰显微图像进行图像处理，获取孔隙的边缘曲线。变分分割方法的做法是：① 提出分割思想，并将其数学化，最终得到相应的最小化问题；② 结合变分方法和水平集方法，应用优化方法对所提最小化问题进行求解。变分分割方法的求解初值为一条曲线，求解过程的中间量对应着演化曲线，所得到的解与最终的分割曲线对应。

定义u0:Ω→R为灰度图像，Ω为矩形区域，表示图像定义域，R为实数域。

文献[18]的出发点是分割出图像的背景部分，便可以得到需要的目标区域，即通过演化迭代得到一个足够大的外部区域，并且该区域内的灰度值是相同的，它可以数学化表示为带约束的优化问题,即

min -S(Ω1)

(1)

式(1)限制条件的目的是保证曲线外部灰度的同一性，但是由于自然图像中背景部分的灰度值很难达到相同，因此采用Lagrange乘子法[20]对限制条件进行弱化，即

(2)

式中：μ>0为Lagrange乘子；l为区域Ω1的边界长度，目的是减小图像噪声带来的影响，即抑制噪声点处产生小曲线，υ为其权重。

式(2)不能直接求解，需要将其数学化表示，应用水平集方法[21]，可以将演化曲线视为水平集函数φ的0-水平集，即集合{x|φ(x)=0}，将外部区域表达为：Ω1={x|φ(x)<0}，同时应用如下一维Heaviside函数，即

(3)

式中：h为自变量。便可最终得到与式(2)对应的最小化问题，即

(4)

关于φ对式(4)进行求导，可以得到其对应的Euler-Lagrange方程，再应用最速下降法，便可得到对应的水平集演化方程为

(5)

式中：δ(φ)为Dirac函数，仅在φ=0处有非零取值。因此在迭代求解式(5)时，两个相邻迭代步的值仅在φ=0处有较小变化，而φ=0对应着曲线的位置，也就是说两个迭代步中的曲线位置变化较小，曲线演化速度较慢，这不利于曲线收敛至目标区域的边界。为了克服δ(φ)对迭代过程带来的影响，将其做近似处理，令δ(φ)≈1，最终得到求解的水平集演化方程为

(6)

求解式(6)，便可得到上述变分分割模型的收敛解φ*。基于此，可以得到孔隙的边缘曲线C*={x|φ*(x)=0}及孔隙区域Ω*={x|φ*(x)>0}。进而结合MATLAB的bwlabel()和regionprops()等命令，对结冰中所有的孔隙进行编号，得到各自对应的面积和周长等相关定量信息。这样就实现了结冰孔隙提取的目的。

3 动态结冰孔隙结构定量结果及分析

在结冰显微图像采集过程中，可见光从下至上，穿透结冰，并进入显微镜镜头部分，从而使成像系统清晰捕捉结冰的放大图像。由于这一成像特点，使得光线经过结冰的气泡孔隙时，发生折射现象，从而在所成图像中呈现出黑色散斑。在结冰微观研究中，孔隙分布的相关特性是分析工作的基础。本文开展了不同温度下的结冰微观孔隙结构定量分析研究工作，其中水滴平均粒径MVD=38 μm，液态水含量LWC=0.75 g/m3，速度v=25 m/s，温度T在试验分析中给出。

3.1 动态结冰显微图像分割结果

图4 阈值法和变分分割方法分割结果对比Fig.4 Comparisons of segmentation results by thresholding and variational methods

为了验证本文分割方法的优越性，将其与传统的阈值法进行了分割效果对比，结果展示于图4中。阈值法是以某一给定像素值为阈值，将待分割区域二值化，从而得到分割结果。本试验中，选择的是彩色图像的第3通道作为分割对象，为了对比的统一性和展示效果，将阈值法的二值化结果转换为类似于变分分割方法的分割曲线形式。阈值法的关键在于阈值的选择，图4(a)～图4(c)给出了3个阈值的分割结果。由于光线传播的特点，某些孔隙中部会显示出亮点，无论阈值如何选择，都没能正确分割孔隙，即将孔隙内部从孔隙中分离出，没有得到完整孔隙(见图4中白色框对比结果)，并且阈值选择不恰当会将焦平面外的孔隙分割出，得到不正确的分割区域，进而直接影响孔隙分布规律的分析结果。而本文分割方法得到了较为准确的分割结果，如图4(d)所示。

结冰条件的不同会直接导致结冰显微图像内孔隙大小和数量的不同，为了研究它们的变化规律，开展了不同温度条件下的结冰试验，并应用变分图像分割方法，对它们的显微图像进行了分割处理，结果如图5所示。从结果可以看出，每个状态下的孔隙均被有效分割出，并且孔隙形状与圆形较为接近。同时，从结果可以发现孔隙数量、大小和分布确实存在着明显不同。

图5 不同温度的显微图像分割结果Fig.5 Segmentation results of microscopic images with different temperatures

3.2 动态结冰孔隙形态

动态结冰内部孔隙的形态是对微观结构进行合理数学建模并分析其结构特征的基础，应用第2节的变分分割方法对动态结冰显微图像进行分割处理，讨论孔隙形态。采用式(7)的圆形度表征量[22]对结冰孔隙的截面形状进行分析，即

C=4πA/L2

(7)

式中：A和L分别为孔隙截面的面积和周长；C值表征的是区域近似于圆形的程度。显然，当区域为圆时，式(7)具有最大值1。当区域越近似于圆，值越趋近于1；当区域越复杂，C值越小。图6中给出了同一放大倍数下，不同温度、不同结冰区域微观图像中的孔隙截面对应的C值，其中横坐标n为孔隙序号。同时表1中给出了C值的相关统计信息。从表中数据可以看出，结冰中C>0.785 4(正方形)的区域占比较大，且所有C值的均值和中值均接近于1，即大多数孔隙截面近似于圆形，这点也可以从图5的分割结果和动态结冰过程孔隙受力情况中得到验证。基于此，结合显微图像采集的广泛性及随机性，可以推广得到孔隙任意截面均为圆形的结论。因此，在对结冰微观结构进行量化描述过程中，将结冰孔隙视为球形是可行的。

图6 结冰孔隙截面圆形度Fig.6 Circularity of pores cross-section of icing

表1 图6数据的统计信息Tabel 1 Statistical information of data in Fig.6

总个数C>0.7854的个数/占比C均值C中值321240/0.74770.83240.8815

3.3 动态结冰孔隙孔径分布

从图5可以看出不同温度下的孔隙孔径的分布特性存在着明显不同，为了更清楚认识这种不同，对其进行定量分析。在孔隙呈球形的结论之上，考虑其区域的复杂性，基于面积和周长数据，通过式(8)计算动态结冰中不同孔隙的直径:

(8)

对大量孔隙直径进行统计，图7给出了不同区间无量纲孔径频率f的直方图。从图中可以看出，在一定区间范围内，孔隙直径均完全覆盖，即其可以取一定范围内的任意值。理论上，在结冰过程中，主要受结冰热/力学的影响，孔隙的形成具有随机性，并且其大小可以为任意实数。因此，孔隙孔径的取值是连续的。

从图7中可以看出不同孔径的分布密度是不同的，为了更清晰、量化认识这种不同，基于结冰显微图像的分割结果，计算直径的分布函数F(d)，并应用式(9)对其进行曲线拟合，即

(9)

式中：k1、k2均为正实数；x∈[0,+∞)。图8给出了不同温度条件下的直径分布函数和拟合函数对应的曲线图(参数k1、k2的取值见表2)，从图中可以看出所有拟合曲线与孔径分布函数吻合较好。因此，孔隙的孔径分布函数可以用式(9)的形式表达。进而，结合孔隙孔径取值是连续的结论，可以推导其概率密度函数(Probability Density Function, PDF)为

图7 结冰孔隙直径的频数直方图Fig.7 Frequency histogram of pore diameters of icing

图8 结冰微观孔径分布曲线及其拟合曲线Fig.8 Distribution curves of micro pore diameter of icing and the fitting curves

(10)

图9中给出了图8对应状态的概率密度函数拟合曲线图，可以看出随着温度的降低，孔径最大值增大，概率密度峰值减小。这是因为在保持其他结冰条件不变的情况下，温度越低，过冷水滴撞击基底表面结冰速率更快，溢流效应减弱，水滴之间以及水滴与结冰之间夹杂的气泡更不易破碎或者逃逸。

式(9)是所提出的关于孔径的分布函数，其在不同结冰条件下的参数取值(k1,k2)不同。为了验证该形式在同一结冰状态下的适用性，分别选取温度为-4.8 ℃和-11 ℃动态结冰不同位置的结冰，进行微观孔径分布分析，结果展示于图10中。从图中可以看出，同一结冰条件不同结冰位置的孔径分布与所提出的分布函数具有较好的吻合度，说明以式(9)和式(10)定量表征结冰孔径分布规律是可行的。

表2 图8中拟合函数参数取值Tabel 2 Parameter values of fitted functions in Fig.8

图9 结冰微观孔径概率密度函数拟合曲线Fig.9 Fitted PDF curves of micro pore diameter of icing

图10 多位置结冰微观孔径分布曲线及 其拟合曲线Fig.10 Distribution curves of micro pore diameter of icing on different locations and the fitting curves

4 结 论

1) 在结冰微观孔隙提取方面,相比于传统方法，文中使用的变分分割方法能够得到更准确的分割结果，能够保持孔隙的完整性，对结冰微观结构的定量分析提供了较好的数据输入。

2) 在结冰微观孔隙结构定量研究方面。通过对结冰内部大量孔隙形态、直径等相关信息进行统计分析，可以得到微观孔隙近似呈球形、孔径分布连续等结论。并在此基础上提出孔径分布函数表达式，试验结果表明所提表达式与结冰内部孔隙孔径分布较为吻合。

本文对动态结冰孔隙形态和孔径的分布进行了初步的定量分析讨论，下一步将基于更完整的结冰风洞试验数据，定量确定参数k1和k2与结冰条件的关系，进而建立动态结冰三维结构数学模型，提取用于表征其微观特征的定量指标，最终建立结冰微观与宏观间的联系。

[1] MESSINGER B L. Equilibrium temperature of an unheated icing surface as a function of airspeed[J]. Journal of the Aeronautical Sciences, 1953, 20(1): 29-42.

[2] MYERS T G, CHARPIN J P F. A mathematical model for atmospheric ice accretion and water flow on a cold surface[J]. International Journal of Heat and Mass Transfer, 2004, 47 (25): 5483-5500.

[3] IKIADES A A, KONSTANTAKI M， CROSSLEY S D. Fiber optic sensors technology for air conformal ice detection[C]∥Optical Technologies for Industrial, Environmental, and Biological Sensing. International Society for Optics and Photonics, 2004: 357-368.

[4] 尹胜生, 叶林, 陈斌, 等. 可识别冰型的光纤结冰传感器[J]. 仪表技术与传感器, 2012(5): 9-12.

YIN S S, YE L, CHEN B, et al. Fiber-optical icing sensor for detecting the icing type[J]. Instrument Technique and Sensor, 2012(5): 9-12 (in Chinese).

[5] IVALL J, RENAULT-CRISPO J S, COULOMBE S, et al. Ice-dependent liquid-phase convective cells during the melting of frozen sessile droplets containing water and multiwall carbon nanotubes[J]. International Journal of Heat and Mass Transfer, 2016, 101: 27-37.

[6] KEITZL T, MELLADO J P, NOTZ D.Impact of thermally driven turbulence on the bottom melting of ice[J]. Journal of Physical Oceanography, 2016, 46(4): 1171-1187.

[7] YE Y, NING N, TIAN M, et al. Nucleation and growth of hexagonal ice by dynamical density functional theory[J]. Crystal Growth & Design, 2017, 17(1): 100-105.

[8] BLAKE J, THOMPSON D, RAPS D, et al. Simulating the freezing of supercooled water droplets impacting a cooled substrate[J]. AIAA Journal, 2015, 53(7): 1725-1739.

[9] LU Y, ZHANG X, CHEN M. Size effect on nucleation rate for homogeneous crystallization of nanoscale water[J]. The Journal of Physical Chemistry B, 2013, 117(35): 10241.

[10] COX S J, RAZA Z, KATHMANN S M, et al. The microscopic features of heterogeneous ice nucleation may affect the macroscopic morphology of atmospheric ice crystals[J]. Faraday Discussions, 2013, 167(1): 389-403.

[11] ZHANG X X, LU Y J, CHEN M. Crystallisation of ice in charged Pt nanochannel[J]. Molecular Physics, 2013, 111(24): 3808-3814.

[12] LUPI L, HUDAIT A, MOLINERO V. Heterogeneous nucleation of ice on carbon surfaces[J]. Journal of American Chemical Society, 2014, 136(8)：3156-3164.

[13] LEI G L, DONG W, ZHENG M, et al. Numerical investigation on heat transfer and melting process of ice with different porosities[J]. International Journal of Heat and Mass Transfer, 2017, 107: 934-944.

[14] 杜雁霞, 桂业伟, 柯鹏, 等. 飞机结冰冰型微结构特征的分形研究[J]. 航空动力学报, 2011, 26(5): 997-1002.

DU Y X, GUI Y W, KE P, et al. Investigation on the ice-type microstructure characteristics of aircraft icing based on the fractal theories[J]. Journal of Aerospace Power, 2011, 26(5): 997-1002 (in Chinese).

[15] FRANK C. Reservoir formation damage: Fundamentals, modeling, assessment and mitigation[M]. 2nd ed. Burlington: Gulf Professional Publishing, 2007.

[16] LEHMANN P, STAHLI M, PAPRITZ A. A fractal approach to model soil structure and to calculate thermal conductivity of soils[J]. Transport in Porous Media, 2003, 52(3): 313-332.

[17] OKABE H, BLUNT M J. Pore space reconstruction of vuggy carbonates using microtomography and multiple-point statistics[J]. Water Resources Research, 2007, 43(12): W12S02.

[18] LI W B, YI X, SONG S H. Convex background removed model for image segmentation using the split Bregman method[J]. Journal of Information and Computational Science, 2015, 12(17): 6641-6652.

[19] 李伟斌, 易贤, 杜雁霞, 等. 基于变分分割模型的结冰形测量方法[J]. 航空学报, 2017, 38(1): 120167.

LI W B, YI X, DU Y X, et al. A measurement approach for ice shape based on variational segmentation model[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(1): 120167 (in Chinese).

[20] BERTSEKAS D P. Constrained optimization and Lagrange multiplier methods[M]∥Computer Science and Applied Mathematics. Boston: Academic Press, 1982: 1.

[21] OSHER S, FEDKIW R. Level set methods and dynamic implicit surfaces[M]. New York: Springer-Verlag, 2002: 51-72.

[22] 王东霞, 宋爱国. 基于三坐标测量机的圆度误差不确定度评估[J]. 东南大学学报(自然科学版), 2014, 44(5): 952-956.

WANG D X, SONG A G. Uncertainty assessment of circularity error based on coordinate measuring machine[J]. Journal of Southeast University (Natural Science Edition), 2014, 44(5): 952-956 (in Chinese).