介觀尺度下多孔介質內水結冰相界面演化機制研究

王文松，楊英英，2，陳周林，楊晴雨，李帥華，武衛東

（1 上海理工大學能源與動力工程學院，上海 200093； 2 上海市動力工程多相流動與傳熱重點實驗室，上海 200093）

引 言

寒區凍土的自然凍融循環[1-3]、食品冷凍保鮮[4-6]、蓄冰制冷[7-8]等都涉及到水在其內部的固液相變問題，研究多孔介質內水凍結過程的傳熱傳質問題在寒區工程、生物和能源等領域中有重要意義。如在寒區建筑工程中，土壤的凍脹和融沉引起的建筑道路破壞給工程建設帶來巨大的經濟損失[9-10]。通過對土壤等多孔介質內的相變過程[11]和溫度場[12]研究，可提出凍害的預報與防治方法，緩解凍脹帶來的危害[13-15]。

目前已有學者通過宏觀實驗來研究多孔介質內的凍結特性，馬欽[16]實驗研究了不同顆粒直徑與表面溫度對多孔介質溫度分布與凍結速率的影響。但大部分的宏觀熱質規律實驗研究難以揭示其微觀的相變熱質傳遞機理。對此，一些學者通過開展孔隙尺度或介觀尺度實驗研究多孔介質內水相變熱質傳遞機理。L?voll 等[17]觀測了隨機二維多孔介質內兩相流驅替現象，驗證了滲流理論。Fen-Chong等[18]使用含水玻璃微珠來構建多孔介質，通過電容法來研究冰-水相變過程中兩相壓力變化。

同時，為了揭示多孔介質內水相變的影響因素及作用機理，一些學者開展了理論與模擬研究。Everett[19]提出毛細理論，認為毛細管在彎曲冰-水界面上產生毛細吸力。Gharedaghloo 等[20]利用孔隙尺度模擬來研究凍融循環時接觸角滯后對多孔介質含水率的影響。劉正明[21]模擬了未凍水在不同的冰-土顆粒組成的孔道微尺度模型中的遷移流動，分析了不同粒度組分下未凍水的遷移速度及凍結緣處的滲透率。Gawin 等[22]采用非平衡法模擬水-冰相變來研究多孔介質中水過冷與含水率滯后現象。Wu 等[23]認為冰晶生長的驅動力是基于熱力學的不同溫度下化學勢梯度。Lin 等[24]對納米顆粒與冰-水界面的相互作用進行了系統的分子動力學模擬，得到由于納米粒子對界面的釘扎作用，冰-水界面形成凹曲率，抵消了晶體生長的驅動力，導致冰晶生長速度變慢。李偉斌等[25]提出三維微觀孔隙結構的建模方法，基于結冰二維定量信息表征三維模型的二維孔隙定量信息，且與實驗結果相吻合。

綜上所述，目前已有學者開展了針對多孔介質凍結的相關實驗與模擬研究，但是從介于微觀與宏觀之間的介觀尺度（μm～mm 級）上對多孔介質內水凍結的相變過程進行的可視化實驗研究非常少。本文利用紅外熱像儀與微距CCD 相機對孔隙尺寸為500 μm 的多孔樹脂內水的凍結過程進行實驗研究，分析其相界面及溫度場的演變規律。并對不同邊界溫度（-15、-10、-5℃），不同孔隙間距（50～500 μm）,不同孔隙結構（圓形、方形）下的凍結過程進行模擬研究，分析不同凍結邊界溫度、孔隙間距和孔隙結構對相界面演變規律和水凍結速率的影響，研究溫度場驅動的凍結演變機制。

1 理論與方法

1.1 毛細理論模型

在微孔道凍結過程中，水相與冰相的轉變與傳熱均滿足連續性方程與導熱能量方程。

式中，ρ為密度；u為速度矢量；c為比熱容；λ為熱導率；qv為熱源，主要來自相變潛熱。

其中未凍水的流動符合不可壓縮Navier-Stokes方程。

式中，Pw為水壓；μ為動力黏度。

當孔道的孔隙半徑處于毛細尺度以下時，相界面還受到毛細壓力的影響。毛細理論可以描述水冰系統相平衡時水壓力、冰壓力與溫度之間的關系，該理論基于廣義Clausius-Clapeyron方程[26]：

其中，Pi為孔隙冰壓；L為水熱相變潛熱，取值為334.88 kJ/kg；Tm是水的凍結溫度，如圖1(a)所示。式(4)表明，凍結過程中的溫度梯度會影響凍結面上的壓力差。在實際情況中，冰的應力張力在各方向上是非均勻性的，式(4)中孔隙冰壓Pi為表面應力張量的法向分量。冰水界面上水壓力Pcl為：

圖1 冰-水相界面示意圖Fig.1 Schematic of ice-water phase interface

在毛細尺度下，由于壓力差導致相界面產生月牙形的彎曲，如圖1(b)，該壓力差由Young-Laplace方程控制[17]：

式中，γiw為冰水界面的表面張力系數，取值為33×10-3N/m；r為冰水相界面的曲率半徑。則孔隙中毛細壓力Pc為：

式中，rp為孔隙半徑，與多孔介質的孔隙結構和大小有關，因此在不同的孔隙結構中，凍結面演變也不同。聯立式（4）和式（6）得到凍結面的溫度Tf：

式（8）為Gibbs-Thomson 方程[27]，冰水界面的曲率和溫度均受該方程約束。當凍結面溫度Tf低于臨界溫度時，冰水界面的曲率半徑小于孔隙半徑，則冰層向孔隙中發展，形成向前彎曲的界面，同時相界面中心位置的平衡溫度高于界面上其他位置的平衡溫度，并且隨著相面曲率半徑的減小而增大，同時相界面溫度并不等于液態水凍結溫度，存在一定的過冷現象[28]。

1.2 單向凍結實驗研究方法

1.2.1 介觀尺度單向凍結實驗臺 本文通過搭建介觀尺度單向凍結實驗臺，研究孔隙半徑為500 μm的微孔道內水凍結過程中溫度與相界面的變化，實驗臺如圖2(a)所示。實驗臺由環境控制及觀測兩部分組成，其中環境控制部分由恒溫水槽、保溫箱與冷腔組成，可以保持環境溫度穩定和提供指定溫度邊界。觀測部分如圖2(b)所示，保溫箱上方設置觀測孔，通過伸入紅外成像儀與微距CCD 鏡頭的方式觀測溫度演變與凍結面演變。

圖2 單向凍結實驗系統圖和實物圖Fig.2 System diagram and object diagram of unidirectional freezing experiment

通過進行重復性單向凍結實驗，分別對凍結過程中的溫度場及凍結面演變進行觀測。使用紅外成像儀記錄整個凍結過程的溫度場變化，微距CCD相機拍攝凍結面變化，邊界溫度Th為-20℃，初始溫度T0為20℃。使用儀器及關鍵參數如表1所示。

表1 設備相關參數Table 1 Parameters of equipment

1.2.2 多孔樹脂微模型 本實驗設計了具有二維多孔介質結構的微模型。該模型由樹脂3D 打印制成，精度達到10 μm，結構如圖3 所示。多孔模型總尺寸為20 mm×10 mm×10 mm，多孔骨架為圓形。如圖3(b)，其中的微網格為直徑1 mm 的微型圓柱，圓柱間距為500 μm。該模型能從二維角度模擬多孔介質內水凍結過程的特征。底部基底用于支撐孔隙骨架，其高度對實驗結果無影響。由于冷板（冷源）溫度恒定，熱通量不會受到基底高度變化的影響。

圖3 被測多孔樹脂微模型結構示意圖Fig.3 Structural diagram of the tested porous resin micro model

1.3 模擬研究方法

1.3.1 物理模型 為研究不同孔隙形狀與孔隙排布對凍結過程的影響，本文設置了四種不同的孔隙結構。圖4為孔道的部分放大圖。水側位于孔隙骨架外側，其中d為直徑或邊長，l為孔間距。

圖4 四種孔隙結構尺寸圖Fig.4 Dimension diagram of four pore structures

模型的幾何尺寸與實驗對象相同，邊界條件設置及模擬取點位置如圖5所示。下側邊界設置為恒定溫度Th，提供凍結過程的冷量。根據上述的假設條件，上側邊界為絕熱條件，兩側邊界為對稱條件。多孔介質邊界為潤濕壁條件。流體初始溫度為5℃，流動狀態為靜止。圖中標記點位置為后續實驗和模擬所選測點。

圖5 孔隙結構物理模型邊界條件Fig.5 Boundary conditions of physical model of pore structure

1.3.2 數學模型 在平面坐標系中，該模型的控制方程為：

二維連續性方程

流體物性來自REFPROP 軟件的材料庫，各物性隨溫度的變化而變化。部分參數如表2所示。

表2 模型物性參數Table 2 Physical parameters of model

1.3.3 模型驗證 本文使用的模擬軟件為Comsol。為驗證模擬結果的準確性，將模擬結果與實驗測試結果進行了對比。模擬條件如下：初始溫度20℃，邊界溫度-20℃。模擬結果與實驗結果對比如圖6所示。

由圖6 可知，模擬結果與實驗結果有較好的吻合性，絕對誤差平均值為1.4℃。實驗值與模擬值的主要偏差集中在凍結初始階段，原因是在實驗過程中低溫恒溫槽輸出的低溫冷凍液到達冷板時會發生能量損耗，實際邊界溫度略低于-20℃，因此導致在凍結初始階段降溫速率與理論模擬有一定的差別。

圖6 凍結過程中相同位置溫度模擬結果與實驗結果對比Fig.6 Comparison of temperature simulation results and experimental results at the same position during freezing

2 結果與討論

2.1 單向凍結實驗結果分析

2.1.1 凍結面演變規律 根據毛細凍脹理論，接觸角不同與微孔道結構變化會使相變過程受到不同程度的毛細作用力。根據Young-Laplace 方程理論計算，可以得到凍結面在經過孔道時，毛細壓力的變化。圖7（a）為孔道中孔隙半徑分布示意圖；圖7（b）展示了在多孔樹脂的孔隙結構中，隨著凍結面推進，孔隙半徑rp與最大毛細壓力Pmax在相對位置h上的變化。相對位置h位于孔隙中軸線處?？梢钥闯?，圓形孔道中，毛細壓力在孔隙半徑最小處達到毛細壓力最大值，之后隨著孔隙半徑增大而減小。

圖7 微孔道中孔隙半徑和毛細壓力隨相對位置的變化Fig.7 Variation of pore radius and capillary pressure with relative position in microporous channel

圓形微孔道中凍結面的相對曲率分析如圖8(a)所示。在凍結面與壁的接觸點處，凍結面受到冰壁表面張力和壁表面張力的影響，形成彎曲相界面現象，其中水與壁面的接觸角為θlw。如圖8(b)所示，在凍結初期，凍結面進入孔隙半徑逐漸減少的孔道時，毛細壓力Pc逐漸變大，邊界毛細作用逐漸明顯。在凍結面上共同受到水壓Pw與冰壓Pi的影響，由于毛細壓力隨凍結面在孔道中發展而變化，因此水壓與冰壓差也隨之變化，如圖8(c)、(d)所示。

圖8 凍結面變化示意圖Fig.8 Schematic of frozen surface change

為量化界面曲率的變化，引入相對曲率概念研究相界面經過孔道時的變化。首先在相界面提取一定量的數據點，根據足夠的數據點進行函數擬合，函數根據圖9(a)中的坐標軸定義，之后把擬合出的函數代入式（14），得到相界面的相對曲率k。

如圖9(b)所示為凍結面相對曲率隨時間的變化。實線為根據實驗中拍攝的凍結面所計算的曲率值，虛線為數值計算模擬值?？梢钥吹?，在凍結面經過一個圓形孔道時，相對曲率呈先上升后下降再上升最后下降的趨勢。相界面從孔道發展結束，準備進入下一排孔道時，由于缺少孔壁約束，相界面逐漸平緩，相對曲率逐漸減小。進入孔道以后，毛細作用逐漸明顯，表面張力從水、冰、壁面接觸點向孔道中軸處影響，導致相對曲率逐漸增加。由于圓形孔道孔隙半徑減小，導致毛細壓力逐漸變大，冰壓與水壓差變大，凍結面會有突起的趨勢（相對于冰側），所以凍結面會從凹形向凸形發展，凍結面相對曲率會先下降至0（點b，平直界面）再上升至點c。最后，隨著凍結面繼續前移，孔隙半徑逐漸變大，凍結面失去孔道的限制逐漸平緩，相對曲率逐漸下降。

圖9 凍結面相對曲率變化規律Fig.9 Change law of relative curvature of freezing surface

2.1.2 溫度場演變規律 圖10 顯示了在邊界溫度Tc為-20℃的情況下，隨著時間的變化，溫度場在孔道間的演變，其中白色部分為多孔樹脂固體骨架，彩色部分為孔道介質溫度。在50 s 時，溫度場沿著孔道向前推移至第三排孔道，至90 s 孔道溫度到達7℃左右，在190 s 到達2℃左右， 在240 s 時，孔道溫度到達0℃，到達310 s時，孔道溫度達到-6℃。溫度場下降速率呈先快后慢的趨勢，與溫度場發展的溫度梯度有關。

圖10 凍結過程中微孔道溫度變化Fig.10 Temperature change of microporous channel during freezing

從圖11 微孔道內測點溫度隨時間變化也可以得到相同結論。從圖中可以看出，環境溫度變化在1℃范圍內，說明保溫箱內溫度基本不變，可視為實驗在環境溫度不變的條件下進行?？椎纼葴囟仁艿綔囟忍荻鹊挠绊懗氏妊杆傧陆?、后緩慢下降的趨勢，分界點在150 s 左右。位置1 的溫度下降不但比位置3 快，300 s 時溫度同樣更低，說明微孔道內在300 s已經形成基本穩定的溫度場，溫度變化逐漸減小。所以如果要保持溫度下降速度，需要改變邊界溫度，保證溫度場形成一定的溫度梯度。

圖11 單向凍結實驗測點溫度隨時間的變化Fig.11 Temperature variation of measuring points with time in unidirectional freezing experiment

圖10(d)中，相界面與0℃等溫線有一定的距離，這是因為在多孔介質中，毛細壓力會造成過冷的現象，因此在到達相變溫度時并不會立即凍結，導致在0℃等溫線與相界面之間出現過冷帶。

為研究過冷帶在孔道間的遷移，同時分析過冷帶所處的溫度場。圖12為單個孔道中的溫度演變，其中紅色代表高溫，藍色代表低溫，黃色代表溫度范圍為-1～0℃。選取時間為228、240、253、266 與285 s（相對時間為0、12、25、38、57 s）?？梢詮膱D中看出，由于交錯排布的孔道結構，導致冷量從孔道兩側向中間匯聚，過冷帶的形狀呈月牙狀。在進入孔道后，由于孔隙半徑的減小，導致溫度梯度在孔隙間迅速增大，因此過冷帶在孔隙間曲率逐漸變小，并在253 s時可以清晰看到過冷帶邊界平直。然后隨著溫度場發展，過冷帶進入下一排孔道。

圖12 單向凍結實驗單孔道溫度場演變Fig.12 Evolution of temperature field in single channel of unidirectional freezing experiment

2.2 單向凍結模擬結果分析

2.2.1 邊界溫度對溫度場及凍結速率的影響 通過模擬可以得到邊界溫度為-15℃的溫度場及對應的相界面云圖，圖13 所示為300 s 時的溫度場及相界面云圖?？梢钥闯?，300 s 時溫度場發展至第六排，且前排溫度梯度比后排高。相界面發展至第四排。

圖13 在300 s時刻的溫度場及相界面云圖（邊界溫度-15℃）Fig.13 The temperature field and the phase interface cloud Contour at 300 s (boundary temperature -15℃)

為分析不同邊界溫度對溫度場演變的影響，取一測試點（圖5）進行研究，圖14 展示了該點在不同邊界溫度時的凍結過程中溫度變化曲線。從圖中可以看出，在前200 s，由于溫度梯度大，溫度迅速下降并接近凍結溫度，在200 s后，溫度梯度趨于平緩，溫度變化速率逐漸減小。邊界溫度越低，溫度下降越快，不同邊界溫度主要影響發生在凍結開始階段的溫降速度與凍結晚期的穩定溫度。不同邊界溫度條件下，該位置到達相變溫度所用的時間也不同。邊界溫度為-5℃的條件下，該點在52 s 時達到相變溫度，邊界溫度為-10℃的條件下需要22 s，而邊界溫度為-15℃的條件下，只需要14 s。

圖14 不同邊界溫度在相同位置的溫度變化Fig.14 Temperature change of different boundary temperatures at the same position

測試點處溫度梯度隨時間變化如圖15 所示?？梢钥吹皆趦鼋Y至200 s左右，溫度梯度發展至最大值，因此在0～200 s 的區間中，測點溫度下降速率最快。而在200 s 之后，隨著溫度場向y方向的發展，溫度梯度逐漸下降，因此在200 s后測點溫度趨于穩定。同時邊界溫度為-5℃的最高溫度梯度比邊界溫度-15℃的最高溫度梯度小，說明溫度越低的邊界溫度帶來的溫度梯度越高，發展至穩定的溫度場所需的時間也越長。

圖15 溫度梯度隨時間的變化Fig.15 Temperature gradient versus time

圖16展示了不同邊界溫度下的平均凍結速率，邊界溫度為-10℃的凍結速率比邊界溫度為-5℃的凍結速率快36.37%，邊界溫度為-15℃的凍結速率比邊界溫度為-10℃的凍結速率快91%，證明邊界溫度可以顯著提高凍結速率，并且邊界溫度越低凍結速率加快越明顯。

圖16 不同邊界溫度下的平均凍結速率Fig.16 Histogram of average freezing rate at different boundary temperatures

2.2.2 孔隙間距對凍結速率的影響 為探究孔隙間距對溫度場演變的影響機理，從溫度梯度角度對邊界溫度為-15℃時，孔隙間距為500、200、100 和50 μm 的圓形交錯孔隙結構進行分析。圖17 為多孔介質中測點溫度隨時間的變化曲線?？梢钥吹娇紫堕g距為500 μm的孔隙結構到達相變溫度最快，需要14 s；孔隙間距為50 μm 的孔隙結構到達相變溫度最慢，需要31 s?？梢宰C明孔隙間距越小，溫度變化速率越慢。

圖17 不同孔隙間距中溫度隨時間的變化Fig.17 Temperature variation with time in different pore spacing

圖18 顯示了不同孔隙間距的孔隙結構中溫度梯度隨相對位置的變化。其中相對位置是指孔道中軸線位置，0 處位于孔道最窄處?？梢悦黠@看出，孔隙間距為500 μm 的溫度梯度變化平穩，最高到1103 K/m，且變化范圍更寬達到1.4 mm?？紫堕g距為50 μm的溫度梯度變化劇烈，最高達到1389 K/m，變化范圍只在0.58 mm 內產生變化，可以量化證明溫度梯度在孔隙中的變化。上述說明了不同孔隙間距對溫度場發展產生影響的原因，由于在孔隙間距最小處發生溫度梯度突躍，會阻礙溫度場發展速度，并且突躍的程度與阻礙溫度場呈正相關，因此孔隙間距越小，溫度場發展越慢。

圖18 單個孔道溫度梯度隨相對位置的變化Fig.18 Variation of temperature gradient of single channel with relative position

平均凍結速率如圖19所示。由圖可知，孔隙間距越小，平均凍結速率越慢。其原因是由于在該孔隙結構中，凍結過程主要由溫度場驅動，溫度場發展越快的孔隙結構，凍結發展速度也越快。結果表明，孔隙間距越小的孔隙結構，其凍結速率越慢。

圖19 不同孔隙間距下的平均凍結速率Fig.19 Histogram of average freezing rate under different pore spacing

2.2.3 孔隙結構對凍結速率的影響 圖20 為四種孔隙結構的平均凍結速率，可以看到凍結速率由快到慢分別是圓形直排、圓形交錯、方形直排與方形交錯。為進一步證明在多孔介質凍結過程中凍結場主要由溫度場進行驅動演變，通過熱流線及等溫線研究四種孔隙結構對溫度場演變影響的機理。圖21為四種孔隙結構的多孔介質在同一時間、位置的溫度場孔道局部圖，其中紅色代表高溫，藍色代表低溫，黑線表示熱流線及等溫線，相鄰兩條等溫線相差0.25 K。

圖20 不同孔隙結構下的平均凍結速率Fig.20 Histogram of average freezing rate under different pore structures

圖21 不同孔隙結構溫度場孔道局部圖Fig.21 Partial diagram of pore channel for temperature field of different pore structures

通過對比圓形孔道與方形孔道，可以發現方形孔道由于受到形狀的限制，當圓形孔道等溫線已經進入孔道時，方形孔道等溫線仍在孔道外的壁面上，因此，圓形孔道比方形孔道更利于傳熱。對比圖中直排與交錯排列的孔隙排布，交錯孔道的熱流通道（整體熱流線）彎曲導致等溫線的彎曲，這意味著一部分冷量提前向后部進行傳熱，導致傳熱面積的減小，減緩了凍結面進入孔道的時間，因此，直排孔道中溫度場演變速率比交錯孔道的溫度場演變速率更快。結果表明，孔隙骨架形狀越光滑，熱流通道越趨于平直，溫度演變越快。

3 結 論

本文實驗研究了多孔介質單向凍結過程中相界面與溫度場的演變，并通過數值模擬方法，研究了不同邊界溫度、孔隙間距及孔隙結構對多孔介質凍結過程及凍結速率的影響，得出如下主要結論。

（1）凍結面在經過孔道時，孔隙半徑先減小后變大，毛細壓力隨之先增大后減小，因此相對曲率先減小后增大。凍結面從凹形向凸形變化，在孔道中軸處相對曲率下降至0，凍結面平緩。

（2）邊界溫度相同時，孔道內同一位置處的溫度梯度先快速增加，隨著凍結的進行逐漸減小，導致該點溫度下降呈現先快后慢的趨勢，最后趨于穩定。邊界溫度越低，形成的溫度梯度越高，進而提高其凍結速率。不同邊界溫度，相同位置處，達到最高溫度梯度的時間相同?？椎纼人Y冰過程存在過冷現象，0℃等溫線和相界面之間存在過冷帶。

（3）溫度梯度在孔道最窄處會發生突躍，阻礙溫度場的發展，孔隙間距越小，突躍越劇烈，凍結速率越低?？紫督Y構對凍結速率的影響與孔隙骨架形狀和排布方式相關，孔隙骨架的形狀越光滑，凍結面進入孔道越快，排布方式越規律，溫度演變越快。