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

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

0 引 言

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

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

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

1.1 物理模型

模擬所選二維模型如圖1所示,左邊界OF為恒溫面,其高度為0.24 m,記為H。為避免水平邊界對豎直邊界層產(chǎn)生影響[10],分別向豎直恒溫面上、下方向延長0.2H?？紤]到黏性邊界層的厚度影響[11],模型的寬度設(shè)置為0.5H。擾動區(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ù);φ0為初始流體溫度;φw為垂直恒溫面溫度。

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

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

(5)

Pr=γ/κ

(6)

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

1.3 擾動設(shè)置

擾動區(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取決于恒溫面和初始流體間的溫差,即模擬通過改變恒溫面的溫度,達到改變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)格獨立性檢測

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

表3 獨立性檢測工況

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

圖2 獨立性檢測

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

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

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

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

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

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

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

1.6 計算求解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖6 R5工況下測點1～測點4處溫度時間序列

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

圖7 R5工況下測點1～測點4處溫度時間序列功率譜

為確定不同Ra的自然對流熱邊界層的特征頻率,對工況R1～工況R7,即對隨機式擾動下,不同Ra狀態(tài)下的自然對流進行模擬,導(dǎo)出不同工況下測點4的溫度時間序列,并對其進行快速傅里葉變換,得到不同Ra下自然對流邊界層的特征頻率,進而得到特征頻率與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的值進行擬合,得到特征頻率fc與Ra之間的關(guān)系式,即fc=0.226Ra0.641。

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

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

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

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

3 結(jié) 論

通過對豎直恒溫面處的自然對流現(xiàn)象進行模擬,探究不同Ra狀態(tài)下的自然對流邊界層的不穩(wěn)定性。模擬首先驗證了豎直恒溫面處自然對流邊界層的不穩(wěn)定性,以及熱邊界層可分為上游低頻區(qū)、過渡區(qū)、下游高頻區(qū)3個區(qū)域。另外,模擬結(jié)果表明在無擾動、隨機式擾動和正弦式擾動狀態(tài)下,隨著Ra增大,豎直恒溫壁面處的自然對流熱邊界層整體厚度均減小;在隨機式擾動和正弦式擾動下,豎直恒溫面處的自然對流熱邊界層沿豎直方向存在振蕩,且相比隨機式擾動,正弦式擾動下的熱邊界層的震蕩幅度更大;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。