郭子漪 趙建福 李 凱 胡文瑞

(中國(guó)科學(xué)院力學(xué)研究所微重力重點(diǎn)實(shí)驗(yàn)室,北京 100190)

(中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

引言

熱毛細(xì)對(duì)流是由溫度梯度導(dǎo)致的界面張力梯度驅(qū)動(dòng)的對(duì)流,是空間自然對(duì)流的基本形式.作為一個(gè)流動(dòng)與傳熱耦合的復(fù)雜非線(xiàn)性過(guò)程,隨著驅(qū)動(dòng)力的增加,熱毛細(xì)對(duì)流會(huì)由穩(wěn)定的層流,經(jīng)歷各種各樣的分岔,發(fā)展為混沌流.熱毛細(xì)對(duì)流的轉(zhuǎn)捩過(guò)程受幾何形狀、熱邊界條件和物性參數(shù)等多種因素的影響[1-2].以往的實(shí)驗(yàn)[3]及數(shù)值模擬中研究了各種工況下的熱毛細(xì)對(duì)流,對(duì)臨界[4-6]及超臨界轉(zhuǎn)捩過(guò)程得到了豐富的結(jié)果,確認(rèn)了熱流體波的振蕩形式[7],發(fā)現(xiàn)了倍周期分岔,準(zhǔn)周期分岔,切分岔等轉(zhuǎn)捩途徑[8-10].然而,大多數(shù)的數(shù)值模擬研究采用直接數(shù)值模擬的方法,需要在不同參數(shù)下進(jìn)行大量數(shù)據(jù)計(jì)算,不但難以獲得完整的分岔圖,而且效率不高.近年來(lái),基于動(dòng)力系統(tǒng)中的分岔理論[11],發(fā)展出了一套數(shù)值分岔方法[12],通過(guò)求解分岔方程,計(jì)算分岔點(diǎn),得到更為完整的分支圖,此方法在流體力學(xué)的許多領(lǐng)域得到了成功的運(yùn)用[13-16].但是,由于流動(dòng)的控制方程是偏微分方程組,空間離散后化為維數(shù)極高的常微分方程組,這為準(zhǔn)確及高效求解帶來(lái)了困難.為了解決這一問(wèn)題,除了嘗試發(fā)展適用于大規(guī)模計(jì)算的數(shù)值方法,還可以考慮對(duì)原方程進(jìn)行降維,在低維模型上進(jìn)行數(shù)值分岔的計(jì)算.

常見(jiàn)的基于流場(chǎng)特征提取的降維方法主要有本征正交分解(POD)和動(dòng)力學(xué)模態(tài)分解方法[17](DMD).POD 方法的本質(zhì)是尋找高階系統(tǒng)的一組正交基,滿(mǎn)足原始流場(chǎng)在這組基上的方差最小,這個(gè)優(yōu)化問(wèn)題最終轉(zhuǎn)化為矩陣的特征值和特征向量的求解[18],其中特征值反映各個(gè)模態(tài)的能量,通常少數(shù)幾階模態(tài)就占絕大部分能量,將流場(chǎng)表示為上述主要模態(tài)的線(xiàn)性組合,從而實(shí)現(xiàn)降維.DMD 方法的本質(zhì)是將流動(dòng)的演化看作線(xiàn)性過(guò)程,假設(shè)前后時(shí)間步的流場(chǎng)相差一個(gè)線(xiàn)性變換.并通過(guò)快照矩陣估計(jì)變換矩陣的特征值和特征向量,構(gòu)造反映流動(dòng)特性的動(dòng)力學(xué)模態(tài).由于DMD 方法對(duì)流動(dòng)演化的過(guò)程進(jìn)行分析,得到的模態(tài)是時(shí)空耦合的,不需要額外建立方程就可以估計(jì)流動(dòng)隨時(shí)間的發(fā)展[19].

與DMD 方法相比,POD 方法無(wú)法計(jì)算流動(dòng)的演化過(guò)程,要想得到低維方程,需結(jié)合Galerkin 投影法,將原系統(tǒng)的解表示為特征模態(tài)的線(xiàn)性組合,利用POD 模態(tài)的正交性,通過(guò)Galerkin 投影將原偏微分方程化為以組合系數(shù)為未知量的常微分方程,上述構(gòu)造低維方程的方法一般被稱(chēng)為POD-Galerkin 降維法.POD-Galerkin 降維法的有效性在許多研究中得到了驗(yàn)證:文獻(xiàn)[20]構(gòu)建了二維矩形腔內(nèi)浮力對(duì)流的POD 降維模型,研究了不同模態(tài)數(shù)對(duì)模型穩(wěn)定性的影響;文獻(xiàn)[21]采取了混合不同參數(shù)下DNS 數(shù)據(jù)的方法,構(gòu)建了三維圓柱繞流的POD 低維模型;Cai 等[22]研究了二維方腔內(nèi)浮力對(duì)流的修正低維模型,并通過(guò)參數(shù)外推,驗(yàn)證了模型的魯棒性;馮俞楷[23]等以?xún)?yōu)化POD 解的梯度為目標(biāo),從溫度場(chǎng)的梯度中提取出最優(yōu)正交基,結(jié)合Galerkin 投影法,建立了二維、常物性、非穩(wěn)態(tài)導(dǎo)熱微分方程的PODGalerkin 降維模型;Jing 等[24]利用三維方腔內(nèi)浮力對(duì)流的低維模型進(jìn)行數(shù)值分岔的計(jì)算,并與直接數(shù)值模擬的結(jié)果做了對(duì)比.

值得注意的是,大多數(shù)工作都是將POD-Galerkin 降維方法用在比較典型的頂蓋驅(qū)動(dòng)流、圓柱繞流、浮力對(duì)流等模型中.降維方法在熱毛細(xì)對(duì)流中的運(yùn)用鮮見(jiàn).熱毛細(xì)對(duì)流與上述流動(dòng)的一個(gè)重要區(qū)別在于邊界條件是速度和溫度強(qiáng)烈耦合,普適的POD 降維方法可能不能直接應(yīng)用.本文工作采用直接數(shù)值模擬方法,獲取二維有限長(zhǎng)液層熱毛細(xì)對(duì)流的流場(chǎng)和溫度場(chǎng)數(shù)據(jù),并利用POD-Galerkin 降維法構(gòu)建低維模型,通過(guò)數(shù)值積分和數(shù)值分岔的計(jì)算,在低維模型上研究流動(dòng)的分岔特性,并與直接數(shù)值模擬的結(jié)果比較,驗(yàn)證模型的準(zhǔn)確性和魯棒性.目的是建立適用于熱毛細(xì)對(duì)流研究的,由直接數(shù)值模擬、POD-Galerkin 降維和數(shù)值分岔三種方法優(yōu)勢(shì)結(jié)合的分岔分析流程.

1 數(shù)值計(jì)算方法

1.1 直接數(shù)值模擬

考慮如圖1 所示體積比(Vr) 分別為0.95,1.00,1.05 的二維有限長(zhǎng)液層中的熱毛細(xì)對(duì)流,其中體積比定義為液層兩端高度固定時(shí),實(shí)際充液體積與界面平直時(shí)液體體積之比.本文假設(shè)流體不可壓,忽略流動(dòng)過(guò)程中界面形狀的變化.

圖1 有限長(zhǎng)二維液層模型Fig.1 Limited liquid layer model

模型的無(wú)量綱參數(shù)為:寬度Ax=3,高度Ay=1,普朗特?cái)?shù)Pr=4.4.

無(wú)量綱控制方程(以體積比為1.00 的情況為例)

其中,V=(u,v),Re為雷諾數(shù).

速度邊界

界面位形通過(guò)求解楊-拉普拉斯方程得到,采用同位網(wǎng)格有限體積法,利用貼體坐標(biāo)變換,在計(jì)算域中離散方程,對(duì)流項(xiàng)用二階中心差分離散,壓力和速度的耦合用SIMPLE 算法,時(shí)間步長(zhǎng)為0.000 1.

1.2 數(shù)值算法驗(yàn)證

采用非均勻網(wǎng)格,網(wǎng)格分布按照如下公式[25]

其中,q為拉伸參數(shù),為均勻的網(wǎng)格點(diǎn).

表1 列舉了Vr=1.00 時(shí)不同網(wǎng)格數(shù)下計(jì)算的臨界雷諾數(shù),當(dāng)網(wǎng)格數(shù)為60 × 20 (非均勻q=2)時(shí),臨界Re與文獻(xiàn)[26]中的誤差在1%左右,所以在之后的計(jì)算中采用的x方向網(wǎng)格數(shù)為60,y方向網(wǎng)格數(shù)為20,參數(shù)q=2.

表1 直接數(shù)值模擬算法驗(yàn)證Table 1 Verification of DNS method

1.3 POD-Galerkin 降維法

1.3.1 本征正交分解

本征正交分解方法利用直接數(shù)值模擬或?qū)嶒?yàn)得到的流場(chǎng)數(shù)據(jù),構(gòu)造一組規(guī)范正交的特征模態(tài),以此反映流動(dòng)的相干結(jié)構(gòu).各個(gè)特征模態(tài)的能量占比不同，一般少數(shù)的幾階模態(tài)就包含了流場(chǎng)的大部分能量,可以反映出主要的流動(dòng)特征.具體分解過(guò)程如下[24]:

由直接數(shù)值模擬的數(shù)據(jù)計(jì)算某一時(shí)段內(nèi)速度場(chǎng)和溫度場(chǎng)的時(shí)間平均量及擾動(dòng)量

利用擾動(dòng)量構(gòu)造快照矩陣

其中,m,n=1,2,···,N,內(nèi)積定義為(a,b)=,D為物理空間,當(dāng)a,b為向量時(shí),點(diǎn)乘表示向量的點(diǎn)乘;當(dāng)a,b為標(biāo)量時(shí),點(diǎn)乘表示普通的乘法.

由于相關(guān)矩陣是實(shí)對(duì)稱(chēng)矩陣,可以求出它的一組規(guī)范正交的特征向量及對(duì)應(yīng)特征值

利用快照矩陣和特征向量構(gòu)造特征模態(tài)

其中

事實(shí)上,這樣得到的特征模態(tài)在前面定義的內(nèi)積意義下,也是規(guī)范正交的.

1.3.2 Galerkin 投影法

Galerkin 投影法是將解用一組規(guī)范正交的基函數(shù)線(xiàn)性表示,代入控制方程,再將方程投影到基函數(shù)上,得到線(xiàn)性組合系數(shù)的常微分方程組.當(dāng)取基函數(shù)為POD 特征模態(tài),就是POD-Galerkin 降維方法,以本文所用模型為例,具體流程如下.

由原方程推出擾動(dòng)方程

將擾動(dòng)量表示為特征模態(tài)的線(xiàn)性組合,代入擾動(dòng)方程

其中,MV,MT分別為截取的速度和溫度的POD 模態(tài)數(shù).

將擾動(dòng)方程投影到各個(gè)模態(tài)上,由正交性,可得以模態(tài)系數(shù)為因變量,時(shí)間 τ 為自變量的常微分方程組,由于特征模態(tài)自然滿(mǎn)足不可壓縮條件,所以這里不考慮連續(xù)性方程.此外,壓力梯度項(xiàng)在模態(tài)上的投影可化為

由于特征模態(tài)在積分區(qū)域內(nèi)滿(mǎn)足連續(xù)性方程,故式(26)第二個(gè)等號(hào)右端第二項(xiàng)等于0,而由于固壁邊界條件為無(wú)滑移,無(wú)穿透,自由面隨時(shí)間不變化,速度特征模態(tài)在各邊界的法向分量均為0,故式(26)第二個(gè)等號(hào)右端第一項(xiàng)也等于0.所以,壓力梯度項(xiàng)的投影為0,從而

此時(shí),方程組的維數(shù)等于截取的模態(tài)數(shù),由此完成了對(duì)原高維方程的降維,可以很容易地求解這個(gè)低維方程,得到模態(tài)系數(shù),重構(gòu)出流場(chǎng),這樣就實(shí)現(xiàn)了流場(chǎng)的快速求解.

1.3.3 修正的POD 降維模型

POD 方法通過(guò)截取前幾階特征模態(tài)實(shí)現(xiàn)降維,忽略了通常起耗散作用的高階模態(tài),由這樣構(gòu)造出的低維方程求得的解往往不穩(wěn)定,隨著時(shí)間的增長(zhǎng),誤差會(huì)逐漸累積.所以,為了提高POD 降維模型的長(zhǎng)期預(yù)測(cè)能力,需要對(duì)模型進(jìn)行修正,常用的修正方法有內(nèi)在穩(wěn)定性方法[27],譜黏性方法[28],渦黏性耗散模型[29]等,這里采用內(nèi)在穩(wěn)定性方法[27,30].

設(shè)原精確方程的動(dòng)力系統(tǒng)和POD 降維系統(tǒng)分別為

具體實(shí)現(xiàn)中,需要把 ε (t) 投影到模態(tài)系數(shù)向量上,將其表示為

其中,M為模態(tài)數(shù),N為快照數(shù),〈 ·,·〉 表示向量的點(diǎn)乘.

最后得到修正的POD 降維方程

其中

1.3.4 速度場(chǎng)和溫度場(chǎng)的耦合處理

需要指出的是,本文中的熱毛細(xì)對(duì)流是速度場(chǎng)和溫度場(chǎng)互相耦合的過(guò)程,因此需要對(duì)上述一般方法做一定修改:

首先,速度和溫度的相關(guān)矩陣不作為獨(dú)立的矩陣處理,而是將兩個(gè)矩陣相加,對(duì)求和后的相關(guān)矩陣取特征向量和特征值

其次,速度和溫度的特征模態(tài)的公式仍如前文所述,但是將速度和溫度表示成特征模態(tài)的線(xiàn)性組合時(shí),設(shè)它們的截?cái)嗄B(tài)數(shù)(M)和組合系數(shù)均相同

這樣處理是考慮到自由面上的邊界條件由速度和溫度耦合表示,其各特征模態(tài)在構(gòu)造低維方程時(shí)也要相對(duì)應(yīng),才能保證滿(mǎn)足邊界條件.

1.4 數(shù)值分岔算法

對(duì)于一個(gè)含參數(shù)p的動(dòng)力系統(tǒng)=F(u,p),可以通過(guò)數(shù)值求解分支解和分岔點(diǎn)來(lái)研究系統(tǒng)的解隨參數(shù)變化的分岔特性.

平衡點(diǎn)滿(mǎn)足方程

一般用牛頓迭代法求解.

周期解滿(mǎn)足方程(其中T為待求的周期)

添加一個(gè)相位的條件,可將其視為兩點(diǎn)邊值問(wèn)題求解.

平衡點(diǎn)和周期解分支可由參數(shù)延續(xù)法(continuation method)得到,各類(lèi)分岔點(diǎn)可通過(guò)對(duì)應(yīng)的方程確定,詳見(jiàn)文獻(xiàn)[13-14].

2 計(jì)算結(jié)果

2.1 直接數(shù)值模擬結(jié)果

利用直接數(shù)值模擬,在不同Re下,計(jì)算了各體積比模型中的熱毛細(xì)對(duì)流.隨著Re的增加,流動(dòng)經(jīng)歷了穩(wěn)態(tài)到振蕩再到穩(wěn)態(tài)的發(fā)展過(guò)程,得到了Rec1,Rec2兩個(gè)臨界Re,如表2 所示

表2 各體積比下的臨界雷諾數(shù)Table 2 Critical Reynolds number under different volume ratios

2.2 POD 降維模型計(jì)算結(jié)果

2.2.1 構(gòu)造修正的降維模型

以Vr=1.00 為例,對(duì)Re=3500 提取振蕩穩(wěn)定后,無(wú)量綱時(shí)間間隔為0.01 的201 個(gè)速度場(chǎng)和溫度場(chǎng)快照,進(jìn)行了本征正交分解,對(duì)特征值按從大到小排序,并計(jì)算前n個(gè)模態(tài)的特征值之和占所有模態(tài)特征值之和的百分比,得到的能量占比隨模態(tài)數(shù)的變化如表3 所示.

表3 不同模態(tài)數(shù)(n)下的能量占比(Re=3500,Vr=1.00)Table 3 The accumulative energy contribution of the first n POD eigenmodes (Re=3500,Vr=1.00)

隨著模態(tài)個(gè)數(shù)的增加,能量占比逐漸接近1,前12 階模態(tài)的能量占比已超過(guò)99.99%,故取前12 階模態(tài)構(gòu)造Re=3500 時(shí)的低維方程,取初值為第一個(gè)快照的模態(tài)系數(shù),對(duì)低維方程數(shù)值積分,得到模態(tài)系數(shù)隨時(shí)間的變化,并將其與原始快照數(shù)據(jù)直接在模態(tài)上投影得到的系數(shù)作對(duì)比.

如圖2(a),未做修正時(shí),隨著時(shí)間的推進(jìn),低維模型得到的第一模態(tài)系數(shù)振蕩的振幅和頻率與DNS數(shù)據(jù)的誤差逐漸增加.而用內(nèi)在穩(wěn)定性方法修正后的低維模型,同樣取初值為第一個(gè)快照的模態(tài)系數(shù),進(jìn)行數(shù)值積分,計(jì)算的結(jié)果與直接數(shù)值模擬的結(jié)果幾乎沒(méi)有差別(圖2(b)).

圖2 DNS 與低維模型得到的第一模態(tài)系數(shù)的時(shí)間序列對(duì)比Fig.2 The comparison of time evolution of first eigenmode coefficient obtained by DNS and low-order model

下面利用修正低維模型進(jìn)行參數(shù)外推,計(jì)算其他Re數(shù) (Re=3400,3300,3200,3100)下的解,并與DNS 計(jì)算結(jié)果對(duì)比(圖3).

圖3 液層中心點(diǎn)溫度振蕩時(shí)間序列對(duì)比Fig.3 The comparison of time evolution of the monitor point at the center of the liquid layer

可以看到,當(dāng)Re數(shù)在3500 附近時(shí),低維模型與DNS 的結(jié)果比較吻合,當(dāng)Re數(shù)逐漸遠(yuǎn)離3500,低維模型得到的振蕩的振幅和頻率與DNS 結(jié)果的誤差逐漸變大.所以,上述方法構(gòu)建的低維模型在一定參數(shù)范圍內(nèi)可以反映原系統(tǒng)的定性和定量特點(diǎn).

為了實(shí)現(xiàn)更大范圍的參數(shù)外推,考慮擴(kuò)大選取的快照數(shù)據(jù)的參數(shù)范圍,即分別取Re=2000,3500,5000 下的201 個(gè)流場(chǎng)和溫度場(chǎng)數(shù)據(jù),無(wú)量綱時(shí)間間隔為0.01,拼接成更大的快照矩陣(共603 列),對(duì)拼接后的快照矩陣取特征模態(tài),構(gòu)建低維方程,這樣可以讓特征模態(tài)包含更大參數(shù)范圍內(nèi)的流動(dòng)特性.簡(jiǎn)便起見(jiàn),后文中將這一組合的情況記為Re=2000,3500,5000.

在此情況下,前20 階模態(tài)的能量占比已超過(guò)99.99%,故取模態(tài)數(shù)為20,對(duì)低維方程分別取Re=2500,3000,4000,4500 進(jìn)行計(jì)算,結(jié)果如圖4 和表4.

表4 低維方程得到的溫度振蕩頻率與DNS 結(jié)果對(duì)比(Vr=1.00)Table 4 The comparison of temperature oscillation frequency obtained by DNS and low-order model

圖4 不同Re 數(shù)下由DNS 和低維模型計(jì)算出的液層中心點(diǎn)溫度振蕩時(shí)間序列、頻譜的對(duì)比 (Vr=1.00)Fig.4 The comparison of time evolution and frequency spectrum of temperature oscillation at the center of the liquid layer obtained by DNS and low-order model for various Re (Vr=1.00)

Re=2500,3000 時(shí),由低維方程計(jì)算出的溫度振幅略小于DNS 得到的溫度振幅,而當(dāng)Re=4000,4500 時(shí),由低維方程計(jì)算出的溫度振幅略大于DNS 得到的溫度振幅.在上述四個(gè)Re處,由低維方程計(jì)算出的振蕩頻率與DNS 得到的振蕩頻率有相似的分布,頻率相對(duì)誤差均在5%左右,說(shuō)明低維方程可以很好地反映原系統(tǒng)的主要振蕩特征.注意到,Re=4500 時(shí),振蕩頻率的相對(duì)誤差超過(guò)了5%,并且低維方程得到的頻譜圖中最高峰對(duì)應(yīng)的頻率與DNS 相差較大.事實(shí)上,由表4 還可以發(fā)現(xiàn),當(dāng)Re=3000,4000 時(shí),頻率誤差較小,而當(dāng)Re接近Rec1和Rec2兩個(gè)臨界值時(shí),誤差增大.這個(gè)現(xiàn)象可能是因?yàn)榭煺站仃囀荝e=2000,3500,5000 下數(shù)據(jù)的拼接,Re=3500 時(shí)的振幅大于Re=2000,Re=5000 的,而特征模態(tài)是根據(jù)能量排序,這就導(dǎo)致Re=3500 的流動(dòng)在特征模態(tài)中占主導(dǎo),從而在Re=3500 附近,低維方程的預(yù)測(cè)效果更好.上述分析表明,取快照矩陣為多個(gè)參數(shù)下數(shù)據(jù)的組合,可以擴(kuò)大低維模型參數(shù)外推的范圍.

2.2.2 在低維模型上計(jì)算分岔

POD-Galerkin 方法降低了原系統(tǒng)的維數(shù),顯著減少了數(shù)值分岔方法的計(jì)算量.對(duì)于Vr=0.95,1.00,1.05 三個(gè)模型,分別構(gòu)建了POD 低維方程,其中模態(tài)數(shù)取為使得前n階模態(tài)能量占比大于99.99%的最小n,并利用低維方程進(jìn)行數(shù)值分岔的計(jì)算,得到如圖5 所示的分岔圖,分支的起點(diǎn)均為Re=1000 下的平衡點(diǎn),其中H 表示Hopf 分岔,黑色為平衡點(diǎn)分支曲線(xiàn),紅色為周期解分支曲線(xiàn),實(shí)線(xiàn)表示穩(wěn)定解,虛線(xiàn)表示不穩(wěn)定解,縱坐標(biāo)定義為

圖5 各體積比下低維方程的分岔圖Fig.5 Bifurcation diagram of reduced-order model under different volume ratios

其中bi(i=1,2,···,M) 為模態(tài)系數(shù).周期解的norm取為一個(gè)周期內(nèi)norm的時(shí)間平均.

對(duì)Vr=0.95,分別取Re=2500,4500,6500 下的201 個(gè)DNS 得到的速度場(chǎng)和溫度場(chǎng)數(shù)據(jù),組合成快照矩陣(603 列),模態(tài)數(shù)為24.對(duì)Vr=1.00,分別取Re=2000,3500,5000 下的201 個(gè)DNS 得到的速度場(chǎng)和溫度場(chǎng)數(shù)據(jù),組合成快照矩陣,模態(tài)數(shù)為20.對(duì)Vr=1.05,分別取Re=1800,2500,3500 下的201 個(gè)DNS 得到的速度場(chǎng)和溫度場(chǎng)數(shù)據(jù),組合成快照矩陣,模態(tài)數(shù)為20.

低維模型的霍普夫分岔點(diǎn),即為起振的臨界值,與DNS 得到的臨界值對(duì)比見(jiàn)表5.

表5 利用低維方程與DNS 計(jì)算的臨界Re 對(duì)比Table 5 The comparison of the critical Reynolds number obtained by DNS and reduced-order model

為了驗(yàn)證低維方程所求周期解的準(zhǔn)確性,以Vr=1.00 為例,將不同Re下周期解的頻率與DNS 的結(jié)果進(jìn)行對(duì)比(表6),可以看到,低維方程上計(jì)算得到的周期解頻率與DNS 的結(jié)果十分接近.

表6 低維方程周期解頻率與DNS 結(jié)果對(duì)比(Vr=1.00)Table 6 The comparison of the frequency of the periodic solution obtained by DNS and reduced-order model (Vr=1.00)

利用低維方程計(jì)算出的分岔圖,在分支計(jì)算和分岔類(lèi)型上,與原系統(tǒng)一致,臨界Re的誤差在15%左右,周期解分支的頻率與原系統(tǒng)振蕩解的頻率較為吻合,說(shuō)明低維模型可以在一定程度上反映原系統(tǒng)解隨參數(shù)變化的轉(zhuǎn)捩過(guò)程.更重要的是,由于流動(dòng)的控制方程是偏微分方程,離散后的維數(shù)通常極高,這使得在原系統(tǒng)上進(jìn)行數(shù)值分岔的計(jì)算十分困難,而利用低維模型計(jì)算分岔,在犧牲一定精度的前提下,極大地提高了計(jì)算效率,低維方程的分岔圖可以為確定流動(dòng)的轉(zhuǎn)捩形式提供初步的指導(dǎo).

2.2.3 對(duì)體積比參數(shù)外推

體積比是液池?zé)崦?xì)對(duì)流分岔過(guò)程的一個(gè)重要臨界參數(shù),下面考慮利用POD 低維方程,對(duì)體積比參數(shù)Vr外推.

用Vr=1.00 時(shí)Re=2000,3500,5000 情況下的特征模態(tài),作為Vr=0.95,1.05 的特征模態(tài),對(duì)Vr=0.95,1.05 分別構(gòu)建低維方程,發(fā)現(xiàn)在一定的體積比和Re范圍內(nèi),低維方程有較好的外推能力,基本反映了流動(dòng)的主要特征(圖6～圖7).

圖6 體積比為0.95 時(shí)不同Re 下由DNS 和低維模型計(jì)算出的液層中心點(diǎn)溫度振蕩時(shí)間序列、頻譜的對(duì)比Fig.6 The comparison of time evolution and frequency spectrum of temperature oscillation at the center of the liquid layer obtained by DNS and low-order model for various Re numbers under Vr=0.95

圖7 體積比為1.05 時(shí)不同Re 下由DNS 和低維模型計(jì)算出的液層中心點(diǎn)溫度振蕩時(shí)間序列、頻譜的對(duì)比Fig.7 The comparison of time evolution and frequency spectrum of temperature oscillation at the center of the liquid layer obtained by DNS and low-order model for various Re numbers under Vr=1.05

3 結(jié)論

作為流動(dòng)與傳熱耦合的非線(xiàn)性過(guò)程,熱毛細(xì)對(duì)流有復(fù)雜的轉(zhuǎn)捩過(guò)程.本文嘗試通過(guò)直接數(shù)值模擬、POD 降維、數(shù)值分岔方法相結(jié)合的手段,高效分析熱毛細(xì)對(duì)流的轉(zhuǎn)捩.本文構(gòu)建了貼體坐標(biāo)系下熱毛細(xì)對(duì)流的POD-Galerkin 低維修正模型,將低維方程的計(jì)算結(jié)果與直接數(shù)值模擬的結(jié)果進(jìn)行對(duì)比,得出以下結(jié)論.

(1) POD-Galerkin 降維方法可以極大地減少所研究系統(tǒng)的維數(shù),流體力學(xué)的控制方程是一組偏微分方程,理論上是無(wú)窮維的,用傳統(tǒng)的數(shù)值方法離散后,維數(shù)通常也非常大,對(duì)于本文中二維的矩形液層,網(wǎng)格數(shù)為1200,有四個(gè)未知量u,v,p,T,如果耦合起來(lái)求解,維數(shù)為4800,通過(guò)POD-Galerkin 降維方法,截取大約20 個(gè)主要模態(tài),將速度與溫度耦合考慮,最后需要求解的常微分方程的維數(shù)約為20.

(2) 對(duì)于熱毛細(xì)對(duì)流這一速度和溫度耦合的模型,POD-Galerkin 降維方法構(gòu)建出的低維方程,可以反映原系統(tǒng)的流動(dòng)特性.進(jìn)一步地,利用內(nèi)在穩(wěn)定性修正方法對(duì)模型進(jìn)行修正,可以提高模型的準(zhǔn)確性;用不同參數(shù)下的直接數(shù)值模擬數(shù)據(jù)拼接成快照矩陣的混合方法可以提高低維方程參數(shù)外推的魯棒性.對(duì)于貼體坐標(biāo)系下的幾何參數(shù),低維方程也在一定參數(shù)范圍內(nèi)有外推能力,可以對(duì)不同物理域內(nèi)的流動(dòng)進(jìn)行快速計(jì)算.

(3) 數(shù)值分岔方法應(yīng)用于低維方程,可以得到與原系統(tǒng)相似的分岔行為,計(jì)算出的周期解振蕩頻率與DNS 計(jì)算出的振蕩頻率相對(duì)誤差在5%以?xún)?nèi),說(shuō)明了在低維模型上用數(shù)值分岔方法來(lái)分析原系統(tǒng)流動(dòng)轉(zhuǎn)捩過(guò)程的可行性.

后續(xù)可能的改進(jìn)及進(jìn)一步工作如下:

(1) 本文中的低維模型從定量角度仍有較大誤差,魯棒性也仍有提高空間;

(2) 本文只是二維液層模型,后續(xù)還可以推廣為擁有更復(fù)雜幾何形狀,更豐富轉(zhuǎn)捩過(guò)程的三維模型;

(3) 本文所分析工況的轉(zhuǎn)捩過(guò)程較為簡(jiǎn)單,PODGalerkin 方法是否可以還原更復(fù)雜的轉(zhuǎn)捩過(guò)程,準(zhǔn)確計(jì)算轉(zhuǎn)捩點(diǎn),得到跨越轉(zhuǎn)捩點(diǎn)前后流動(dòng)的變化特性,仍有待探究.