水平微電極表面氣泡生長過程Marangoni對流數值模擬

詹水清，袁睿，黃雨捷，張偉，王軍鋒

(江蘇大學 能源與動力工程學院，江蘇 鎮江，212013)

發展氫能是我國向可再生零碳能源結構轉型，實現“碳達峰、碳中和”目標的必由之路。可再生能源電解水制氫技術具有技術成熟、原料豐富、清潔低碳和安全高效等優點，備受研究者廣泛關注[1]。電解水制氫電極表面氫氣泡生長及脫離等行為引發氣泡局部微擾動，促進局部傳質。同時，氫氣泡在電極表面連續黏附，減小實際電化學活性反應面積，阻礙電極界面電子/離子傳輸，提高電極過電位和歐姆壓降等，增大電解能耗、降低效率[2]。因此，深入研究析氫電極表面氣泡動力學規律對于突破電解水制氫能耗高、能量轉化效率偏低等瓶頸問題至關重要。

為探明電解水制氫電極表面氣泡動力學行為，諸多學者主要采用高速攝影可視化、電化學測量和理論分析等方法研究氣泡生長、黏附和脫離等行為。當溶解氫分子濃度達到過飽和濃度條件時，電極表面成核點處產生氣泡，隨后氣泡繼續在電極表面生長和黏附。學術界主要采用理論分析結合電化學實驗的方法研究電極表面氣泡生長過程中的氣液傳質行為[3-6]。當電極尺寸明顯大于氣泡直徑時(常指宏觀電極)，目前普遍認為氣泡生長行為受氣液界面擴散主導控制，并提出了相關氣液界面擴散傳質速率模型[7-8]。而對于電極尺寸非常小的微電極表面，目前僅從理論上推測氣泡生長行為受微液層直接注入控制。

宏觀平面電極表面氣泡成核點分布的隨機性導致實驗測量較為困難。近年來相關研究聚焦于微電極表面單個氫氣泡行為[9-12]，獲得了較多微電極表面氫氣泡動態生長過程中的氣泡生長直徑、生長時間及相關電化學實驗數據，而針對氣泡底部微液層結構及演變等數據非常匱乏。CHEN等[13]通過理論分析，建立了基于相界面潤濕動力學的氣泡底部微液層模型；BASHKATOV等[14]通過微電極實驗，捕捉到氣泡底部存在不穩定的微液層結構，認為微液層主要是由多個微小子氣泡組成，微液層內的微小氣泡群之間的聚并行為對氣泡生長及脫離過程起到關鍵作用。此外，電解液歐姆熱阻產生焦耳熱，導致氣液界面呈現電解液溫度梯度分布，引發Marangoni對流效應[15]，從而影響氣泡生長和脫離等行為。微電極表面析氫氣泡尺度非常小，目前僅靠實驗研究難以深入揭示微液層結構和Marangoni對流效應的影響機制，而采用數值模擬手段能夠細致描述和求解微尺度氫氣泡周圍溫度場和速度場等信息。以往針對電極表面氣泡生長行為數值模擬研究主要采用固定氣泡直徑模型[16-18]。實際的氣泡動態生長速度相對于Marangoni對流速度來說非常小，基于此探索若干不同氣泡生長時刻對應不同直徑的單個氣泡周圍多物理場信息及影響因素。CHEN等[19-20]采用固定氣泡直徑模型研究電極表面氣泡生長過程中Marangoni對流行為，但沒有深入探究微電極體系考慮微液層結構對傳熱過程和溫度場的影響，也沒有深入研究單個連續氣泡生長周期內Marangoni對流結構演變及影響規律。

本文作者基于前期獲得的水平微電極表面單個氫氣泡生長行為實驗數據，采用固定氣泡直徑模型對實驗電流密度條件下連續多個氣泡生長時刻的不同直徑氣泡周圍的Marangoni對流效應進行數值模擬，開發基于微液層模型的電場分布和傳熱過程計算方法，揭示氣液界面溫度場及Marangoni對流結構的演變規律。

1 計算模型及求解方法

1.1 物理模型及網格劃分

本文采用前期設計的水平微電極電解水制氫實驗系統為參考物理模型[12]，將水平鉑微電極(直徑為100 μm)鑲嵌在高硬度的環氧樹脂玻璃平面中心，實驗采用1.5 mol/L 的H2SO4溶液作為電解液，在常溫25 ℃環境下進行電解，在水平微電極表面周期性地產生單個氫氣泡。為了盡量減少計算機耗時，并保證計算模擬的準確性，建立了以底部圓形微電極中心的法線為中心軸的三維圓柱幾何模型，其直徑為5 mm，高度為5 mm，實際電解實驗系統尺寸遠大于本文設計的計算區域尺寸。電極反應誘導產生的焦耳熱主要在氣泡表面附近電解液區域內傳輸，導致氣液界面周圍、特別是氣泡底部區域產生溫度梯度，從而僅主要在氣液界面引發Marangoni對流效應。因此，本文僅對氣泡周圍局部區域的流動和傳熱過程進行模擬是可行的，設計的三維幾何模型尺寸是合理的。實驗中獲得的氣泡動態生長速度為0.01～1.00 mm/s，明顯小于Marangoni 對流速度(后文詳細展示)，因此，本文采用固定氣泡直徑模型進行計算也是合理可靠的。

圖1所示為局部放大顯示的微電極表面單個氫氣泡生長示意圖。本文模擬設計的氣泡半徑Rb和氣泡接觸半徑Rc來源于恒電流密度為5 A/cm2條件下的氣泡幾何尺寸實驗，定義氣泡弧長從氣泡接觸點開始到氣泡頂點結束(圖1 中深藍色弧線)，氣泡接觸角為θ，氣液界面圓周角度為φ(φ≥θ)。計算區域的氣泡底部微液層為圓柱結構，微液層底部直徑與微電極直徑相同。為簡化求解，參考文 獻[14]，本文設計氣泡底部平均微液層厚度為 5 μm。在實際電解過程中，氣液界面和氣泡底部微液層區域的物理場參數變化較大，為保證計算結果的準確性，需要對這2個區域進行網格加密處理。圖2所示為三維圓柱模型中的氣泡垂直中心截面和陰極底部面的網格劃分和加密處理情況，從圖2可見：氣液界面和氣泡底部微液層區域的網格非常小，氣液界面最小網格尺寸為1.2 μm，氣泡底部微液層最小網格尺寸為0.5 μm。

圖1 微電極表面單個氫氣泡生長示意圖Fig.1 Diagram of single hydrogen bubble growth on microelectrode surface

圖2 2種關鍵截面網格劃分Fig.2 Local mesh at two key sections

1.2 數學模型及控制方程

電解水制氫體系內氣泡周圍的電解液流速較小，流動形態為明顯的層流。電解液為不可壓縮流體，描述單個氫氣泡生長過程中的流體流動方程包括連續性方程和動量方程，采用Boussinesq假設來考慮電解液溫度變化引起的浮力和自然對流過程，

式中：u為電解液速度；ρ為電解液密度；P為流體壓力；μ為電解液黏度；g為重力加速度；t為時間；β為流體熱膨脹系數；T為電解液溫度；T0為初始環境溫度(常溫25 ℃)，相關模擬結果主要用電解液溫度與初始環境溫度的差值(即溫度差)進行描述和討論。

電解液歐姆熱阻產生焦耳熱傳輸過程，引發氣液界面產生熱Marangoni對流效應，描述單個氣泡周圍電解液溫度變化的能量方程為

式中：λ為電解液導熱系數；cp為電解液定壓比熱容；Sheat為焦耳熱量源項，主要與局部電流密度和電導率有關。

電解過程的電勢分布通過拉普拉斯方程進行求解，電流密度分布由電勢梯度和局部電導率進行求解，

式中：φ為電解電勢；j為電解電流密度；κm為局部混合電導率，與局部區域的純電解液電導率κ和氣體含量αg有關，其計算公式為

氣泡周圍區域和微液層區域的電解液電導率不同，其中氣泡周圍為純電解液，而微液層內幾乎為大量的微/納米子氣泡組成，僅有少量的電解液，因此微液層區域的混合相電導率非常小，不考慮電解液電導率的溫度補償效應的影響。根據計算的局部電流密度和混合相電導率，計算得到焦耳熱源項

1.3 邊界條件及求解方法

對于考慮Boussinesq假設的電解液自然對流過程，為了提高計算收斂穩定性，采用非穩態迭代方法進行計算，當監測的氣液界面若干位置處的溫度和速度保持不變，即認為計算過程達到收斂。定義微電極表面為實際實驗的電流密度，定義陽極表面為壓力出口。將氣液界面定義為熱Marangoni對流邊界條件，表面張力隨溫度變化的系數為-1.5×10-4N/(m·K)[15]，其他壁面定義為無滑移壁面邊界條件。采用用戶自定義標量函數UDS求解槽內電勢分布，混合相電導率和焦耳熱源項均采用了用戶自定義函數UDF 進行編程計算。計算過程中的連續性方程和動量方程均采用二階迎風格式，能量方程和UDS 方程采用Quick 格式，壓力-速度耦合采用Phase-Coupled SIMPLE 算法，所有計算方程的松弛因子均采用默認值，時間步長為1 ms。模擬所需的電解液物性參數如表1 所示，微液層內氣體的體積分數約為94%[11]，計算得到微液層混合相電導率為2.28 S/m。為了系統研究氣泡生長周期內氣泡周圍Marangoni對流效應的演變規律，對不同氣泡生長時刻對應的不同氣泡尺寸模型進行數值模擬，結果如表2所示。

自古以來，數學家們都致力于揭示現象背后的本質，牛頓作為人類歷史上最偉大的數學家和物理學家之一，他利用數學解釋物理現象，并且創立了微積分。數學模型可以解釋事物背后的隱蔽模式，今天數學家和應用者們從實際中提煉出數學問題，再尋找合適的數學算法來解題，從而建立模型，這些模型可以應用到復雜、多變的自然現象、人類行為、社會系統等問題，微積分讓我們能夠更加深刻認識實數的性質，認識世界的本質。微積分的誕生極大地推動了力學、光學、熱學等各個領域的科技發展，促進了現代學科專業的發展[1]。

表1 電解液物性參數Table 1 Main physical properties of electrolyte

表2 實驗中不同氣泡生長時刻的氣泡直徑[12]Table 2 Bubble diameter at different bubble growth time from electrolysis experiment

2 網格無關性計算

圖3 所示為氣泡直徑為789 μm 時不同網格數量下的氣液界面Marangoni對流速度變化曲線。在不同網格數量下，Marangoni對流速度隨氣泡弧長增大均呈現先急劇增大、后緩慢減小的趨勢。當氣泡弧長大于0.1 mm時，Marangoni對流速度幾乎相同，因此，主要對比氣泡接觸點附近的氣泡弧長區域的Marangoni對流速度分布曲線。當網格數量增大時，Marangoni對流速度逐漸增大；當網格數量大于313 萬時，Marangoni 對流速度曲線幾乎相同。為提高計算效率和保證計算精度，后續相關數值模擬均采用313萬的網格數量進行網格劃分。

圖3 網格無關性計算Fig.3 Mesh independent calculation

3 溫度場和Marangoni對流結構

微電極表面氣泡生長過程中，電極附近的電解液產生焦耳熱，特別是微液層內產生的熱量比較大，導致氣泡周圍形成電解液溫度分布特性。如前所述，氫氣泡周圍的物理場參數呈現以圓形微電極中心的法線為中心軸的對稱結構，因此，圖4僅顯示單個氫氣泡垂直中心右邊截面的電解液溫度場分布與Marangoni 對流分布結果。由圖4(a)可知：氣泡底部微液層內溫度較大，氣液界面周圍電解液溫度較小。這是由于氣泡電絕緣和電流線集中扭曲效應所導致的，微液層內電流密度較大和電導率較小，導致電解液焦耳熱主要產生于微液層內。這些熱量隨著對流擴散逐漸影響到更多的氣泡周圍電解液區域，從氣泡接觸點到氣泡頂點的氣泡弧長方向，電解液溫度逐漸降低，溫度最大值在氣泡接觸點位置，導致氣液界面附近區域形成明顯的溫度梯度分布，距離氣液界面較遠區域內的溫度近似保持為環境溫度25 ℃。

圖4 單個氫氣泡周圍溫度場與Marangoni對流分布Fig.4 Temperature field and Marangoni convection distribution around single hydrogen bubble

氣液界面溫度梯度分布導致形成表面張力梯度分布，其中氣泡接觸點附近表面張力較小，沿著氣泡弧長方向表面張力逐漸增大，從而促使電解液從低表面張力向高表面張力方向運動，形成明顯的不穩定的Marangoni 對流旋渦結構，如圖4(b)所示。氣液界面附近區域的電解液沿著氣泡弧長方向流動，Marangoni對流速度逐漸減小，最大速度出現在氣泡接觸點附近的右上側位置。由于微液層內不考慮溫度梯度引起的Marangoni對流結果，因此，微液層內的電解液流動速度幾乎為零，近似認為微液層邊緣區域為固體邊界條件，導致Marangoni 對流到達接觸點位置受到一定的阻礙，使最大速度出現在氣泡接觸點附近的右上側位置。Marangoni對流旋渦的外邊緣區域速度較大，旋渦內部區域的速度較小，具體Marangoni對流旋渦結構及演變規律與整個氣泡生長周期不同時刻的氣泡直徑有關。

4 Marangoni對流演變規律

為了深入探明氫氣泡單個生長周期內氣液界面溫度梯度誘導產生的Marangoni對流結構的演變規律，本文對表2所示的不同氣泡直徑和接觸直徑下的Marangoni 對流效應進行了數值模擬。圖5 所示為微電極表面電極中心點到氣泡接觸點區域內的電解液溫度差變化規律，每條溫度差曲線最右側端點的橫坐標為氣泡接觸半徑Rc。由圖5 可見：隨著氣泡不斷生長，氣泡直徑持續增大，氣泡接觸直徑先增大、再近似保持不變、最后再減小且近似保持不變。從電極中心點到氣泡接觸點位置區域內，不同氣泡生長時刻的電解液溫度差均呈現先緩慢減小、然后快速減小的趨勢。微電極邊緣附近區域和接觸點區域內的溫度差減小趨勢非常明顯，這與圖4(a)一致。這是由于靠近氣泡底部的氣液界面的Marangoni對流速度較大，電解液焦耳熱能夠快速傳輸到更大范圍的電解液本體區域內。當氣泡直徑逐漸增大時，電解液溫度差逐漸增大，但是溫度差增大的趨勢逐漸減小。氣泡生長引起氣泡直徑增大，進而導致氣泡電絕緣和電流線集中扭曲效應也更加明顯，促使微液層內電流密度進一步增大，產生更多的焦耳熱量。當氣泡直徑從189 μm增大到467 μm時，氣泡直徑和接觸直徑增大明顯，上述解釋可進一步闡述溫度差增大幅度較大的原因。而當氣泡直徑大于578 μm時，氣泡生長過程中的接觸直徑變化趨勢非常小，且氣泡直徑增大的幅度也較小，因此，總體上電解液溫度差的增大幅度較小。

圖5 不同氣泡直徑(不同氣泡生長時刻)下微液層內溫度差變化Fig.5 Variation of temperature difference of microlayer with different bubble diameters (different growth time)

圖6所示為不同氣泡直徑下氣液界面電解液溫度差和Marangoni 對流速度的變化規律。由圖6(a)可知：在任意氣泡直徑下，氣液界面電解液溫度差均沿著氣泡弧長方向逐漸減小。氣泡直徑越大，氣液界面電解液溫差也較大，這主要與電極中心到氣泡接觸點區域內溫度分布及演變規律有密切關聯。在相同氣泡弧長下，大直徑氣泡的氣液界面溫度差比較大，意味著氣泡生長后期的大直徑氣泡周圍的溫度梯度較大。不同氣泡直徑下，明顯溫度梯度均主要體現在氣泡接觸點附近大約 0.3 mm 距離內的氣泡弧長界面上，距離氣泡接觸點稍遠的氣液界面上幾乎不存在溫度梯度。因此，當表面張力隨溫度變化的系數保持不變時，大直徑氣泡界面的表面張力梯度也較大，從而導致氣液界面Marangoni對流速度近似隨著氣泡直徑增大而增大，但這種變化規律對于不同直徑的氣泡來說，體現在不同的氣泡弧長區域范圍內。由圖6(b)可知：當氣泡直徑較小時，如氣泡直徑小于 733 μm，上述變化規律主要體現在除氣泡接觸點附近極小區域之外的其他絕大部分氣泡弧長界面內。當氣泡直徑僅為189 μm，氣液界面Marangoni對流速度反而較大，這可能與氣泡生長初期極小氣泡的接觸直徑及氣泡周圍的Marangoni 對流旋渦結構有關。而當氣泡直徑較大時，如氣泡直徑大于844 μm，上述變化規律主要體現在接觸點附近區域的部分有限氣泡弧長界面內。不同氣泡生長時刻下的氣液界面Marangoni 對流速度大于 10 mm/s，數量級明顯大于實驗獲得的氣泡生長速度0.01～1.00 mm/s，這再次說明本文建立的固定氣泡直徑模型是可靠合理的。

圖6 不同氣泡直徑(不同氣泡生長時刻)下溫差和Marangoni對流速度Fig.6 Temperature difference and Marangoni convection velocity with different bubble diameters (different growth time)

圖7 所示為不同氣泡直徑下氣液界面Marangoni 對流旋渦結構的變化規律。相對于圖4(b)的局部顯示結果，圖7給出了不同直徑氣泡周圍的更加直觀的Marangoni 對流旋渦形態及演變過程。當氣泡直徑較小時，旋渦結構中心大概在氣泡的右上端邊緣區域。隨著氣泡生長變大，氣泡周圍的Marangoni對流旋渦形態發生演變，特別是旋渦結構中心逐漸上移。對于任意氣泡直徑下，相對于氣泡弧長的上端頂點所在的高度而言，旋渦結構中心位置的高度較小。隨著氣泡直徑增大，氣液界面最大Marangoni對流速度所在的氣泡弧長圓周角度φ逐漸減小，但都大于對應的氣泡接觸角θ，如圖8 所示。最大電解液溫差出現在氣泡接觸角θ所在的氣泡弧長位置上，進一步說明沿著氣泡弧長方向，最大Marangoni對流速度的位置稍微偏離于氣泡接觸點位置。

圖7 不同氣泡直徑(不同氣泡生長時刻)下Marangoni對流結構演變Fig.7 Variation of Marangoni convection structure with different bubble diameters (different growth time)

圖8 不同氣泡生長時刻下最大電解液溫差和Marangoni對流速度的角度位置Fig.8 Angle positions of maximum temperature difference and Marangoni convection velocity with different growth time

5 結論

2) 氣液界面溫度梯度分布導致形成表面張力梯度分布，促使電解液從氣泡底部沿著氣泡弧長方向流動，形成不穩定的Marangoni 對流旋渦結構。沿著氣泡弧長方向的Marangoni對流速度逐漸減小，最大速度出現在氣泡接觸點附近的右上側界面位置。

3) 從電極中心點到氣泡接觸點位置區域內，電解液溫度差先緩慢減小、后快速減小。當氣泡直徑逐漸增大時，微液層內電解液溫度差逐漸增大，但是溫度差增大的趨勢逐漸減小。

4) 在相同氣泡弧長下，氣泡直徑越大，氣液界面溫度差越大。氣液界面Marangoni對流速度近似隨著氣泡直徑增大而增大，但這種變化規律對于不同直徑的氣泡來說，體現在不同的氣泡弧長區域范圍內。具體Marangoni對流旋渦結構及演變規律與整個氣泡生長周期不同時刻的氣泡尺寸有關。