熱毛細液層非模態(tài)穩(wěn)定性分析*

鄭晟 胡開鑫

(寧波大學機械工程與力學學院力學系 寧波 315211)

0 引言

熱毛細對流是由液體表面溫差引起表面張力梯度驅動的流動，是空間等微重力條件下流體的主要運動形式，廣泛存在于微流控[1，2]、熔焊[3]、薄膜鍍層[4]等工業(yè)應用中。鑒于其在晶體生長中的重要性[5]，近幾十年來圍繞熱毛細對流已開展了大量研究[6-7]，包括理論分析、數值模擬以及地面和空間實驗。Davis[8]和Schatz等[9]分別在1987年和2001年對相關內容進行了綜述。

熱毛細流動的穩(wěn)定性直接影響晶體生長的質量，因此一直受到廣泛關注。Smith等[10]采用線性穩(wěn)定性方法研究了無限大平板上的熱毛細液層，發(fā)現其失穩(wěn)機制取決于Prandtl 數。Chan等[11]應用上述液層模型研究了重力效應對熱毛細穩(wěn)定性的影響，其臨界參數與空間實驗結果[12]相符。2016年Kang等[13]在中國 實踐十號返回式科學實驗衛(wèi)星中對環(huán)形液池的熱毛細對流進行了空間實驗，結果顯示流體從穩(wěn)態(tài)到振蕩態(tài)過程中存在定向傳播的熱流體波.

以往熱毛細對流失穩(wěn)的理論研究多采用模態(tài)分析方法，即假定小擾動隨時間呈指數變化。雖然該方法可預測擾動的長時間變化，但是對于許多情況，例如管道流動，模態(tài)分析得到的臨界Reynolds 數與實驗觀測不符。流動轉捩對初始狀態(tài)和外加擾動十分敏感。因此有研究發(fā)展了非模態(tài)理論，研究流動失穩(wěn)的短期效應[14-17]。對于熱毛細對流，在一些情況下模態(tài)分析得到的臨界Reynolds數Rec量級達到Rec≈O(104)或O(105)[18]，遠大于一般情況下的轉捩Reynolds數，因此可能存在亞臨界條件下的流動失穩(wěn)。此外非模態(tài)理論已應用于一些熱毛細力驅動的流動，研究結果與實驗相符[19]。本文采用非模態(tài)方法研究熱毛細液層的流動穩(wěn)定性對初始擾動和外加激勵的敏感性，并 討論Reynolds 數和Prandtl 數的影響。

1 控制方程

采用圖1 所示熱毛細液層模型，即厚度為d的無限長液體薄層因液面溫度梯度而產生對流。x、y、z分別表示流向、展向和法向。流動為平行剪切流（包括回流和線性流）。

圖1 熱毛細液層模型Fig.1 Schematic of thermocapillary liquid layers

流動的無量綱控制方程組分別為連續(xù)性方程、動量方程和能量方程，即

邊界條件如下。

其中：u、p、T、t分別為速度、壓強、溫度和時間；τ為應力張量；T∞為無窮遠處溫度；Q?為外加熱通量。在剛性平面處（z=0），要求速度無滑移且豎直方向無熱通量；在自由面處（z=1），前兩式表示熱毛細力提供流動剪切力，第三式表示法向速度為0，第四式表示豎直方向熱通量與T ?T∞呈線性關系。

Reynolds數Re、Marangoni數Ma和Prandtl數Pr分別定義為

其中，ρ為流體密度，U為流體特征速度，d為液層厚度，μ為動力粘度，χ為熱擴散系數，h為傳熱系數，κ為熱導率。為了簡略起見，在這里只討論Bi=0的情況。流體為牛頓流體，則本構方程為τ=S，應變率張量S=?u+(?u)T。

任宏等[6]認為，該區(qū)域的礦床形成主要經歷了海相熱鹵水噴流-沉積交代階段、變質疊加改造階段和近地表淺部氧化淋濾次生富集階段三個階段，銅鈷元素的初始富集主要發(fā)生在同生沉積交代成礦階段。熱水同生沉積巖主要由海底熱水發(fā)生同生沉積作用而形成，有少量陸源碎屑物質的混入。巖石發(fā)育紋層狀、厚層塊狀等構造，硅質巖中金屬硫化物呈浸染狀、層紋狀等產出[24]。

考慮如下兩種形式的基本流[8]。線性流：

回流：前者速度呈線性分布，而后者法向質量流量為0。

以下分析流動穩(wěn)定性。首先采用模態(tài)分析方法，在基本流場中疊加正則小擾動，有

其中，σ為復頻率，α和β分別表示在x和y軸上的空間波數?？偛〝岛筒ǖ姆较蚪欠謩e用

表示。σ可以通過Chebyshev 配點法求解。下標0 表示基本流，以下方程中無下標量表示擾動。

采用非模態(tài)分析方法研究初始狀態(tài)和外加擾動對流動穩(wěn)定性的影響。設擾動演化方程具有如下形式：

其中，ψ=(u,v,w,T)為擾動量，L為演化算子（可以由模態(tài)分析得到），Θ為外加激勵。對于無外加激勵的情況，可用瞬態(tài)增長函數G(t)反映流場對初始擾動的放大，即t時刻擾動能量與初始時刻擾動能量之比；對于外加激勵為頻率ω的簡諧函數情況，可用反饋函數?(ω)反映流場對外加信號的放大，即輸出擾動能量與輸入擾動能量之比。計算公式分別為[20]

范數‖·‖2表示擾動能量E，可定義為[21]

2 數值結果

圖2 顯示了不同Re數下擾動的瞬態(tài)增長函數G(t)?？梢钥吹骄€性流小Pr下的亞臨界流動中（線性流Pr=0.002 的臨界Marangoni數Mac≈3.41）存在較大的瞬態(tài)增長，擾動經過一段短暫增長之后開始衰減。因此可以定義最大瞬態(tài)增長

圖2 線性流Pr=0.002,k=4,?=90°不同Reynolds 數下瞬態(tài)增長函數G(t)隨時間t 的變化Fig.2 Variation of G(t) with time of the linear flow at Pr=0.002, k=4,?=90° with various Reynolds numbers

隨著Re增大，Gmax及對應的時間tmax明顯增大。對于反饋函數，可以類似定義最大擾動放大?max=max?(ω)。圖3 表明在小Pr下Gmax和?max近似與Re2成正比。圖4 顯示了大Pr下不同Pr的反饋函數?(ω)?？梢钥吹剑亓鞔驪r下的亞臨界流動中（回流Pr=150 的臨界Mac≈383.7）存在較大的擾動放大。圖5顯示了回流大Pr下?max隨Pr和Re的變化，可以看到?max分別隨Re5和Pr5呈線性增長，漸近線分別為?max≈3.28Re5+12.05和?max≈1.92×10?9Pr5+15.45。能量分析表明大Pr下擾動能量主要來自表面熱毛細力做功。比較瞬態(tài)增長和擾動放大的量級發(fā)現|?max|2?Gmax，表明熱毛細液層對外加噪聲比初始狀態(tài)更為敏感。

圖3 線性流Pr=0.002,k=4,?=90°最大瞬態(tài)增長Gmax 和最大擾動放大?max隨Reynolds 數的變化Fig.3 Maximum transient growth function Gmax and the maximum response function ?max versus Reynolds numbers of the linear flow at Pr=0.002,k=4,?=90°

圖4 回流Re=2,k=2.3,?=0°不同Prandtl 數下擾動放大 ?隨實頻率 ω 的變化Fig.4 Response function ? versus the real frequency ω of the return flow at Re=2,k=2.3,?=0°with various Prandtl numbers

圖5 回流k=2.3,?=0°最大擾動放大?max隨Prandtl 數和Reynolds 數的變化Fig.5 Maximum response function ?max versus the Prandtl numbers and Reynolds numbers at k=2.3,?=0°

圖6 與圖7 分別顯示了線性流和回流對不同頻率和波數外加激勵的放大情況。定義最優(yōu)擾動放大可以發(fā)現，相比大Pr，小Pr下的流動存在更大的擾動放大，并且小Pr數下的流動在靠近90°時取到?opt，而大Pr下在0°時取到該值。在平面管道流動中，?opt均存在于90°附近[16]。隨著頻率ω的增大，線性流小Pr下的最優(yōu)擾動從β軸（β ≈4）變?yōu)棣?≈0.45,β ≈3；大Pr回流的最優(yōu)擾動從α ≈2.3減小到α ≈1.5。此結果說明隨著外加激勵頻率ω的增大，最優(yōu)擾動波數逐漸減小。

圖6 線性流Pr=0.002、Ma=3 在不同實頻率下α ?β平面內 lg?的等值線Fig.6 Level lines of the logarithm of the response lg?in the α ?β plane for the linear flow at Pr=0.002,Ma=3 with various real frequency

圖7 回流Pr=150，Ma=300 在不同實頻率下α ?β平面內 lg?的等值線Fig.7 Level lines of the logarithm of the response lg?in the α ?β plane for the return flow at Pr=150,Ma=300 with various real frequency

圖8 和圖9 顯示了不同Pr下輸入和輸出的擾動流場及溫度場?？梢钥吹綌_動輸出速度和溫度的量級遠大于輸入的量級。在小Pr下輸入和輸出的流場中均有流向條紋（沿流速度大于或小于平均流速的狹窄區(qū)域）。此放大機制為“舉起”機制[16]，即流向條紋是由流向渦放大產生的。輸入和輸出流場在豎直方向上的渦數量分別為1 和2，而在平面管道流動中，渦數量都是1[16]。輸入溫度在液層底部為0，而輸出溫度在全場均有分布。在大Pr下，輸入速度和溫度只分布于靠近液面處，輸出溫度分布在液層內部，而在平面管道流動中，輸入速度分布于全流場[16]。大Pr下Re的量級僅為O(1)，此時溫度場十分關鍵，放大主要是熱毛細效應。Smith[22]發(fā)現大Pr下熱毛細液層的不穩(wěn)定性主要是由垂直方向熱對流引起的。當Pr增大時，熱對流比熱傳遞更重要，因此?max隨Pr的增大而增大（見圖3 和圖4）。

圖8 線性流Pr=0.002,Ma=3,k=2,?=90° 擾動放大對應的擾動場。(a) 輸入速度場，(b) 輸出速度場，(c) 輸入溫度場，(d) 輸出溫度場Fig.8 Perturbation fields corresponding to the response for the linear flow at Pr=0.002, Ma=3, k=2,?=90°.(a) Input velocity field,(b) output velocity field,(c) input temperature field and(d) output temperature field

圖9 回流Pr=150,Ma=300,k=2.3,?=0° 擾動放大對應的擾動場。(a) 輸入速度場，(b) 輸出速度場，(c) 輸入溫度場，(d) 輸出溫度場Fig.9 Perturbation fields corresponding to the response for the return flow at Pr=150, Ma=300,k=2.3,?=0°.(a) Input velocity field,(b) output velocity field,(c) input temperature field and(d) output temperature field

3 結論

利用非模態(tài)穩(wěn)定性方法研究了熱毛細液層對初始擾動和外加激勵的放大，得到以下結論。

（1）線性流和回流小Pr下的亞臨界流動對初始擾動和外加激勵均十分敏感，最大瞬態(tài)增長Gmax和最大擾動放大?max與Re2近似成正比。

（2）回流大Pr下的亞臨界流動中存在對外加激勵的顯著放大，?max分別隨Re5和Pr5呈線性增長。

（3）相比于大Pr，小Pr的流動存在更大的擾動放大。隨著外加激勵頻率ω的增大，最優(yōu)擾動波數逐漸減小。

（4）流場和溫度場表明，擾動輸出速度和溫度的量級遠大于輸入的量級。在小Pr下，輸入和輸出的流場中均有流向條紋，但豎直方向的渦數量及溫度分布不同；在大Pr下輸入的速度和溫度只分布于靠近液面處，而輸出溫度分布在液層內部。兩者放大機制不同。