動態(tài)結(jié)冰孔隙結(jié)構(gòu)三維建模方法

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

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

引 言

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

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

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

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

1 動態(tài)結(jié)冰孔隙分布的基本假設

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

圖1 具有氣泡的動態(tài)結(jié)冰Fig.1 Dynamic icing with pores

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

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

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

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

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

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

2 動態(tài)結(jié)冰微觀結(jié)構(gòu)的三維建模方法

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

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

2.1 三維微觀結(jié)構(gòu)的區(qū)域

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

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

2.2 區(qū)域內(nèi)孔隙數(shù)量

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

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

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

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

2.3 孔隙位置

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

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

2.4 孔徑分布

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

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

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

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

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

2.5 三維微觀孔隙結(jié)構(gòu)

更新零矩陣T，生成動態(tài)結(jié)冰三維微觀孔隙結(jié)構(gòu)，其對應值為1。

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

3 建模結(jié)果及可靠性分析

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

3.1 二維定量信息

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

3.2 三維微觀結(jié)構(gòu)生成

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

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

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

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

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

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

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

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

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

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

圖5 所生成1 cm×1 cm×1 cm結(jié)冰的圖像顯示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 二維信息定量對比 為了直觀顯示所生成的三維結(jié)構(gòu)內(nèi)部形態(tài)，圖9 對比展示了結(jié)冰孔隙實驗測量結(jié)果及所生成結(jié)冰二維截面。圖中可以看出生成結(jié)冰的內(nèi)部孔隙分布與實驗較為相似，表現(xiàn)在孔隙大小、分布、隨機性等方面。

實驗圖片的分割結(jié)果中，二維氣泡的孔隙率約為0.0512，為了對比驗證本文建模方法的準確度，對所生成結(jié)冰三維微觀結(jié)構(gòu)的所有x,y,z 方向的截面進行孔隙率計算，得到的結(jié)果以柱狀圖展示于圖10中?？梢钥闯鏊鼈兊闹抵饕性?.05 左右，變化范圍有限，并且它們的平均值為0.0511，與二維實驗值相近，這說明三維建模得到的孔隙率結(jié)果與二維定量結(jié)果有較高的一致性。

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

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

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

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

3.4 不同分辨率下的建模結(jié)果

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

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

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

4 結(jié) 論

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

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

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

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