蔡文豪，許雄文，2

（1 華南理工大學(xué)電力學(xué)院，廣東 廣州 510640； 2 廣東省能源高效清潔利用重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510640）

引 言

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

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

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

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

本文擬在兩銅壁面之間構(gòu)建納米尺寸的冰立方模型，采用銅板原子靜態(tài)和脈沖荷電的方式來(lái)改變?cè)械臏?zhǔn)液體層厚度，通過(guò)對(duì)水分子施加垂直于壁面且大小一定的拉力，使得冰塊恰好脫離壁面，以測(cè)定其黏附強(qiáng)度。本研究擬量化同一溫度不同壁面荷電量時(shí)的準(zhǔn)液體層厚度及其對(duì)應(yīng)的冰黏附強(qiáng)度，以期提出減小冰黏附強(qiáng)度的有效方法。

1 模型和模擬

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

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

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

式中，E為兩粒子之間的勢(shì)能；ε為勢(shì)阱深度，描述了兩個(gè)粒子之間的相互作用強(qiáng)度；σ為兩個(gè)分子處于平衡位置時(shí)的距離；C是能量轉(zhuǎn)換常數(shù)；qi和qj是i原子和j原子上的電荷；r是兩個(gè)分子之間的距離；rc為對(duì)應(yīng)的截?cái)嗑嚯x。

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

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

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

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

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

1.2 平衡模擬和冰水判定

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

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

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

0.1 ps，每1 ps進(jìn)行一次數(shù)據(jù)采集。

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

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

1.3 脫附模擬和黏附強(qiáng)度的測(cè)量

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

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

式中，F(xiàn)是基底對(duì)冰塊的吸附力或拉力附加給冰塊的拉力；A是冰-基底接觸面積。

本文采用給所有水分子施加恒定拉力的方法進(jìn)行黏附強(qiáng)度測(cè)量。在每種工況下選擇10 個(gè)不同時(shí)刻的構(gòu)型，測(cè)定每個(gè)構(gòu)型中冰塊在2.5 ns 內(nèi)脫附的最小拉力，再將10個(gè)構(gòu)型中的最小脫離應(yīng)力作為該工況下的黏附強(qiáng)度。

2 結(jié)果與分析

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

2.1 冰塊生長(zhǎng)和準(zhǔn)液體層形成

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

圖6 類冰水分子數(shù)目隨時(shí)間的變化及對(duì)應(yīng)時(shí)刻的模型快照（黃銅色為銅原子，紅色為類液水分子，藍(lán)色為類冰水分子，白色為氫原子）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 “l(fā)iquid” water molecules, blue are “solid” water molecules, and white are hydrogen atoms)

2.2 冰的黏附強(qiáng)度

在大塊的水中結(jié)冰需要漫長(zhǎng)的形核過(guò)程[47]，這提供了充分的時(shí)間來(lái)測(cè)量過(guò)冷水的黏附強(qiáng)度。因此，本文對(duì)4 種工況均進(jìn)行了10 次初始構(gòu)型不同（體系平衡后取連續(xù)的10.0 ns，每隔1.0 ns 選取一個(gè)構(gòu)型）的相互獨(dú)立的脫附模擬。

以T=255 K，Qperiod=±0.1123 e/nm2時(shí)為例，取t=33.5 ns 的構(gòu)型為初始構(gòu)型。脫離應(yīng)力的測(cè)量需要不斷更改給定的拉力大小，根據(jù)冰塊在上一測(cè)試?yán)ο率欠衩摳綄?duì)拉力進(jìn)行調(diào)整，每次調(diào)整的拉力大小間隔約使外加法向應(yīng)力增加或減少1.0 MPa。圖7記錄了在測(cè)定該構(gòu)型脫離應(yīng)力過(guò)程中，其中6種不同拉力下冰塊的質(zhì)心（center of mass，COM）縱坐標(biāo)隨時(shí)間的變化。下銅板的荷電情況也標(biāo)在圖7 中。可以觀察到，脫附模擬過(guò)程中電荷產(chǎn)生的電場(chǎng)方向發(fā)生了兩次變化。冰塊質(zhì)心坐標(biāo)的迅速升高標(biāo)志著冰塊從銅板上剝離。在外加法向應(yīng)力σpull=135.4和140.3 MPa 時(shí)，冰塊質(zhì)心保持在同一水平，沒(méi)有發(fā)生脫附的現(xiàn)象。將外加法向應(yīng)力增加到141.2 MPa時(shí)，冰塊在拉力作用2.4 ns 后發(fā)生脫附。然而將外加法向應(yīng)力增加到143.1 MPa 時(shí)，冰塊卻意外地沒(méi)有脫附。當(dāng)進(jìn)一步將外加法向應(yīng)力增加到144.1 MPa 和147.0 MPa 時(shí)，冰塊脫附的時(shí)間點(diǎn)有所提前，且在較大外力作用時(shí)具有更大的提前幅度，這是因?yàn)橥饧臃ㄏ驊?yīng)力增大，脫附概率增加，可以在更短的時(shí)間內(nèi)實(shí)現(xiàn)脫附。綜上，在T=255 K、Qperiod=±0.1123 e/nm2的工況下，t=33.5 ns 時(shí)刻的脫離應(yīng)力為141.2 MPa。

圖7 T=255 K，Qperiod=±0.1123 e/nm2，t=33.5 ns時(shí)，不同外加法向應(yīng)力下的脫附模擬中冰塊質(zhì)心的位置變化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 時(shí)間段內(nèi)的10 個(gè)構(gòu)型（每1.0 ns 選取一個(gè)構(gòu)型）的脫離應(yīng)力測(cè)定結(jié)果如圖8 所示。在樣本范圍內(nèi)，t=35.5 ns時(shí)的脫離應(yīng)力最小，為130.6 MPa。綜上，T=255 K、Qperiod=±0.1123 e/nm2的工況下，冰的黏附強(qiáng)度為130.6 MPa。

圖8 T=255 K，Qperiod=±0.1123 e/nm2，t=33.5～42.5 ns時(shí)，各個(gè)冰塊構(gòu)型在不同外加法向應(yīng)力下的脫附情況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次獨(dú)立脫附模擬得到的脫離應(yīng)力。以這10 個(gè)脫離應(yīng)力最小值作為這種工況下的黏附強(qiáng)度，如圖9(b)所示。

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

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

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

圖10 體系平衡時(shí)的冰塊靜態(tài)結(jié)構(gòu)（虛線為區(qū)分冰和水的閾值，虛線之上為冰，虛線之下為水）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 結(jié) 論

本文使用全原子模型的分子動(dòng)力學(xué)模擬研究了光滑銅表面上黏附冰的脫附行為，采用序參量進(jìn)行水分子的冰水判定，以冰塊質(zhì)心位置和脫離應(yīng)力大小對(duì)冰塊的脫附過(guò)程進(jìn)行了定量表征，量化了255 K 時(shí)不同壁面條件下冰塊的黏附強(qiáng)度，得到如下結(jié)論。

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

符 號(hào) 說(shuō) 明

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

C——能量轉(zhuǎn)換常數(shù)

E——兩粒子之間的勢(shì)能，kcal/mol

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

lOM——四點(diǎn)TIP4P/ICE 水模型中無(wú)質(zhì)量的電荷點(diǎn)與O原子的距離，?

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

Ql——鍵序參數(shù)

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

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

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

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

r——兩個(gè)分子之間的距離，?

rc——對(duì)應(yīng)的截?cái)嗑嚯x，?

r0——四點(diǎn)TIP4P/ICE水模型中O—H鍵鍵長(zhǎng)，?

T——熱力學(xué)溫度，K

t——平衡時(shí)間，ns

ε——?jiǎng)葳迳疃?，kcal/mol

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

σ——兩個(gè)分子處于平衡位置時(shí)的距離，?

下角標(biāo)

period——脈沖荷電

static——靜態(tài)荷電