壁面電荷對銅表面冰黏附的影響研究

蔡文豪，許雄文，2

（1 華南理工大學電力學院，廣東 廣州 510640； 2 廣東省能源高效清潔利用重點實驗室，廣東 廣州 510640）

引 言

過冷水式動態冰蓄冷技術制備具有相變潛熱的冰漿作為儲能介質，具有高能量密度的特點，且制備系統簡單，傳熱效率高[1]。然而該方法面臨的一個問題就是所形成的冰晶極易黏附在低溫換熱表面[2]，降低系統的換熱效率，甚至引發堵塞流道的問題。目前常采用的解決方案是提高過冷水的流速[3-4]，利用流體擾動產生的分離效果剝離過冷壁面上的黏附冰，但大大增加了水泵功耗。減小冰與壁面之間的黏附強度[5-7]是降低過冷水流速的有效途徑，因此有必要分析固體表面冰的脫附機制，找到減小冰黏附強度的有效方法。另外，冰黏附在固體表面還會對很多工業運行安全造成威脅。如風機葉片和光伏面板上附冰會降低可再生能源基礎設施的性能[8]，熱泵翅片結霜會影響制冷設備的穩定運行[9]，機翼附冰會威脅航空安全[10]，等等。因此減小冰與壁面之間的黏附強度研究具有很好的科學和工程意義。

冰在自然狀態下最常見的是六角形形態，根據Bernal-Fowler 冰規則[11]，每個水分子必須接受和提供兩個氫鍵，才能維持規則的冰晶體結構，然而對處于界面區域的水分子來說，這顯然是不可能的，其結構必然會由于兩側的受力不同以及固體壁面無法提供氫鍵而發生重構，形成定向非隨機無序的結構，使該層中水分子具有可移動性，最終導致該層黏度與液態水相當，稱為準液體層[12]。準液體層已被一系列實驗證明存在，并在剪切和剝離過程中具有潤滑作用。為了降低表面黏附，人們利用這一特性，將與水不互溶的潤滑劑注入多孔表面(slippery liquid-infused porous surface，SLIPS)[13]或保持在粗糙納米結構表面涂層[14-16]中，使固體表面結冰時，潤滑劑能填充到冰和固體表面之間，降低冰的黏附強度，從而實現冰的自然力（例如風、重力和振動）去除[17]。一些設計力求獲得真正的準液體層而非使用潤滑劑。它們采用高吸水性涂層[18]、吸濕性聚合物[19]、水合性聚合物[20]以及聚電解質涂層[21]，利用高的吸水性和水合作用在壁面區域形成并維持較厚的準液體層，以減小冰的黏附強度。相當于通過表面處理來保持比原有狀態下更厚的準液體層來實現低黏附強度的目的，這種采用潤滑層來促進黏附冰脫附的策略在實際應用中具有很好的應用前景[22-23]。

在動態蓄冰的換熱過程中，流體的剝蝕力較小，重力和振動基本可忽略，對于冰的脫附要求更高。因此，希望通過其他物理手段來增加準液體層厚度并減小冰的黏附力。固體表面荷電，可以改變近壁面水分子的偶極方向，從而改變水的結冰過程。在帶電Pt（qPt=0.12 e，1e=1.6×10-19C）納米通道中，吸附在帶電表面的水分子之間會形成一個二維的氫鍵網絡，干擾水分子層之間的氫鍵形成，最終導致壁面附近無法結冰，且冰在接近表面時生長速度減慢[24]。在帶電石墨烯上的水覆蓋層會隨著碳原子攜帶電荷量（qC=0～0.18 e）的增加依次經歷冰到液體和液體到冰的轉變[25]。因此壁面荷電有望增加準液體層厚度。由于準液體層厚度在納米尺度，為了驗證這個想法，需要進行低于冰點的水分子體系的分子動力學模擬實驗。

分子動力學模擬憑借其原子級別的分辨率在準液體層的模擬中占有一席之地[12,26]。對于黏附冰脫附的行為，分子動力學模擬包含了模擬空間中所有的力學關系，可以方便地導出固體壁面對冰塊的黏附力。Xiao 等[27]在原子尺度上完成了冰黏附力學的分子動力學研究，并且給定準液體層厚度，研究了硅表面和石墨烯表面準液體層的存在對冰黏附強度和剪切強度的潤滑作用。然而人為給定的液層厚度無法反映準液體層對外界環境的依賴性。Afshar等[28-29]改進了上述不足，他們使用全原子分子動力學模擬，分析了平衡后冰塊的靜態結構，關注了不同溫度下準液體層的厚度變化，并研究了表面粗糙度和溫度對石墨襯底上的冰黏附強度和剪切強度的影響，相對全面地在原子層面解釋了冰的脫附機理。

本文擬在兩銅壁面之間構建納米尺寸的冰立方模型，采用銅板原子靜態和脈沖荷電的方式來改變原有的準液體層厚度，通過對水分子施加垂直于壁面且大小一定的拉力，使得冰塊恰好脫離壁面，以測定其黏附強度。本研究擬量化同一溫度不同壁面荷電量時的準液體層厚度及其對應的冰黏附強度，以期提出減小冰黏附強度的有效方法。

1 模型和模擬

1.1 水分子-銅壁面體系模型

1.1.1 模擬參數 本研究采用分子動力學軟件LAMMPS[30]進行納米尺寸冰塊脫附的模擬，采用數據可視化軟件OVITO[31]觀察冰塊的平衡情況。冰塊的建模采用TIP4P/ICE[32]剛性非極化水模型，因為該模型產生的六角形冰的熔點為272.2 K，密度為0.909 g/cm3（250.0 K，1 bar，1 bar=0.1 MPa），能夠很好還原水的固液相圖，其模型參數見表1。在模擬過程中使用 SHAKE 算法[33]固定 H—O 鍵的鍵長以及H—O—H 的鍵角，以保持水分子的結構。采用晶格常數為3.615的面心立方單元對具有1 nm 厚度的銅板進行建模，得到由六層銅原子組成的銅板。銅原子的相對位置在之后的模擬中不進行更新，這相當于忽略銅原子的熱運動。這一假設在冰黏附的分子模擬中被廣泛采用[27-29]，在研究水在固體表面上的接觸角時也有采用[34-35]。由于本文并不對銅原子的相對位置進行更新，所以銅原子相互作用勢的選擇并不是唯一的，甚至可以說是任意的。為滿足所有種類的原子之間都具有對相互作用類型的要求，銅原子的相互作用采用嵌入原子模型（eam）。水分子與水分子之間和水分子與銅壁面之間的相互作用采用Lennard-Jones12-6 勢函數和庫侖相互作用來描述，設置L-J 勢的截斷半徑為10 ?（1?=1×10-10m），庫侖相互作用的截斷距離為8.5 ?，超過這一距離的庫侖相互作用采用PPPM/TIP4P 在倒數空間中計算。

表1 水分子模型參數（TIP4P/ICE）Table 1 Model parameters of water molecule（TIP4P/ICE）

式中，E為兩粒子之間的勢能；ε為勢阱深度，描述了兩個粒子之間的相互作用強度；σ為兩個分子處于平衡位置時的距離；C是能量轉換常數；qi和qj是i原子和j原子上的電荷；r是兩個分子之間的距離；rc為對應的截斷距離。

水在紫銅表面的靜態接觸角θ為86.4°[36]，根據靜態接觸角與勢阱深度的線性關系θ=178.57°-497.41°×εO-Cu[35]，本文選用εO-Cu=0.1853。Cu和O之間的平衡距離參數根據Lorentz-Berthelot 混合規則[37]確定。

銅原子之間平衡距離σCu-Cu=2.3377 ?[38]，得到σO-Cu=2.7523 ??？紤]到重力在納米尺度上遠不如相互作用勢重要，因此本文忽略了重力的影響。模擬體系中的相互作用參數列于表2，其單位選用LAMMPS中的 real 單位制。

表2 各組分之間的相互作用參數Table 2 Interaction parameters between components

1.1.2 物理模型 冰塊由六角形冰晶胞在x、y、z三個方向上復制得到，包含2160 個水分子，如圖1（a）所示。擴大模擬盒子并沿冰晶體基面{0001}[39]引入銅板，如圖1（b）所示。融化界面區域的冰塊，形成1.5 nm 厚的液層，以保證在不同工況下進行結冰模擬時形成對應的準液體層，如圖1（c）所示。為保證荷電工況下體系呈現出電中性，在冰塊上方引入一塊可以上下活動的銅板，該銅板帶有與下銅板等量相反的電荷，如圖1（d）所示。整個模擬空間尺寸為4.13 nm×3.88 nm×11.50 nm，x和y方向采用周期性邊界。由于本文中z方向的周期性并沒有意義且為了節省計算資源，將z方向設為固定邊界條件。同時將模擬盒子z方向的周期性圖像變為空(體積因子設置為 3)[40]，以滿足長程庫侖力求解器要求的周期性邊界條件。

圖1 物理模型建模過程Fig.1 Modeling process of physical model

1.2 平衡模擬和冰水判定

本文中模擬分為平衡和脫附兩個階段。在平衡階段，以上述準備好的物理模型為初始結構，保持溫度T=255 K，銅板分別加載不同電荷量（Q=0 e/nm2、脈 沖 交 變Qperiod=±0.1123 e/nm2、靜 電Qstatic=±0.1123 e/nm2）進行三個相互獨立的結冰模擬。本文中壁面帶電量與文獻[40]中的數值在同一量級。當銅壁面荷電時，電荷平均地加載到與冰塊直接接觸的那一層銅原子上[40-41]，如圖2所示。銅板靜態荷電和脈沖荷電時，下銅板帶電荷銅原子（數量N=242）的帶電量隨時間變化如圖3（a）、（b）所示，上銅板帶電荷銅原子時刻保持帶等量的相反電荷。模擬過程中，冰塊采用Nose-Hoover 方法來控制溫度恒定在指定溫度，模擬時間步長為1 fs，耦合時間常數為

圖2 銅板荷電時的模型圖（紅色銅原子帶正電荷，藍色銅原子帶負電荷）Fig.2 Model diagram of charged copper plate (copper atoms in red are positively charged, and copper atoms in blue are negatively charged)

圖3 下銅板原子的荷電量時變關系Fig.3 The time-varying relationship of the charge of the lower copper plate atoms

0.1 ps，每1 ps進行一次數據采集。

采用鍵序參數Ql來區分水分子是處于“固態”還是“液態”[42-43]。根據文獻[44-45]的方法，計算了溫度T=255 K 下六角形冰和液態水的鍵序參數Q6，其概率密度函數(probability density function，PDF)如圖4所示。同溫度下冰和水的Q6曲線的交點作為區分水分子狀態的分界線，高于該值水分子處于“固態”，否則處于“液態”?？梢缘玫皆?55 K 時的代表分界線的Q6值為0.3992。如果冰分子的數目隨時間沒有明顯的上升，則可判定冰塊和壁面體系已經達到平衡，界面區域的準液體層形成，冰塊穩定地附著在銅板上，可以進行準液體層的分析以及脫附過程的模擬。

圖4 冰和過冷水Q6的概率密度函數圖Fig.4 Probability density function of Q6 of ice and water

1.3 脫附模擬和黏附強度的測量

在脫附階段，為保證體系電荷守恒而引入的上銅板可以上下活動，下銅板固定，給所有水分子施加一個沿銅板法線方向的拉力，使其具有遠離下銅板的趨勢，如圖5 所示。通常，拉力的大小從零開始，隨時間步長線性增加，使基底和冰塊之間相互作用的z方向分量表現為吸引力，直到拉力足夠大，使冰塊脫離壁面[29,46]。記錄過程中吸引力大小與相應拉力的值，再利用冰-基體接觸面積進行換算，將吸引力和拉力的值轉化為界面法向應力和外加法向應力。

圖5 外加拉力使冰塊脫附示意圖Fig.5 Schematic diagram of applying tension to make ice block desorb

式中，F是基底對冰塊的吸附力或拉力附加給冰塊的拉力；A是冰-基底接觸面積。

本文采用給所有水分子施加恒定拉力的方法進行黏附強度測量。在每種工況下選擇10 個不同時刻的構型，測定每個構型中冰塊在2.5 ns 內脫附的最小拉力，再將10個構型中的最小脫離應力作為該工況下的黏附強度。

2 結果與分析

本文討論了溫度T=255 K 時的銅板脈沖荷電以及靜電荷電的工況，同時以壁面不帶電荷時冰和過冷水的黏附強度為參照，表征了黏附強度降低的幅度。過冷水由冰塊升溫至350 K 融化后，再降溫到目標溫度得到。

2.1 冰塊生長和準液體層形成

圖6（a）顯示了平衡階段被判定為冰的水分子數目變化。相對于冰塊生長的平衡，大塊水的平衡要快的多，因而平衡時間較短。在溫度T=255 K 且壁面不帶電荷（Q=0 e/nm2）時，組成冰塊的水分子數目在上升到一定程度后，會在一定的范圍內波動，這歸因于分子的熱運動。在溫度T=255 K 且銅板帶脈沖電荷Qperiod=±0.1123 e/nm2時，尤其是平衡時間t=30.0～40.0 ns這一區間，冰分子數目的變化呈現出明顯的鋸齒狀，即在銅板電荷產生的電場方向轉變時，會在一定程度上打亂規則的六角形冰結構，而在隨后的時間里冰塊繼續生長，直到下一次電場方向轉變的時刻再次被破壞。溫度T=255 K 且壁面帶靜電Qstatic=±0.1123 e/nm2時，并沒有出現冰分子數目呈現鋸齒狀的現象?？梢酝茢?，在冰分子數目不再明顯增長時，體系已經處于平衡狀態。圖6（b）展示了一系列不同時刻的模型快照，并在圖6（a）中的相應位置用圓圈數字進行了標記。對應冰分子數量的增長，從模型快照中很容易觀察到冰塊的生長，同時也可以觀察到冰塊主體部分規則的冰雙層結構以及界面區域的準液體層。準液體層的分析放在后面進行，首先研究不同工況下的黏附強度。

圖6 類冰水分子數目隨時間的變化及對應時刻的模型快照（黃銅色為銅原子，紅色為類液水分子，藍色為類冰水分子，白色為氫原子）Fig.6 Time-dependent changes in the number of ice-like water molecules and model snapshots at corresponding moments (brass are copper atoms, red are “liquid” water molecules, blue are “solid” water molecules, and white are hydrogen atoms)

2.2 冰的黏附強度

在大塊的水中結冰需要漫長的形核過程[47]，這提供了充分的時間來測量過冷水的黏附強度。因此，本文對4 種工況均進行了10 次初始構型不同（體系平衡后取連續的10.0 ns，每隔1.0 ns 選取一個構型）的相互獨立的脫附模擬。

以T=255 K，Qperiod=±0.1123 e/nm2時為例，取t=33.5 ns 的構型為初始構型。脫離應力的測量需要不斷更改給定的拉力大小，根據冰塊在上一測試拉力下是否脫附對拉力進行調整，每次調整的拉力大小間隔約使外加法向應力增加或減少1.0 MPa。圖7記錄了在測定該構型脫離應力過程中，其中6種不同拉力下冰塊的質心（center of mass，COM）縱坐標隨時間的變化。下銅板的荷電情況也標在圖7 中?？梢杂^察到，脫附模擬過程中電荷產生的電場方向發生了兩次變化。冰塊質心坐標的迅速升高標志著冰塊從銅板上剝離。在外加法向應力σpull=135.4和140.3 MPa 時，冰塊質心保持在同一水平，沒有發生脫附的現象。將外加法向應力增加到141.2 MPa時，冰塊在拉力作用2.4 ns 后發生脫附。然而將外加法向應力增加到143.1 MPa 時，冰塊卻意外地沒有脫附。當進一步將外加法向應力增加到144.1 MPa 和147.0 MPa 時，冰塊脫附的時間點有所提前，且在較大外力作用時具有更大的提前幅度，這是因為外加法向應力增大，脫附概率增加，可以在更短的時間內實現脫附。綜上，在T=255 K、Qperiod=±0.1123 e/nm2的工況下，t=33.5 ns 時刻的脫離應力為141.2 MPa。

圖7 T=255 K，Qperiod=±0.1123 e/nm2，t=33.5 ns時，不同外加法向應力下的脫附模擬中冰塊質心的位置變化Fig.7 The position change of the ice mass center in the stripping simulation under different applied normal stress when T=255 K, Qperiod=±0.1123 e/nm2,t=33.5 ns

該 工 況 下（T=255 K，Qperiod=±0.1123 e/nm2）t=33.5～42.5 ns 時間段內的10 個構型（每1.0 ns 選取一個構型）的脫離應力測定結果如圖8 所示。在樣本范圍內，t=35.5 ns時的脫離應力最小，為130.6 MPa。綜上，T=255 K、Qperiod=±0.1123 e/nm2的工況下，冰的黏附強度為130.6 MPa。

圖8 T=255 K，Qperiod=±0.1123 e/nm2，t=33.5～42.5 ns時，各個冰塊構型在不同外加法向應力下的脫附情況Fig.8 Stripping of ice under different applied normal stresses when T=255 K，Qperiod=±0.1123 e/nm2，t=33.5～42.5 ns

圖9(a)記錄了各種工況下10次獨立脫附模擬得到的脫離應力。以這10 個脫離應力最小值作為這種工況下的黏附強度，如圖9(b)所示。

圖9(b)顯示，在255 K 時，Q=0 e/nm2的工況下的冰塊和銅板之間的黏附強度145.1 MPa 與過冷水與銅板之間的黏附強度99.6 MPa 之差決定了黏附強度可以減小的最大限度約為45.5 MPa。以Q=0 e/m2工況下，黏附強度可減少的最大限度為參照，Qperiod=±0.1123 e/nm2和Qstatic=±0.1123 e/nm2時的黏附強度分別為130.6 和148.0 MPa，降低的幅度分別為31.9%和-6.4%。這表明電荷的存在會增強壁面和水之間的庫侖相互作用，從而增加冰的黏附強度，然而交變電場的存在卻使冰塊的黏附強度降低。造成這一現象的原因是脫附模擬冰塊初始結構的差異，即平衡階段對應時刻的冰塊的準液體層厚度差異。

圖9 不同工況下的脫離應力以及黏附強度Fig.9 Detaching stress and adhesion strength under different conditions

圖10展現了平衡階段平衡后的冰塊靜態結構。沿銅表面法線方向上以2 ? 的厚度進行分層，將同一層內水分子的鍵序參數Q6求取平均值，得到Q6(Z)，以判斷該層水分子處于固態還是液態，其中鍵序參數Q6每1 ps 采集一次，采集的范圍是與脫附模擬取樣時間段對應的10.0 ns。圖中虛線與實線交點的橫坐標即為冰水界面的坐標，盡管實際上的冰水界面有一個過渡的區域。Q6(Z)的分布再次表明，在銅表面一定高度范圍內水分子的晶體狀結構遭到破壞，存在水分子無序排列的準液體層。由于各種工況下銅板最上層原子的坐標（Z=-1.8075 ?）相同，因此冰水界面的坐標直接反映了該工況下的準液體層厚度。這與期望的一樣，準液體層的厚度正是隨著冰塊黏附強度的降低而增厚，因此壁面加載脈沖電荷是增加準液體層厚度的有效方式，從而成為一種減小冰黏附強度的有效手段。

圖10 體系平衡時的冰塊靜態結構（虛線為區分冰和水的閾值，虛線之上為冰，虛線之下為水）Fig.10 The static structure of the ice cube when the system is in equilibrium (the dotted line is the threshold for distinguishing ice and water； above the dotted line is ice, and below the dotted line is water)

3 結 論

本文使用全原子模型的分子動力學模擬研究了光滑銅表面上黏附冰的脫附行為，采用序參量進行水分子的冰水判定，以冰塊質心位置和脫離應力大小對冰塊的脫附過程進行了定量表征，量化了255 K 時不同壁面條件下冰塊的黏附強度，得到如下結論。

在壁面施加靜態電荷無法降低冰的黏附強度，但壁面施加脈沖電荷可降低冰在壁面的黏附強度。壁面靜態荷電Qstatic=±0.1123 e/nm2對準液體層厚度影響不大，由于壁面電荷的存在增大了壁面與水分子之間的庫倫相互作用，導致冰的黏附強度由壁面不帶電荷時的145.1 MPa增大到148.0 MPa。壁面脈沖荷電Qperiod=±0.1123 e/nm2會使冰塊與固體壁面之間的準液體層加厚，從而使冰的黏附強度降低到130.6 MPa。以過冷水與銅板之間的黏附強度99.6 MPa 為參照，降低幅度達到31.9%。本文從原子尺度上揭示了準液體層在冰脫附過程中起到的關鍵作用，證明了壁面脈沖荷電是一種有效促進壁面上黏附冰脫附的方法。

符 號 說 明

A——冰-基底接觸面積，nm2

C——能量轉換常數

E——兩粒子之間的勢能，kcal/mol

F——基底對冰塊的吸附力或外加給冰塊的拉力，N

lOM——四點TIP4P/ICE 水模型中無質量的電荷點與O原子的距離，?

MO,MH——分別為氧原子和氫原子的分子量，g/mol

Ql——鍵序參數

Qperiod,Qstatic——分別為脈沖荷電的電荷密度和靜態荷電的電荷密度，e/nm2

qPt,qC——分別為單個Pt 原子和單個C 原子的荷電量，e

qO,qH——分別為氧原子和氫原子上的電荷，e

qi,qj——i原子和j原子上的電荷，e

r——兩個分子之間的距離，?

rc——對應的截斷距離，?

r0——四點TIP4P/ICE水模型中O—H鍵鍵長，?

T——熱力學溫度，K

t——平衡時間，ns

ε——勢阱深度，kcal/mol

θ0——四點TIP4P/ICE 水模型中H—O—H 鍵角，（°）

σ——兩個分子處于平衡位置時的距離，?

下角標

period——脈沖荷電

static——靜態荷電