瑞利數(shù)對豎直恒溫面處自然對流邊界層不穩(wěn)定性的影響

王春雨,劉澤勤,劉聰聰,趙鵬鵬

(1.天津商業(yè)大學(xué) 機(jī)械工程學(xué)院,天津 300134; 2.天津商業(yè)大學(xué) 冷凍冷藏技術(shù)教育部工程研究中心,天津 300134; 3.天津商業(yè)大學(xué) 天津市制冷技術(shù)重點(diǎn)實(shí)驗(yàn)室,天津 300134; 4.揚(yáng)州大學(xué) 水利與能源動(dòng)力工程學(xué)院,江蘇 揚(yáng)州 225009)

0 引 言

邊界層不穩(wěn)定性特征是自然對流研究領(lǐng)域中的重要課題,認(rèn)識邊界層不穩(wěn)定性,有助于理解層流至湍流的轉(zhuǎn)化、強(qiáng)化對流換熱等問題。早在20世紀(jì)50年代,有關(guān)邊界層不穩(wěn)定性的研究已經(jīng)展開。文獻(xiàn)[1]最早通過馬赫-曾德爾干涉儀(Mach-Zehnder interferometer)觀察受到擾動(dòng)的層流邊界層,并發(fā)現(xiàn)熱邊界層能夠放大該擾動(dòng),并最終過渡至湍流流動(dòng),該現(xiàn)象即為邊界層的不穩(wěn)定性;文獻(xiàn)[2-3]研究了使用絲帶在豎直面附近產(chǎn)生震動(dòng),進(jìn)而觀察擾動(dòng)對邊界層的影響,實(shí)驗(yàn)結(jié)果表明,豎直邊界層僅能放大特定頻率內(nèi)的擾動(dòng)。當(dāng)邊界層同時(shí)受到低頻和高頻擾動(dòng)時(shí),邊界層的不穩(wěn)定性取決于高頻擾動(dòng)[4-5],而非低頻擾動(dòng)。模擬中常用隨機(jī)變化的熱流區(qū)域或溫度邊界使邊界層產(chǎn)生不穩(wěn)定性,文獻(xiàn)[6]探究在豎直恒溫面的上游位置引入隨機(jī)式變化的熱流下,觀測自然對流邊界層內(nèi)部溫度時(shí)間序列的功率譜。模擬結(jié)果表明在邊界層上游區(qū),低頻帶的功率值遠(yuǎn)高于高頻帶,即低頻擾動(dòng)占主導(dǎo)作用,隨著邊界層向下游發(fā)展,低頻帶不斷衰減,高頻帶對應(yīng)功率值升高。根據(jù)功率譜的頻帶分布,自然對流豎直邊界層可劃分為上游低頻區(qū)、過渡區(qū)(既有高頻帶也有低頻帶)以及下游高頻區(qū)[7-8]。利用邊界層不穩(wěn)定性特征能夠?qū)崿F(xiàn)共振強(qiáng)化對流換熱,文獻(xiàn)[9]研究在豎直邊界層上游引入正弦式變化的熱流區(qū)域,其變化頻率與邊界層的特征頻率相同,從而實(shí)現(xiàn)共振強(qiáng)化換熱。結(jié)果發(fā)現(xiàn),相比于無擾動(dòng)狀態(tài)下的對流換熱,在頻率為特征頻率的正弦擾動(dòng)下,對流換熱量提高了44%。

目前,關(guān)于豎直恒溫面處邊界層不穩(wěn)定性的研究相對較少,有關(guān)瑞利數(shù)Ra對邊界層不穩(wěn)定性影響的研究還需進(jìn)一步進(jìn)行。因此本文對不同Ra下,普朗特?cái)?shù)Pr=6.24的豎直恒溫面處的自然對流進(jìn)行模擬,探究Ra對自然對流邊界層不穩(wěn)定性的影響。

1 數(shù)學(xué)與物理模型

1.1 物理模型

模擬所選二維模型如圖1所示,左邊界OF為恒溫面,其高度為0.24 m,記為H。為避免水平邊界對豎直邊界層產(chǎn)生影響[10],分別向豎直恒溫面上、下方向延長0.2H?？紤]到黏性邊界層的厚度影響[11],模型的寬度設(shè)置為0.5H。擾動(dòng)區(qū)域設(shè)置在OE處,其范圍為0

圖1 模擬模型

表1 邊界條件

1.2 數(shù)學(xué)模型

假設(shè)流體為不可壓縮流體,且采用Boussinesq假設(shè),自然對流的二維無量綱控制方程組為:

(1)

(2)

(3)

(4)

其中:U=u/(γH-1);V=v/(γH-1);X=x/H;Y=y/H;t=τ/(H2γ-1);P=p/(ργ2H-2);T=(φ-φ0)/(φw-φ0);U、V、X、Y、t、P、T分別為u、v、x、y、τ、p、φ的無量綱形式;u、v分別為x、y方向的速度;p為壓力;τ為時(shí)間;ρ為密度;φ為溫度;γ為流體動(dòng)力黏性系數(shù);φ0為初始流體溫度;φw為垂直恒溫面溫度。

規(guī)定頻率為F,無量綱頻率為f,頻率按照f=FH2γ-1進(jìn)行無量綱化,本文所述頻率均為無量綱頻率。

自然對流的狀態(tài)取決于無量綱瑞利數(shù)Ra和普朗特?cái)?shù)Pr,計(jì)算公式為:

(5)

Pr=γ/κ

(6)

其中:g為重力加速度;β為熱膨脹系數(shù);κ為熱擴(kuò)散系數(shù);γ為流體動(dòng)力黏性系數(shù);Δφ=φw-φ0。

1.3 擾動(dòng)設(shè)置

擾動(dòng)區(qū)設(shè)置在范圍為0

1.4 模擬工況

模擬中,不同工況下流體的初始溫度均為298.15 K,該狀態(tài)下γ的值為0.892 92×10-6m2/s,β的值為0.000 25,κ的值為0.143×10-6m2/s,通過(6)式可以算得Pr為6.24。Ra的值與恒溫面和流體間的溫差、恒溫面的高度有關(guān),由于恒溫面的無量綱高度恒定為1(量綱長度為0.24 m),因此Ra取決于恒溫面和初始流體間的溫差,即模擬通過改變恒溫面的溫度,達(dá)到改變Ra的目的。模擬工況所選用的Ra分別為5.31×108、7.97×108、1.06×109、1.33×109、2.66×109、3.98×109、5.31×109,上述Ra對應(yīng)的壁面溫度分別為300、301、302、303、308、313、318 K,對應(yīng)的恒溫面與初始流體間的溫差分別為2、3、4、5、10、15、20 K。模擬工況見表2所列。

表2 模擬工況

1.5 網(wǎng)格獨(dú)立性檢測

分別對網(wǎng)格、時(shí)間步長和擾動(dòng)振幅進(jìn)行獨(dú)立性檢測,獨(dú)立性檢測工況見表3所列。選用不同網(wǎng)格數(shù)量、時(shí)間步長和擾動(dòng)振幅進(jìn)行模擬,記錄不同工況下相同測點(diǎn)(X=0.004 2,Y=0.830 0)處的溫度時(shí)間序列,再將測點(diǎn)的溫度時(shí)間序列導(dǎo)入MATLAB中進(jìn)行快速傅里葉變換,將時(shí)域內(nèi)的溫度信號處理至頻域的頻率信號,通過MATLAB導(dǎo)出溫度時(shí)間序列的功率譜,功率譜顯示了不同頻率所對應(yīng)的功率值,其中對應(yīng)功率值最大的頻率即為峰值頻率。對比不同工況下的峰值頻率和功率譜分布趨勢,完成獨(dú)立性檢測。

表3 獨(dú)立性檢測工況

首先對網(wǎng)格大小為336×131、448×164的2種網(wǎng)格(工況a、工況b)進(jìn)行模擬,模擬選用時(shí)間步長為Δt,其值為1.3×10-6,擾動(dòng)振幅為A,其值為1,獨(dú)立性檢測如圖2所示。

圖2 獨(dú)立性檢測

測點(diǎn)(X=0.004 2,Y=0.830 0)處溫度時(shí)間序列的功率譜如圖2a所示。功率譜的橫坐標(biāo)f為溫度變化的頻率,縱坐標(biāo)P(f)為每個(gè)頻率所對應(yīng)的功率值。

由圖2a可知,在2種不同網(wǎng)格大小下,功率圖分布趨勢大致相同,且對應(yīng)功率值最大的無量綱頻率相同,均為23 676(對應(yīng)的非無量綱頻率為0.33 Hz),規(guī)定該頻率為峰值頻率fc。為保證計(jì)算的精確度,且避免計(jì)算量過大,下面選用網(wǎng)格大小為336×131的模型。

選用網(wǎng)格大小為336×131的模型,分別使用Δt、Δt/2、2Δt的時(shí)間步長(工況a、工況c、工況d)進(jìn)行模擬,擾動(dòng)振幅為A,測點(diǎn)(X=0.004 2,Y=0.830 0)處溫度時(shí)間序列的功率譜如圖2b所示。

由圖2b可知,在3種不同的時(shí)間步長下,功率譜的分布規(guī)律相似,且峰值頻率fc相同,均為23 676。為保證計(jì)算的精確度,選用Δt的時(shí)間步長。

選用網(wǎng)格大小為336×131的模型,分別使用振幅為A、A/2、3A/2擾動(dòng)(工況a、工況e、工況f)進(jìn)行模擬,時(shí)間步長為Δt,得到測點(diǎn)7的溫度功率譜如圖2c所示。

由圖2c可知,在不同擾動(dòng)振幅下,功率譜的分布規(guī)律相似,且得到的峰值頻率fc相同,均為23 676。為保證計(jì)算準(zhǔn)確度,選擇振幅為A的擾動(dòng)。

綜上所述,本文選用網(wǎng)格大小為336×131、時(shí)間步長為1.3×10-6、振幅為1的模型進(jìn)行模擬。

1.6 計(jì)算求解

在ANSYS FLUENT定義運(yùn)行文件,模擬為非穩(wěn)態(tài)模擬,選用壓力-速度耦合以及層流模型。流體采用Boussinesq假設(shè),應(yīng)用有限體積法和SIMPLE算法求解控制方程,不穩(wěn)定項(xiàng)采用二階向后差分格式進(jìn)行積分,空間導(dǎo)數(shù)采用二階中心差分格式進(jìn)行離散。

根據(jù)獨(dú)立性檢測選用網(wǎng)格大小為336×131的模型,時(shí)間步長選用0.083 25 s,其對應(yīng)的無量綱時(shí)間步長為1.3×10-6。為記錄模擬開始后100 s內(nèi)的自然對流邊界層的變化情況,時(shí)間步長的步數(shù)為1 200,每一時(shí)間步長的最大迭代次數(shù)為20次,收斂方案中最大殘差值設(shè)為10-4,按上述內(nèi)容設(shè)置好后進(jìn)行求解工作。

2 模擬結(jié)果及分析

2.1 豎直恒溫面處自然對流邊界層的溫度特性

對表2中T1～T7工況,即不同瑞利數(shù)Ra下,豎直恒溫面處的自然對流進(jìn)行模擬。以T1、T5、T7工況為例, 3種工況下豎直恒溫面處的自然對流熱邊界層如圖3所示。

圖3 T1、T5、T7工況下的豎直熱邊界層

圖3中Y軸代表豎直恒溫面,即模型中的OF邊。由圖3可知,熱邊界層內(nèi)溫度分層均勻,邊界層流動(dòng)處于層流狀態(tài),且在3種不同工況下,T1工況下熱邊界層厚度最大,T7工況下熱邊界層厚度最小。該現(xiàn)象證明,隨著Ra增大,邊界層厚度減小,邊界層內(nèi)溫度梯度增大。

2.2 隨機(jī)擾動(dòng)下的自然對流邊界層不穩(wěn)定性

為探究自然對流邊界層的不穩(wěn)定性,確定熱邊界層的特征頻率,在擾動(dòng)區(qū)加入隨機(jī)式擾動(dòng),使擾動(dòng)區(qū)的溫度呈2A(Rand(0,1)-0.5)的數(shù)學(xué)形式變化。對工況R1～工況R7,即隨機(jī)式擾動(dòng)下,不同Ra下的自然對流進(jìn)行模擬。以R1、R5、R7工況為例, 3種工況下豎直恒溫面處的自然對流熱邊界層如圖4所示。由圖4a可知,在R1工況下,熱邊界層出現(xiàn)明顯的波動(dòng),非線性增強(qiáng),且隨著熱邊界層由上游發(fā)展至下游,其波動(dòng)程度增強(qiáng)。由圖4可知,R1工況下熱邊界層的厚度最大,R7工況下熱邊界層厚度最小,該現(xiàn)象證明隨著Ra增加,隨機(jī)式擾動(dòng)下的邊界層厚度逐漸減小,該現(xiàn)象與無擾動(dòng)狀態(tài)下的邊界層規(guī)律相同。

圖4 R1、R5、R7工況下的豎直熱邊界層

以工況T5 、工況R5為例,記錄并對比2個(gè)工況下測點(diǎn)(0.004 2,0.830 0)處的溫度時(shí)間序列,結(jié)果如圖5所示。將無擾動(dòng)與隨機(jī)式擾動(dòng)狀態(tài)下測點(diǎn)處的溫度時(shí)間序列劃分為增長階段、過渡階段和穩(wěn)定階段[10](圖5),在增長階段和過渡階段,測點(diǎn)的溫度變化趨勢相同,溫度值相差不大,在穩(wěn)定階段存在明顯差異。無擾動(dòng)狀態(tài)下,測點(diǎn)的溫度穩(wěn)定維持在定值處,隨機(jī)式擾動(dòng)下,測點(diǎn)的溫度則呈現(xiàn)出不規(guī)則波動(dòng)。該現(xiàn)象證明,邊界層上游的隨機(jī)式擾動(dòng)對邊界層在增長階段和過渡階段的影響不大,其影響主要體現(xiàn)于穩(wěn)定階段。

以R5工況為例,選取測點(diǎn)1(0.004 2,0.042 0)、測點(diǎn)2(0.004 2,0.170 0)、測點(diǎn)3(0.004 2,0.500 0)、測點(diǎn)4(0.004 2,0.830 0),并記錄在隨機(jī)式擾動(dòng)下測點(diǎn)處的溫度時(shí)間序列,結(jié)果如圖6所示。通過對比4個(gè)測點(diǎn)處的溫度時(shí)間序列,發(fā)現(xiàn)隨著Y坐標(biāo)的增大,即隨著邊界層高度的增加,測點(diǎn)在穩(wěn)定階段內(nèi)的溫度波動(dòng)增加,該現(xiàn)象證明豎直自然對流邊界層能夠放大上游位置的擾動(dòng)。其中,測點(diǎn)4在穩(wěn)定階段的溫度波動(dòng)最大,測點(diǎn)1在穩(wěn)定階段的溫度波動(dòng)最小。

為確定熱邊界層的特征頻率,將上述4個(gè)測點(diǎn)在穩(wěn)定階段的溫度時(shí)間序列(5×10-4

圖5 T5 、R5工況下測點(diǎn)(0.004 2,0.830 0)處溫度時(shí)間序列

圖6 R5工況下測點(diǎn)1～測點(diǎn)4處溫度時(shí)間序列

圖7中,頻率均為無量綱頻率,縱坐標(biāo)P(f)為功率值。觀察測點(diǎn)1(0.004 2,0.042 0)處溫度時(shí)間序列的功率譜,如圖7a所示。在熱邊界層內(nèi)測點(diǎn)1處,頻率范圍0～20 000對應(yīng)的功率值較大,且峰值頻率為10 762,該現(xiàn)象說明在邊界層上游,低頻帶占主導(dǎo)作用。觀察測點(diǎn)2(0.004 2,0.170 0)處溫度的功率譜,如圖7b所示。相對測點(diǎn)1而言,在熱邊界層內(nèi)測點(diǎn)2的位置,頻率范圍為0～10 000對應(yīng)功率值下降,頻率范圍為20 000～30 000對應(yīng)的功率值上升,且峰值頻率為23 676,該現(xiàn)象說明隨著邊界層高度的升高,占主導(dǎo)作用的頻率帶向高頻區(qū)移動(dòng)。觀察測點(diǎn)3(0.004 2,0.500 0)、測點(diǎn)4(0.004 2,0.830 0)處溫度的功率譜,如圖7c、圖7d所示,頻率范圍為15 000～30 000對應(yīng)的功率值較大,峰值頻率均為23 676,該現(xiàn)象說明在邊界層的中下游,高頻帶占主導(dǎo)作用。上述現(xiàn)象證明豎直邊界層的發(fā)展可分為上游低頻區(qū)、過渡區(qū)、下游高頻區(qū)3個(gè)階段。由于熱邊界層的不穩(wěn)定性取決于高頻擾動(dòng)[12],規(guī)定下游邊界層的高頻峰值頻率為特征頻率,即工況R5下的自然對流邊界層的特征頻率為23 676。

圖7 R5工況下測點(diǎn)1～測點(diǎn)4處溫度時(shí)間序列功率譜

為確定不同Ra的自然對流熱邊界層的特征頻率,對工況R1～工況R7,即對隨機(jī)式擾動(dòng)下,不同Ra狀態(tài)下的自然對流進(jìn)行模擬,導(dǎo)出不同工況下測點(diǎn)4的溫度時(shí)間序列,并對其進(jìn)行快速傅里葉變換,得到不同Ra下自然對流邊界層的特征頻率,進(jìn)而得到特征頻率與Ra之間的擬合關(guān)系式,結(jié)果如圖8所示。Ra為5.31×108、7.97×108、1.06×109、1.33×109、2.66×109、3.98×109、5.31×109狀態(tài)下的豎直自然對流邊界層所對應(yīng)的無量綱特征頻率分別為6 457、12 914、15 066、16 143、23 676、32 286、38 743,由此可知,在Ra為5.31×108～5.31×109的范圍內(nèi),隨著Ra增加,自然對流邊界層的特征頻率增加。

圖8 特征頻率fc與Ra的擬合關(guān)系式

應(yīng)用Origin軟件,對特征頻率fc與Ra的值進(jìn)行擬合,得到特征頻率fc與Ra之間的關(guān)系式,即fc=0.226Ra0.641。

2.3 正弦擾動(dòng)下自然對流邊界層不穩(wěn)定性

為了解頻率為特征頻率的正弦式擾動(dòng)對熱邊界層不穩(wěn)定性的影響,在擾動(dòng)區(qū)設(shè)置頻率為特征頻率的正弦式擾動(dòng),使擾動(dòng)區(qū)的溫度呈T=Asin(2πfct)的數(shù)學(xué)形式變化。對工況S1～工況S7進(jìn)行模擬,并導(dǎo)出各工況下自然對流熱邊界層圖。以S1、S5、S7工況為例,3種工況下豎直恒溫面處的自然對流熱邊界層如圖9所示。

由圖9a可知,工況S1下的熱邊界層產(chǎn)生劇烈波動(dòng),且隨著熱邊界層由上游發(fā)展至下游,波動(dòng)程度增強(qiáng)。對比圖4、圖9可知,相比于隨機(jī)式擾動(dòng),正弦式擾動(dòng)下熱邊界層的波動(dòng)程度更大。由圖9可知,隨著Ra增加,熱邊界層的整體厚度減小,整體溫度梯度增加,該現(xiàn)象與無擾動(dòng)狀態(tài)下的溫度梯度邊界層梯度規(guī)律相同。

圖9 S1、S5、S7工況下的豎直熱邊界層

3 結(jié) 論

通過對豎直恒溫面處的自然對流現(xiàn)象進(jìn)行模擬,探究不同Ra狀態(tài)下的自然對流邊界層的不穩(wěn)定性。模擬首先驗(yàn)證了豎直恒溫面處自然對流邊界層的不穩(wěn)定性,以及熱邊界層可分為上游低頻區(qū)、過渡區(qū)、下游高頻區(qū)3個(gè)區(qū)域。另外,模擬結(jié)果表明在無擾動(dòng)、隨機(jī)式擾動(dòng)和正弦式擾動(dòng)狀態(tài)下,隨著Ra增大,豎直恒溫壁面處的自然對流熱邊界層整體厚度均減小;在隨機(jī)式擾動(dòng)和正弦式擾動(dòng)下,豎直恒溫面處的自然對流熱邊界層沿豎直方向存在振蕩,且相比隨機(jī)式擾動(dòng),正弦式擾動(dòng)下的熱邊界層的震蕩幅度更大;Ra為5.31×108、7.97×108、1.06×109、1.33×109、2.66×109、3.98×109、5.31×109狀態(tài)下的豎直自然對流邊界層所對應(yīng)的無量綱特征頻率分別為6 457、12 914、15 066、16 143、23 676、32 286、38 743，Ra在5.31×108～5.31×109的范圍內(nèi),自然對流邊界層的特征頻率與Ra間的耦合公式為fc=0.226Ra0.641。