動態結冰孔隙結構三維建模方法

李偉斌，宋超，易賢，馬洪林，杜雁霞

（1 中國空氣動力研究與發展中心計算空氣動力研究所，四川綿陽621000； 2 中國空氣動力研究與發展中心飛行器結冰與防除冰重點實驗室，四川綿陽621000）

引 言

動態結冰是指結冰所需的液態水由氣流輸運而來，并在液態水運動撞擊基底的過程中不斷凍結而動態形成的。相比于靜態結冰，正是由于這一“動態”特點，使得在凍結過程中，水膜與冰層間的空氣以及部分水滴中的空氣來不及逃逸，被“囚困”于結冰中，形成內部孔隙結構[1]。這種孔隙結構特征決定了結冰的許多重要物理性質，比如密度、熱傳導能力、強度等[2-4]。動態結冰的這一孔隙結構因結冰條件的不同而具有可見的差異，快速獲取動態結冰孔隙結構定量信息是根本上認識其動態形成過程的內在機制，也是結冰熱/力學相關方向研究的基礎[5]。

目前關于結冰微觀特性的研究多集中于數值模擬。在結冰生長微觀研究方面，動態結冰被認為主要包括晶核的形成和晶體生長兩個階段[6]，其中第一階段是過冷水滴熱力學非平衡態到平衡態轉變的過程[7]；第二階段是液/固界面推進直至完全凝固的過程[8-9]。在飛機結冰中，第一階段明顯快于第二階段[10]，這也導致了早期飛機結冰研究忽略形核這一重要現象。同時，應用分子動力模擬結冰的全過程也迅速發展[11-12]，表明形核是足夠數量的氫鍵結合從而形成第一階段的初始核，進而受熱力學行為影響，初始核改變形狀和大小，快速生長直至最終的全結晶。在基底表面接觸結冰方面，針對不同表面、不同粗糙度的接觸表面，現有研究從分子和納米尺度，對形核開展了數值模擬，結果表明相對于光滑表面，粗糙表面不易形核[13]；同樣地，曲率大的表面也會減少形核速率[14-16]。這一點在結合電子顯微成像系統開展的超疏水材料涂層實驗中得到了驗證[17]，因為這種材料被認為是在原有基底表面形成了一層微納米粗糙的特征，從而延緩了形核凍結的時間。正因為如此，超疏水材料對于減小結冰表面的黏附作用具有重大貢獻[18-19]。另外，大水滴靜止結冰[20]和碰撞結冰[21-22]的熱力學和動力學行為也得到了研究，表明過冷態溫度越大、水滴尺度越小或接觸面溫度越低，水滴的相變速率越高；提高表面疏水性，不僅能夠減小水滴碰撞后的鋪展濕潤面積，而且能夠縮短水滴同壁面的接觸時間，從而降低水滴內部的相變速率，這點也在大水滴結冰實驗中得到了較好的驗證[23]。

現有針對動態結冰微觀特性的實驗研究和模擬研究已經較多，研究表明影響動態結冰微觀結構的結冰條件主要有速度、溫度、粒徑、液態水含量等，其中溫度越低，孔徑最大值越大，孔隙率越大[3]；速度越高，來流對氣泡的粉碎作用越強，因而氣泡體積也相對越小，孔隙率越小[24]。這些研究對結冰微觀特性有了一定的定性理解，然而缺少類似多孔材料的細致定量研究。在二維孔隙率分析方面，主要做法是低溫環境下獲取結冰的截面顯微圖像，進而進行相關分析討論[3,24]。在多孔材料孔隙結構定量信息提取研究領域，目前較為有效的途徑是通過建立三維模型進行定量分析，具體做法是利用CT 或者Micro-CT 獲取對象的二維序列圖像，并對其進行圖像處理，再進行三維重建，結果可視化軟件包括OpenGL、Matlab、Tecplot 等，相關研究已經遍歷土壤學[25-26]、材料學[27-28]、醫學[29-31]等。然而由于動態結冰的材料特點，并不能進行相同的掃描建模，因此，其定量信息的提取是動態結冰的難點問題之一[32]。

動態結冰孔隙結構相關研究已經取得了一定的定量信息，然而對于結冰的三維模型建立以及三維信息定量提取方面，缺少必要的研究手段及方法。本文針對動態結冰孔隙結構三維建模問題，基于結冰孔隙形狀、孔徑、分布等二維圖像定量信息，在合理假設的基礎上，研究孔隙的定位及定量方法，探索孔隙模型的量化表征及顯示方法，并驗證方法的可靠性。以期為動態結冰三維刻畫以及特征參數的表征提供新途徑。

1 動態結冰孔隙分布的基本假設

動態結冰具有孔隙結構，宏觀上表現為圖1 的不透明狀，微觀上含有大量氣泡孔隙。李偉斌等[3]在中國空氣動力研究與發展中心0.3 m×0.2 m 的結冰風洞中開展了結冰實驗，并使用Olympus CX31型顯微鏡獲取結冰微觀圖像，采用圖像分割的方法提取結冰微觀圖像中的二維定量信息。基于此，文章從微觀角度入手，通過對不同結冰條件下的結冰、同一實驗狀態不同位置結冰的定量分析，最終得到了動態結冰如下的定量分析結果。

圖1 具有氣泡的動態結冰Fig.1 Dynamic icing with pores

(1) 孔隙呈球形。利用圓形度表征量[33]對孔隙截面進行統計分析，該值大于0.7854(正方形)的區域占比較大，結合顯微圖像采集的廣泛性及隨機性，可以推廣得到孔隙任意截面均為圓形、結冰孔隙為球形的結論。

(2)孔隙直徑取值是連續的。統計大量孔徑數據，得出其可以取一定范圍內的任意值。并且理論上，在結冰過程中，主要受結冰熱力學的影響，孔隙的形成具有一定隨機性，其半徑可以為任意實數。因此，孔隙孔徑的取值是連續的。

(3)孔隙直徑服從特定分布。文獻[3]中所使用的顯微鏡在40倍放大率，對不同結冰狀態下的結冰直徑分布進行曲線擬合，其服從式(1)的特定分布。

其中，x ∈[0,+ ∞)，k1、k2均為正實數，不同來流條件下結冰微觀孔隙的直徑分布對應不同的k1和k2。

(4)孔隙分布位置的隨機性?？紫兜男纬蛇^程受結冰條件和熱力學的影響，雖然呈現一定規律，但具有較強的隨機性。

本文以此作為結冰三維孔隙結構建模的假設，該假設是本文建模方法的重要基礎。

2 動態結冰微觀結構的三維建模方法

為更直觀重建結冰內部孔隙結構，將結冰微觀結構抽象化三維矩陣的形式，基于隨機的形式和孔徑分布函數式(1)確定矩陣各坐標點處的取值，其中以0表示結冰，1表示結冰內部的孔隙。三維建模是問題數學化過程，表示為求刻畫結冰微觀結構的矩陣T:Ω →{0,1}，其中Ω ?Z3，其求解過程如圖2 所示。建模過程主要分為五個步驟：①確定需要生成的結冰的大小(尺寸)；②根據實驗采集到的結冰二維微觀圖像，測量其孔隙率，進而確定三維建模時孔隙的數量N；③根據第1 節中結論(4)隨機生成N個孔隙的位置信息；④根據第1節中結論(3)，生成服從分布式(1)的孔隙直徑信息；⑤最終基于以上步驟生成結冰微觀結構，并予以顯示。這五步對應的詳細過程見2.1～2.5節。

圖2 三維建模方法流程Fig.2 Flow diagram of 3-D modeling method

2.1 三維微觀結構的區域

給定區域Ω =[1,I]×[1,J]×[1,K],Ω ?R3，并令T(Ω) = 0，即以大小為I × J × K 的三維0-矩陣(矩陣元素值均為0)T表示待確定的結冰三維微觀結構區域，該矩陣也可視為是三維圖像。

在微觀結構確定的基礎上，必須給定矩陣大小與實際結冰尺寸的變換關系。將矩陣元素的坐標與長度間建立關系，兩個元素間距離為1，記為1 像素，建模中以比例r(像素/毫米)表示，即r個矩陣元素長度對應真實結冰的1 mm。比例尺r可以由結冰二維顯微圖像中的標定結果確定。如果需要獲取更精細化的三維微觀結構，可以增大r。

2.2 區域內孔隙數量

動態結冰內部由一定數量的孔隙組成，記其數量為N?？紫斗植季哂须S機性，又具有一定的規律性，基于此文中假設三維每個界面的孔隙分布是相同的。那么，三維的孔隙體積占比就等于二維的孔隙面積占比λ2D，稱它們為孔隙率[34]，并以此作為三維建模的數量特征，數量N的確定方式基于式(2)

式中，F-1(1)為式(1)分布函數逆函數在1處的取值；F'(x)為分布函數的導數，即概率密度函數。

由于在計算二維孔隙率時，并未考慮孔隙交融的情況，而在三維孔隙結構建模時會存在這種現象。因此，在建模時需要考慮這點帶來的三維體積損失和二維面積損失，經過多次建模實驗，文中經驗性地選取3%進行修正補償，即

式中，int()表示求整運算。

2.3 孔隙位置

基于孔隙為球形的假設，孔隙的位置可以由球心位置確定。本文球心坐標以隨機的形式確定，記其為Oi(xi,yi,zi),i = 1,…,N。

在區間[1，I×J×K]上以均勻分布方式生成N 個隨機實數{Xi|i = 1,2,…,N}，以式(4)生成孔隙的球心坐標。

2.4 孔徑分布

由式(1)孔隙直徑分布可以直接得到孔隙率，而孔隙率決定了結冰的密度、致密程度及其熱力學特性等，是決定結冰宏觀特性的主要和關鍵因素，是關乎防除冰策略優化的重要因子。因此，獲取準確的孔徑分布(孔隙率)是建立結冰微觀結構三維模型的關鍵步驟，是開展后續研究的前提。

建模中采用的孔徑分布由式(1)得到，而其中的參數k1和k2由二維圖像中得到，具體做法是在低溫環境下制作結冰的切片，在顯微鏡下獲取結冰微觀圖像，進而采用圖像分割的方法獲取孔隙部分，采用Matlab 自帶命令得到各孔隙的周長和面積，進而以孔隙呈球形的假設得到其直徑，最后根據多個孔隙直徑信息確定式(1)中的兩個參數k1和k2。

式(1)是特定實驗條件下提取的分布公式，需在其基礎上轉化為本文的孔徑d，轉換為式(5)

式中，η(像素/毫米)是結冰顯微圖像中像素與真實尺寸之間的比例尺。

具體地，按照均勻分布規則，在區間(0,1)上生成N 個隨機數{Yi|i = 1,2,…,N}，并求出分布函數式(1)所對應的逆{F-1(Yi)|i = 1,2,…,N}，則本文的孔隙直徑為

2.5 三維微觀孔隙結構

更新零矩陣T，生成動態結冰三維微觀孔隙結構，其對應值為1。

以{(Oi,di)}(i = 1,…,N)為對，依次對矩陣中的元素取值進行更新。令

3 建模結果及可靠性分析

以實驗狀態粒徑MVD=38 μm，液態水含量LWC=0.75 g/m3，v=25 m/s，T=-6℃下的結冰二維顯微圖像為例，構造其對應的三維微觀結構。程序在Matlab 7.6中編寫，PC配置為8核CPU 3.4 GHz，內存為16 GB。

3.1 二維定量信息

圖3(a)是該狀態下的顯微圖像，大小為382×288，其與真實尺寸之間的比例尺η ≈27(像素/毫米)。圖3(b)是變分分割方法的分割結果。圖4 是該結冰狀態下直徑分布函數和擬合函數式(1)對應的曲線圖，其中k1= 0.0075，k2= 0.46，且求得孔隙率λ2D= 0.0512。

3.2 三維微觀結構生成

指定結冰矩陣大小I × J × K 為270×270×270，建模中矩陣元素長度與真實結冰尺度之間的比例r = 27(像素/毫米)，即給定大小為270×270×270 的矩陣T = 0，其代表1 cm×1 cm×1 cm。

在結冰二維定量信息和擬生成結冰尺度基礎之上，結合式(2)和式(3)可以確定孔隙的數量為N ≈1663。

在區間[0,2703]上按均勻方式生成N(=1663)個隨機數Xi，i = 1,2,…,N，按照式(4)可以得到孔隙的位置信息。

圖3 結冰顯微圖像及其分割結果Fig.3 Microscopic image of ice and its segmentation results

圖4 結冰微觀孔徑分布曲線及其擬合曲線Fig.4 Distribution curves of micro pore diameter of ice and corresponding curve-fits

在區間[0,1]上按均勻方式生成N 個隨機數Yi，i = 1,2,…,N，按 照 式(1)數 值 求 解 隨 機 數Yi的 逆{F-1(Yi)|i = 1,2,…,N}，以此作為孔隙直徑。

在得到孔隙球心位置及其直徑的基礎上，按照2.5 節掃描的方式對270×270×270 個矩陣元素進行賦值，得到最終的動態結冰三維微觀結構模型，共耗時8.8 s。

圖5 是在Matlab 中對該模型進行圖像顯示，黑色部分表示的是氣泡孔隙；在Tecplot 中對其進行顯示，并按照直徑大小著色，結果展示于圖6(a)，對模型0.25 等值面進行顯示，結果展示于圖6(b)。從圖中可以看出，孔隙分布較隨機，也較均勻，說明本文建模方法具有可行性。

3.3 方法的驗證及可靠性分析

3.3.1 隨機性驗證 圖7是確定孔隙位置和孔徑分別所需的N(=1663)個隨機數Xi和Yi,i = 1,2,…,N 隨序號變化的示意圖，其中橫坐標是隨機數序號i，縱坐標分別是隨機數Xi和Yi的值。可以看出其滿足第1 節中假設(3)和(4)的建模要求，并且具有較好的均勻性和隨機性，這說明以本文方法模擬結冰內部孔隙的隨機性是可行的。

圖5 所生成1 cm×1 cm×1 cm結冰的圖像顯示Fig.5 Visional display of modeling ice with size of 1 cm×1 cm×1 cm

3.3.2 孔徑分布吻合度 圖8是生成的孔徑分布與所服從的分布(1)之間的對比，它們之間的最大絕對誤差為0.02，這說明吻合度較好，進而說明本文孔徑生成方法是切實可行的。

3.3.3 二維信息定量對比 為了直觀顯示所生成的三維結構內部形態，圖9 對比展示了結冰孔隙實驗測量結果及所生成結冰二維截面。圖中可以看出生成結冰的內部孔隙分布與實驗較為相似，表現在孔隙大小、分布、隨機性等方面。

實驗圖片的分割結果中，二維氣泡的孔隙率約為0.0512，為了對比驗證本文建模方法的準確度，對所生成結冰三維微觀結構的所有x,y,z 方向的截面進行孔隙率計算，得到的結果以柱狀圖展示于圖10中。可以看出它們的值主要集中在0.05 左右，變化范圍有限，并且它們的平均值為0.0511，與二維實驗值相近，這說明三維建模得到的孔隙率結果與二維定量結果有較高的一致性。

圖6 所生成1 cm×1 cm×1 cm結冰的孔隙分布Fig.6 Pores distribution of modeling ice with size of 1 cm×1 cm×1 cm

圖7 孔隙位置與直徑對應的隨機數Fig.7 Random numbers for locations and diameters of pores

圖8 生成孔徑分布及其所服從的分布函數Fig.8 Distribution curves of modeling diameter and obeyed function

3.3.4 重復性實驗 為了驗證方法的可重復性，在保持所有輸入條件相同的情況下，開展了100 次建模實驗，并計算了它們的三維孔隙率，結果如圖11所示。可以看出計算結果幾乎全部在區間[0.047,0.054]上，并且它們的平均值為0.0506，與二維孔隙率0.0512 接近，這說明本文方法重復性較好，可靠性較高。

3.4 不同分辨率下的建模結果

考慮實際建模情況，需要生成不同分辨率下的結冰微觀結構模型，即需要對比例因子r 進行適度調整。為了直觀對比分辨率引起的不同，在本文比例尺r=27 的基礎上[圖12(b)]，分別增加了0.5 倍和2倍的縮放模型，結果分別展示于圖12(a)、(c)中，對應真實結冰大小約為1.85 mm×1.85 mm×1.85 mm。從圖中可以看出，隨著分辨率的提高，即r 的增大，所建模中的球形表面“顆粒感”減弱，即氣固界面更加光滑真實，并且更接近于真實球面。因此，如果需要更精細化的建模結果，提高比例尺r；反之則減小r。另外，相同參數設置下，CPU時間隨著r的增大而增加。

圖10 x,y,z三個方向所有截面孔隙率Fig.10 Porosities of all cross-sections in x,y and z directions

圖11 重復實驗得到的結冰三維孔隙率Fig.11 3-D porosities of ice by repetitive experiments

4 結 論

對動態結冰二維顯微圖像進行分割處理及定量分析，并基于此提出了結冰三維微觀孔隙結構的建模方法，得到如下結論。

(1)建模方法是可行的。建模過程中包括孔隙類球形、孔隙隨機性、孔徑分布規律性等在內的定量依據均來自二維圖像信息，建模方法的基礎可靠；建模過程中，將二維定量信息推廣至三維，遵循氣泡孔隙的隨機性、孔徑分布等，實現方式有據、特征統一；實驗部分對比驗證了所生成結冰與實驗所得結冰的二維微觀定量信息，結果表明它們具有較高的相似度。

(2) 建模方法為結冰精細化研究提供新途徑。在本文建?；A上，可以進行孔表面積、孔體積等三維孔隙空間結構特征參數的表征工作；可拓展至結冰生長模擬的微觀展示方面，提高冰形計算精度；開展小尺度結冰融化過程模擬研究，改善現有結冰融化模型。

(3)建模方法尚有完善空間。為了數學化展示結冰三維微觀形貌，建模方法必須基于一定假設，而文中假設來自于有限的研究工作，存在完善空間；二維定量信息取自遠離基底表面的穩定生長段的實驗結冰，實驗結果不受基底表面物性參數影響，結果不能直接應用于基底表面結冰特征研究，且對結冰黏附力研究無指導意義。