附加不同豎直顫動下表面內能變化導致的熱毛細對流界面響應行為研究

楊 碩，周曉峰，鄧喬聲，王國峰

（沈陽工程學院 遼寧省潔凈燃燒發電與供熱技術重點實驗室，沈陽 110136）

與地面環境相比，許多地球上難以獲得或無法獲得的深空環境包含了大量有待探索的物理、化學和生命科學問題，并由此形成了一批前沿學科——微重力科學。這一新興學科和前沿課題的發展又孕育了新型空間高技術產業的誕生。人類在“深空”極端環境方面的涉足為制備高質量、大尺寸半導體元件提供了可能。如何掌握深空極端環境材料生長過程背后的微重力流體物理傳熱與傳質規律，成為多相流物理學基礎研究中亟待解決的重要課題。2020年前后，我國將建成載人空間站核心艙與實驗艙，其中已規劃的微重力流體物理實驗平臺為空間微重力極端環境下毛細對流和界面現象的研究創造了有力條件。雖然深空環境下，重力水平被極大削弱（表面張力梯度引起的對流運動突顯），但由于空間輻射、飛船姿態調整等因素產生的重力跳動和殘余應力仍然存在，周期性顫動的體積力使得胞元流產生相應的運動和熱力學參數變化，并影響毛細對流系統界面形狀，流動平衡和運動規律。已有空間實驗表明：重力跳動和殘余應力效應與振蕩熱毛細對流頻率相關。高頻重力顫動作用下，熱毛細對流系統的響應時間遠大于重力顫動的固有周期，在一定控制參數下，其對熱毛細對流自發振蕩行為將產生積極作用。因此，研究重力顫動對熱毛細對流系統的影響及其響應機制，有效利用這種作用減緩非平衡態熱毛細對流不穩定性具有重要意義。

2005年，Kanashima等［1］研究了重力顫動對液橋中振蕩熱毛細對流引起的動態表面形變的影響。在壓電驅動器驅動的振動實驗臺上模擬重力顫動施加的振幅和頻率，并利用微成像位移計測量重力顫動下的界面變形。結果表明：重力顫動引起的振蕩熱毛細對流自由面變形按頻率可劃分為諧波和非諧波。Ichikawa等［2］基于“massspring-damper”模型，考察水平橫向顫動下液橋毛細對流及界面響應行為，數值結果與“SpaceLab D-1”空間實驗結果吻合。Liu等［3］采用線性不穩定性分析方法，研究了在高頻顫動下2層流體系統中Rayleigh-Marangoni-Benard對流的不穩定性。結果表明：高頻熱動可以改變2層流體系統中Marangoni-Benard對流的不穩定性和Rayleigh-Marangoni-Benard振蕩間隙。垂直高頻振動可以延遲對流不穩定性，抑制對流向下流動。Liang等［4］用數值方法研究了絕熱條件下單向以及多向顫動對三維液橋表面波的影響，外加橫向顫動的頻率達到反響頻率時，可以完全控制液橋表面水平方向的振動。Kawaji等［5］利用level set方法數值捕捉附加重力顫動條件下三維液橋毛細對流的界面行為，并進行液橋共振頻率的預測，數值結果表明：共振頻率是隨著液橋高度、直徑以及密度的增加而減小，隨表面張力的增加而增大，這一結果與Ichikawa等預測結論相一致。2014年，Lyubimova等［6］基于廣義Boussinesq近似，研究高頻豎直顫動對軸對稱半浮區熱毛細流動的影響及其三維擾動的線性穩定性。采用有限差分法在不同Weber數（We）下計算振動幅值和普朗特數對系統的影響。結果表明：豎直顫動可用來控制并降低流動強度，使流動穩定。2016年，Lyubimova等［7］研究外加高頻微小豎直顫動對具有恒定垂直熱通量的2層不相溶流體穩定性的影響。結果表明：長波不穩定性不受小、中等強度顫動的影響。當顫動瑞利數足夠大時，單調有限波長擾動會加劇破壞2層流體系統的穩定性。Vjatkin等［8］在實驗和理論上研究了轉動頻率和振動頻率接近時，橫向顫動對旋轉圓柱空腔內非等溫流體對流的影響。研究表明：橫向顫動是控制旋轉系統熱對流的有效手段。當顫動頻率接近旋轉頻率時，對流受到激發并加劇熱傳遞。在一定頻率下，非等溫流體柱內發生二維方位角共振，并導致中心軸處的溫度明顯降低。熱傳遞明顯依賴于共振液體的振蕩及頻率。Kovskaya等［9］研究高頻水平簡諧顫動對帶有自由表面的液層熱毛細對流的影響。結果表明：顫動的縱向分量對熱毛細對流不穩定性沒有影響。如果穩定持續的顫動同時包含橫向分量，則豎直顫動具有抑制熱毛細對流不穩定性的很好效果。Kovskaya等［10］研究豎直顫動對具有自由表面不可壓縮均勻流體液層Marangoni對流的影響。得到3個主要失穩類型的臨界失穩參數方程，相應地同步參數共振域、次諧波，以及高頻漸近形式的頻率值被確定。研究表明：在豎直顫動下長波振蕩Marangoni對流的Ma數只依賴于Pr數和Biot數（Bi）。

將空氣域與液橋區形成的兩相自由界面動力學形變及周圍空氣的影響考慮進來，模擬殘余重力顫動及微重力環境，數值研究高Pr數液橋熱毛細對流胞元流受迫運動規律，對于大密度比及高黏度比為流體介質的液橋熱毛細對流，采用質量完全守恒水平集方法追蹤氣液界面響應行為。

1 物理模型

1.1 控制方程

本文研究對象的幾何模型如圖1所示，液橋半徑和高度分別為R和H，置于2個軸對稱圓盤的中間，上、下圓盤的溫度分別為高溫Tt及低溫Tb，溫差為ΔT＝Tt－Tb。液橋周圍被空氣包圍，氣相區域的外徑為2R。FGX為附加重力顫動力，微重力條件下豎直方向重力場加速度g＝0 m／s2（液橋系統內浮力效應導致的浮力對流被忽略）。為考察自由界面在豎直顫動下的反響行為，參考系原點設置在液橋左側空氣域最左側邊界處。

圖1 帶有豎直顫動的半浮區液橋坐標系統

由以下無量綱化的質量守恒方程，Navier-Stokes方程和能量守恒控制方程描述本文的研究對象。

式中：u＝（u，v）為流體流速；ρ＝ρ（x，t）為流體密度；μ＝μ（x，t）為流體動力黏度；D是黏性應力張量；κ是自由界面的曲率；d為計算質點到界面的法向距離；δ為狄拉克函數（δ（x）＝0，x≠0），；n為界面處的單位法向矢量；t為時間項；P為壓力項。

表面張力系數被看作是溫度的線性函數，表示為σ＝σc－σT（T－Tb）；其中，σc是環境溫度下（T0＝25℃）表面張力的參考值；σT是表面張力的溫度系數；σT＝?σ／?T，T為界面流溫度。FGX是由附加重力顫動引起的外力是y軸方向上的外部加速度，A是顫動幅度，ω角頻率（ω＝2πf）；f為顫動頻率，為特征長度，定義特征長度＝D，其中D為液橋的初始直徑。U∞是特征速度則微重力下，有為無量綱時間，表示為。此外，本課題組認為液橋介質為不可壓縮牛頓流體，對連續性方程、動量方程和能量方程采用Boussinesq近似。

在表面內能變化數值研究方面，將界面形狀特性和溫度特性差異完整引入界面內能變化，

式中：σ為液橋自由界面表面張力；Θliq為非等溫液橋自由表面溫度；Θair為環境側溫度；Q0為由非等溫液橋側向環境側的熱通量。（4）中等號右側第1項表征界面溫度引起的表面內能變化項；第2項表征自由面形變導致的表面內能變化項。其他無量綱參數定義如下：和分別對應環境空氣區域內的無量綱密度比和動力黏度比；ρl和μl分別是液橋區域流體介質的密度和動力黏度；ρg和μg分別是環境空氣域內流體介質的密度和動力黏度。液橋內部流體介質的無量綱密度比（ρl／ρl）和無量綱黏度比（μl／μl）均等于1。雷諾數定義為；韋伯數定義為We＝；普朗特數定義為Pr＝μl／ρla；熱毛細數定義為Ca＝μlU∞／σ，熱Marangoni數定義為Ma＝＝RePr，絕對溫度的無量綱形式為Θ＝（T－Tb）／ΔT；a是熱擴散率。

1.2 水平集方法及重整化處理

水平集方法最初是由Osher等［11］提出，用于數值預測兩相流體介質之間的運動界面。Level set方法在整個計算域內隱式跟蹤運動界面。在水平集方法中，?（x，t）代表兩相界面函數，液橋界面空氣域?（x，t）＞0，液橋界面處為?（x，t）＝0，液橋界面內部?（x，t）＜0。通過求解Level set函數?（x，t）可以預測界面的運動位置。

基于上面對液橋界面內外介質密度和黏度的處理，結合Level set函數界面位置，流場中的密度函數ρ（?）、黏度函數μ（?）和Heaviside函數可以分別表示為：

其中，為避免求解過程中界面附近密度和黏度階躍性變化帶來的計算不穩定性，在計算中選取α＝2.0Δx，α表示液－氣界面厚度，Δx是x方向上的網格距離。

為避免Level set函數求解過程中自身帶來的質量失衡，Sussman等［12］提出采用方程（9）（10）來實現Level set函數重新初始化。

式中：?0（x）和?（x）有著相同的水平集數。然而即使采用重新初始化程序，總質量在時間尺度上仍然不能完全滿足完全的質量守恒。這主要是由于水平集函數數值離散過程中自身的缺陷，為克服計算失真問題，質量守恒過程必須通過求解如下面積補償（二維模型）方程來實現。

式中：A（t）是水平集函數?（t）在時間t時所對應的液橋面積（二維模型），A0是初始條件下液橋的初始面積。面積約束函數F（c）的定義如下：

F（c）隨著h和n的變化而變化，當h＝0和n＝0時，數值計算過程的收斂速度會更快，收斂性的判別條件如下：

1.3 邊界條件和初始條件

計算區域所有的邊界均假設為絕熱壁面（除液橋的熱盤、冷盤和自由界面之外），熱盤和冷盤的初始溫度分別保持為Tt和Tb（隨后保持恒定的升溫速率），溫差為ΔT，滿足邊界條件如下：

本文所研究的兩相系統中，需要考慮初始靜止狀態和液橋周圍空氣的初始速度。

無滑移邊界條件同樣適用于液橋上、下盤。

采用質量守恒水平集法分析液橋的自由界面運動，在交錯網格上求解Navier-Stokes方程和能量守恒方程，對平流項采用QUICK迎風插值法進行離散，對其他項采用中心有限差分法進行離散，對時間積分項采用2階Adams-Bashforth方法離散。

對于表面張力的求解計算中采用CSF（continuum surface force）方法。假設在液橋自由界面上表面張力連續分布，由Level Set函數可以將表面張力項表示為：

其中，界面曲率κ（φ）由κ（φ）＝－▽·（n）和▽·n＝▽·確定。

2 對比驗證

為了驗證本文模型的正確性，將液橋內部等溫線計算結果與Shevtsova等［13］的數值結果進行比較。在零重力場條件下，以10 cSt硅油液橋作為驗證對象（R＝3 mm，H＝4 mm，ΔT＝40 K，Re＝166，Pr＝105）。如圖2所示，計算結果與Shevtsova等的結果（b）較為吻合，但本文的計算方法考慮了液橋自由界面的動態變形過程。

圖2 數值方法的對比驗證（Pr＝105，Re＝166，Γ＝R／h＝4／3）

3 豎直顫動對胞元流渦心位置及自由界面行為的影響

3.1 豎直顫動下胞元流渦心在軸向上的位置變化

為了研究微重力條件下豎直顫動對液橋表面以及內部熱毛細對流胞元流產生的作用，設定附加垂直顫動力為20 mg（FGX），并分別選取顫動頻率f＝0 Hz，f＝5 Hz，f＝20 Hz進行模擬。在數值計算中，液橋介質選用5cSt硅油：Re＝320，Ma＝22 037，We＝3.5，Pr＝66.93，Ca＝5，上、下圓盤的初始溫差為ΔT＝25℃，氣相、液相的密度比和黏度比分別為ρg／ρl＝1.35×10－3和μg／μl＝5.528×10－3，其他物性參數見表1。

表1 5cSt硅油物性參數

圖3給出了數值結果分析中液橋幾何位置關系的主要信息。本文建立的計算域原點并不在液橋的中軸線上（如圖1所示），但研究對象為穩定熱毛細對流，熱毛細對流的胞元流型態呈現軸對稱特征，下文中為了簡化表達只給出了一半液橋的數值結果。

圖3 計算模型的幾何位置關系

以液橋區域流函數絕對值的最大值作為捕捉渦心位置的判據，流函數絕對值的變化意味著渦心位置的變化。在沒有施加豎直顫動的情況下（f＝0 Hz），渦心位置隨熱毛細對流發展，渦心存在自發性的擺動，如圖4所示（f＝0 Hz）。在不同頻率外加顫動的作用下（f＝5 Hz，f＝20 Hz），液橋熱毛細對流胞元流渦心位置在軸向上的變化如圖4所示，施加的豎直顫動加速度水平相同，均為20 mg。圖4（a）（b）分別表示初始階段胞元流渦心位置的變化和顫動后期胞元流渦心位置的變化，y軸表示渦心距冷端（下盤）的縱向無量綱距離。由于熱毛細對流首先在熱角發起，所以初始時刻渦心高度位于熱角附近，且總體擺動幅度較大（最大擺動幅度（Amax＝0.25）。外加顫動對流動在軸向方向上產生影響，使熱毛細對流在軸向上受到抑制。在經過初始階段后，熱角區的活躍性逐漸減弱，胞元流渦心逐步向冷端移動，并處于穩定的小振幅狀態（最大擺動幅度Amax＝0.10），顫動狀態下的平衡位置如圖4（b），且圍繞著某一平衡位置為上限（y＝0.55）進行周期性波動。但無論是在“初始階段”還是“穩定階段”，隨著豎直顫動頻率的增加，熱毛細對流胞元流渦心位置軸向擺動幅度逐漸減小。

圖4 豎直顫動頻率對胞元流渦心軸向位置的影響曲線（Pr＝66.93，Ma＝22 037，FGX＝20 mg）

3.2 豎直顫動下胞元流渦心在徑向上的位置變化

如圖5所示，熱毛細對流胞元流的渦心在豎直顫動力的作用下，不僅影響渦心在軸向上的位置分布，還會影響其在徑向上的位置分布。圖5（a）和圖5（b）分別為胞元流渦心徑向位置在“初始階段”和“穩定階段”的分布。在初始階段，當沒有附加豎直顫動力時（f＝0 Hz）渦心徑向位置穩定的位于R＝0.6處（R為液橋無量綱半徑）。隨著豎直顫動頻率的增加（f＝5 Hz，f＝20 Hz），渦心徑向位置圍繞著R＝0.6處進行周期性不規則擺動。豎直顫動頻率為f＝20 Hz時，渦心徑向位置擺動幅度加劇。在穩定階段，渦心徑向位置呈現周期性的擺動特點，并且擺動的下限維持在R＝0.55處，與未施加豎直顫動力相比，渦心徑向位置擺動的上限仍然維持在R＝0.60。從圖5可知，附加豎直顫動會迫使胞元流渦心徑向位置呈現穩定性的周期擺動，且擺動周期隨顫動頻率的增加而減小。

圖5 豎直顫動頻率對胞元流渦心徑向位置的影響曲線（Pr＝66.93，Ma＝22 037，FGX＝20 mg）

3.3 附加豎直顫動下液橋自由界面波動的響應行為

豎直顫動會引起液橋自由界面產生微弱的高頻振蕩，并與胞元流渦心軸向和徑向上的擺動共同構成豎直顫動下復雜熱毛細對流流動結構演化行為，在以往的研究中這種變形很難從實驗中直接定量觀察。在豎直正弦顫動頻率f＝0、5、20 Hz下，分別選取液橋表面上距離冷端（下盤）垂直高度1／40H、1／4H、1／2H、3／4H處作為監測點，定量考察自由界面橫向位置擺動（如圖6所示），圖中y軸表示液橋自由界面橫向位置偏離液橋半徑的變化量。

圖6 豎直顫動頻率對自由界面不同高度處表面橫向位置振幅的影響曲線（Pr＝66.93，Re＝320，Ma＝22 037，FGX＝20 mg）

從圖6（a）中可以看出：未加入豎直顫動情況下（f＝0 Hz），液橋自由界面橫向位置存在微弱的大周期擺動，但是隨著時間的推移，液橋表面橫向位置振幅逐漸穩定，且A1／2H＝1.5×10－4＞A3／4H＝1.2×10－4＞A1／4H＝1.0×10－4＞A1／40H＝0（中間高位置處自由界面橫向位置偏離原始位置R最大），這與熱毛細對流胞元流的發展、壯大以及渦心游動有關；施加豎直顫動頻率為5、20 Hz情況下（如圖6（b）（c）所示），自由界面橫向位置呈現出周期性明顯的擺動。在豎直顫動頻率f＝5 Hz情況下，液橋表面高度1／40H處自由界面橫向位置沒有發生顫動。液橋界面高度1／4H處自由面橫向位置的振幅ΔAm＝0.5，周期為TAm＝25。液橋表面高度1／2H處自由界面橫向位置振幅最小。液橋表面高度3／4H處自由界面橫向位置的振幅均大于其他位置（ΔAm＝0.7，TAm＝25）；與f＝5 Hz相比，在豎直顫動頻率f＝20 Hz情況下，液橋自由界面橫向位置脈動性擺動明顯，自由界面橫向位置振幅的大小關系如下，ΔA3／4H＞ΔA1／4H＞ΔA1／2H＞ΔA1／40H。綜上所述，隨著豎直顫動頻率的增加，自由界面橫向位置的擺動周期減小，振幅基本沒有改變。

3.4 附加豎直顫動下液橋自由界面軸向和徑向速度的響應行為

液橋自由界面表面流的軸向速度分量如圖7所示。軸向速度分量v為正，表示速度方向與y軸正向一致；軸向速度分量v為負，表示與y軸負向一致，液橋表面流由熱角（上盤）流向冷角（下盤）。

圖7 直顫動頻率對自由界面表面流軸向速度分量的影響曲線（Pr＝66.93，Re＝320，Ma＝22 037，FGX＝20 mg）

在沒有施加豎直顫動情況下（f＝0 Hz，如圖8（a）所示），初始階段t＝8～25 s的軸向速度v值較大，隨著熱毛細對流的發展和渦心的游動，t＝100～200 s時刻的軸向速度逐漸衰減；相對于初始階段，t＝100～200 s階段液橋內的熱毛細對流已接近穩定流動。在施加豎直顫動頻率f＝5 Hz情況下（如圖7（b）所示），液橋表面流軸向速度v明顯衰減，最大衰減幅度達Δv＝3.5×10－3，這說明豎直顫動力的施加抑制了液橋自由界面表面流軸向速度向冷端（下盤）的發展。

液橋自由界面表面流的徑向速度分量如圖8所示，徑向速度分量u為正，表示速度方向與x軸正向一致，說明表面流由液橋中心流向界面；徑向速度分量u為負，表示速度方向與x軸負向一致，說明表面流由界面流向液橋中心。

圖8 豎直顫動頻率對自由界面表面流徑向速度分量的影響曲線（Pr＝66.93，Re＝320，Ma＝22 037，FGX＝20 mg）

如圖8（a）所示，液橋自由界面表面流徑向速度u在距冷端無量綱距離y＝0.55以上為正，在液橋自由界面高度y＝0.0～0.55范圍內，表面流徑向速度u為負，即液橋表面流徑向速度在靠近熱角處為正（上盤），靠近冷角處為負（下盤），且徑向速度峰值分別位于y＝0.9和y＝0.2。同表面流軸向速度表現相近，初始階段（t＝8～25 s）的速度波動范圍較大，而在穩定階段（t＝100～200 s）速度波動范圍明顯變小。在附加豎直顫動力頻率為f＝5 Hz情況下（如圖8（b）所示），隨著頻率增加徑向速度減小，且穩定階段（t＝100～200 s）徑向速度明顯小于未施加豎直顫動下的徑向速度，最大衰減幅度達Δu＝1.0×10－3。

3.5 附加豎直顫動下液橋自由界面位形的響應行為

圖9為附加豎直顫動20 mg、20 Hz條件下，液橋自由界面位形的變化情況，選取液橋右半部分作為分析對象（如圖3所示），左半部分為液橋域，右半部分為空氣域，從圖9中可以看出：在無量綱時間t＝860～864 s，液橋自由界面高度y＝0.45～1.0部分的表面逐漸向液側收縮；在時間段t＝866～870 s，液橋自由界面高度y＝0.45～1.0部分的表面逐漸向氣側擴張。與此同時，液橋冷端處表面（y＝0～0.45）向液側收縮。因而，液橋自由界面在附加豎直顫動力作用下，液橋近熱端和近冷端的表面呈現交替狀的擴張和收縮，且上、下2部分的轉折點距冷端的無量綱距離為y＝0.45，不在液橋的中間高位置y＝0.5。

圖9 豎直顫動對液橋自由界面位形的影響曲線（Pr＝66.93，Re＝320，Ma＝22 037，FGX＝20 mg，f＝20 Hz）

4 結論

1）豎直顫動對熱毛細對流胞元流的作用可以由胞元流渦心軸向和徑向位置的變化體現。無論是在“初始階段”還是“穩定階段”，隨著豎直顫動頻率的增加（f＝0 Hz→f＝5 Hz→f＝20 Hz），渦心存在的自發性擺動受到抑制，熱毛細對流胞元流渦心位置軸向擺動幅度逐漸減小，最大擺動幅度由ΔAmax＝0.25衰減至ΔAmax＝0.10。受豎直顫動的影響，胞元流渦心徑向位置由原來無規則擺動呈現出穩定性的周期擺動，且擺動周期逐漸減小。

2）隨著豎直顫動頻率的增加，液橋自由界面橫向位置的擺動周期減小，振幅基本沒有改變，脈動特性明顯。在豎直顫動頻率20 Hz情況下，自由界面橫向位置振幅的大小呈如下排布，ΔA3／4H＞ΔA1／4H＞ΔA1／2H＞ΔA1／40H。

3）豎直顫動對液橋自由界面表面流速度的影響：豎直顫動力能夠抑制液橋自由界面表面流軸向速度向冷端（下盤）發展。隨著豎直顫動頻率的增加，穩定階段（t＝100～200 s）徑向速度明顯減小，且最大衰減幅度達Δu＝1.0×10－3。

4）在豎直顫動力的作用下，液橋近熱端（y＝0.45～1.0）和近冷端（y＝0～0.45）的表面呈現交替性的擴張和收縮。