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

鄭晟 胡開鑫

(寧波大學(xué)機(jī)械工程與力學(xué)學(xué)院力學(xué)系 寧波 315211)

0 引言

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

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

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

1 控制方程

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

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

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

邊界條件如下。

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

Reynolds數(shù)Re、Marangoni數(shù)Ma和Prandtl數(shù)Pr分別定義為

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

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

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

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

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

其中，σ為復(fù)頻率，α和β分別表示在x和y軸上的空間波數(shù)。總波數(shù)和波的方向角分別用

表示。σ可以通過Chebyshev 配點(diǎn)法求解。下標(biāo)0 表示基本流，以下方程中無下標(biāo)量表示擾動。

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

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

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

2 數(shù)值結(jié)果

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

圖2 線性流Pr=0.002,k=4,?=90°不同Reynolds 數(shù)下瞬態(tài)增長函數(shù)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及對應(yīng)的時間tmax明顯增大。對于反饋函數(shù)，可以類似定義最大擾動放大?max=max?(ω)。圖3 表明在小Pr下Gmax和?max近似與Re2成正比。圖4 顯示了大Pr下不同Pr的反饋函數(shù)?(ω)。可以看到，回流大Pr下的亞臨界流動中（回流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下擾動能量主要來自表面熱毛細(xì)力做功。比較瞬態(tài)增長和擾動放大的量級發(fā)現(xiàn)|?max|2?Gmax，表明熱毛細(xì)液層對外加噪聲比初始狀態(tài)更為敏感。

圖3 線性流Pr=0.002,k=4,?=90°最大瞬態(tài)增長Gmax 和最大擾動放大?max隨Reynolds 數(shù)的變化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 數(shù)下擾動放大 ?隨實(shí)頻率 ω 的變化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 數(shù)和Reynolds 數(shù)的變化Fig.5 Maximum response function ?max versus the Prandtl numbers and Reynolds numbers at k=2.3,?=0°

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

圖6 線性流Pr=0.002、Ma=3 在不同實(shí)頻率下α ?β平面內(nèi) 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 在不同實(shí)頻率下α ?β平面內(nèi) 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下輸入和輸出的擾動流場及溫度場。可以看到擾動輸出速度和溫度的量級遠(yuǎn)大于輸入的量級。在小Pr下輸入和輸出的流場中均有流向條紋（沿流速度大于或小于平均流速的狹窄區(qū)域）。此放大機(jī)制為“舉起”機(jī)制[16]，即流向條紋是由流向渦放大產(chǎn)生的。輸入和輸出流場在豎直方向上的渦數(shù)量分別為1 和2，而在平面管道流動中，渦數(shù)量都是1[16]。輸入溫度在液層底部為0，而輸出溫度在全場均有分布。在大Pr下，輸入速度和溫度只分布于靠近液面處，輸出溫度分布在液層內(nèi)部，而在平面管道流動中，輸入速度分布于全流場[16]。大Pr下Re的量級僅為O(1)，此時溫度場十分關(guān)鍵，放大主要是熱毛細(xì)效應(yīng)。Smith[22]發(fā)現(xiàn)大Pr下熱毛細(xì)液層的不穩(wěn)定性主要是由垂直方向熱對流引起的。當(dāng)Pr增大時，熱對流比熱傳遞更重要，因此?max隨Pr的增大而增大（見圖3 和圖4）。

圖8 線性流Pr=0.002,Ma=3,k=2,?=90° 擾動放大對應(yīng)的擾動場。(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° 擾動放大對應(yīng)的擾動場。(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 結(jié)論

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

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

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

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

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