基于POD-Galerkin 降維方法的熱毛細對流分岔分析1)

郭子漪 趙建福 李 凱 胡文瑞

(中國科學院力學研究所微重力重點實驗室,北京 100190)

(中國科學院大學工程科學學院,北京 100049)

引言

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

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

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

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

1 數值計算方法

1.1 直接數值模擬

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

圖1 有限長二維液層模型Fig.1 Limited liquid layer model

模型的無量綱參數為:寬度Ax=3,高度Ay=1,普朗特數Pr=4.4.

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

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

速度邊界

界面位形通過求解楊-拉普拉斯方程得到,采用同位網格有限體積法,利用貼體坐標變換,在計算域中離散方程,對流項用二階中心差分離散,壓力和速度的耦合用SIMPLE 算法,時間步長為0.000 1.

1.2 數值算法驗證

采用非均勻網格,網格分布按照如下公式[25]

其中,q為拉伸參數,為均勻的網格點.

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

表1 直接數值模擬算法驗證Table 1 Verification of DNS method

1.3 POD-Galerkin 降維法

1.3.1 本征正交分解

本征正交分解方法利用直接數值模擬或實驗得到的流場數據,構造一組規范正交的特征模態,以此反映流動的相干結構.各個特征模態的能量占比不同，一般少數的幾階模態就包含了流場的大部分能量,可以反映出主要的流動特征.具體分解過程如下[24]:

由直接數值模擬的數據計算某一時段內速度場和溫度場的時間平均量及擾動量

利用擾動量構造快照矩陣

其中,m,n=1,2,···,N,內積定義為(a,b)=,D為物理空間,當a,b為向量時,點乘表示向量的點乘;當a,b為標量時,點乘表示普通的乘法.

由于相關矩陣是實對稱矩陣,可以求出它的一組規范正交的特征向量及對應特征值

利用快照矩陣和特征向量構造特征模態

其中

事實上,這樣得到的特征模態在前面定義的內積意義下,也是規范正交的.

1.3.2 Galerkin 投影法

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

由原方程推出擾動方程

將擾動量表示為特征模態的線性組合,代入擾動方程

其中,MV,MT分別為截取的速度和溫度的POD 模態數.

將擾動方程投影到各個模態上,由正交性,可得以模態系數為因變量,時間 τ 為自變量的常微分方程組,由于特征模態自然滿足不可壓縮條件,所以這里不考慮連續性方程.此外,壓力梯度項在模態上的投影可化為

由于特征模態在積分區域內滿足連續性方程,故式(26)第二個等號右端第二項等于0,而由于固壁邊界條件為無滑移,無穿透,自由面隨時間不變化,速度特征模態在各邊界的法向分量均為0,故式(26)第二個等號右端第一項也等于0.所以,壓力梯度項的投影為0,從而

此時,方程組的維數等于截取的模態數,由此完成了對原高維方程的降維,可以很容易地求解這個低維方程,得到模態系數,重構出流場,這樣就實現了流場的快速求解.

1.3.3 修正的POD 降維模型

POD 方法通過截取前幾階特征模態實現降維,忽略了通常起耗散作用的高階模態,由這樣構造出的低維方程求得的解往往不穩定,隨著時間的增長,誤差會逐漸累積.所以,為了提高POD 降維模型的長期預測能力,需要對模型進行修正,常用的修正方法有內在穩定性方法[27],譜黏性方法[28],渦黏性耗散模型[29]等,這里采用內在穩定性方法[27,30].

設原精確方程的動力系統和POD 降維系統分別為

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

其中,M為模態數,N為快照數,〈 ·,·〉 表示向量的點乘.

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

其中

1.3.4 速度場和溫度場的耦合處理

需要指出的是,本文中的熱毛細對流是速度場和溫度場互相耦合的過程,因此需要對上述一般方法做一定修改:

首先,速度和溫度的相關矩陣不作為獨立的矩陣處理,而是將兩個矩陣相加,對求和后的相關矩陣取特征向量和特征值

其次,速度和溫度的特征模態的公式仍如前文所述,但是將速度和溫度表示成特征模態的線性組合時,設它們的截斷模態數(M)和組合系數均相同

這樣處理是考慮到自由面上的邊界條件由速度和溫度耦合表示,其各特征模態在構造低維方程時也要相對應,才能保證滿足邊界條件.

1.4 數值分岔算法

對于一個含參數p的動力系統=F(u,p),可以通過數值求解分支解和分岔點來研究系統的解隨參數變化的分岔特性.

平衡點滿足方程

一般用牛頓迭代法求解.

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

添加一個相位的條件,可將其視為兩點邊值問題求解.

平衡點和周期解分支可由參數延續法(continuation method)得到,各類分岔點可通過對應的方程確定,詳見文獻[13-14].

2 計算結果

2.1 直接數值模擬結果

利用直接數值模擬,在不同Re下,計算了各體積比模型中的熱毛細對流.隨著Re的增加,流動經歷了穩態到振蕩再到穩態的發展過程,得到了Rec1,Rec2兩個臨界Re,如表2 所示

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

2.2 POD 降維模型計算結果

2.2.1 構造修正的降維模型

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

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

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

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

圖2 DNS 與低維模型得到的第一模態系數的時間序列對比Fig.2 The comparison of time evolution of first eigenmode coefficient obtained by DNS and low-order model

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

圖3 液層中心點溫度振蕩時間序列對比Fig.3 The comparison of time evolution of the monitor point at the center of the liquid layer

可以看到,當Re數在3500 附近時,低維模型與DNS 的結果比較吻合,當Re數逐漸遠離3500,低維模型得到的振蕩的振幅和頻率與DNS 結果的誤差逐漸變大.所以,上述方法構建的低維模型在一定參數范圍內可以反映原系統的定性和定量特點.

為了實現更大范圍的參數外推,考慮擴大選取的快照數據的參數范圍,即分別取Re=2000,3500,5000 下的201 個流場和溫度場數據,無量綱時間間隔為0.01,拼接成更大的快照矩陣(共603 列),對拼接后的快照矩陣取特征模態,構建低維方程,這樣可以讓特征模態包含更大參數范圍內的流動特性.簡便起見,后文中將這一組合的情況記為Re=2000,3500,5000.

在此情況下,前20 階模態的能量占比已超過99.99%,故取模態數為20,對低維方程分別取Re=2500,3000,4000,4500 進行計算,結果如圖4 和表4.

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

圖4 不同Re 數下由DNS 和低維模型計算出的液層中心點溫度振蕩時間序列、頻譜的對比 (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 時,由低維方程計算出的溫度振幅略小于DNS 得到的溫度振幅,而當Re=4000,4500 時,由低維方程計算出的溫度振幅略大于DNS 得到的溫度振幅.在上述四個Re處,由低維方程計算出的振蕩頻率與DNS 得到的振蕩頻率有相似的分布,頻率相對誤差均在5%左右,說明低維方程可以很好地反映原系統的主要振蕩特征.注意到,Re=4500 時,振蕩頻率的相對誤差超過了5%,并且低維方程得到的頻譜圖中最高峰對應的頻率與DNS 相差較大.事實上,由表4 還可以發現,當Re=3000,4000 時,頻率誤差較小,而當Re接近Rec1和Rec2兩個臨界值時,誤差增大.這個現象可能是因為快照矩陣是Re=2000,3500,5000 下數據的拼接,Re=3500 時的振幅大于Re=2000,Re=5000 的,而特征模態是根據能量排序,這就導致Re=3500 的流動在特征模態中占主導,從而在Re=3500 附近,低維方程的預測效果更好.上述分析表明,取快照矩陣為多個參數下數據的組合,可以擴大低維模型參數外推的范圍.

2.2.2 在低維模型上計算分岔

POD-Galerkin 方法降低了原系統的維數,顯著減少了數值分岔方法的計算量.對于Vr=0.95,1.00,1.05 三個模型,分別構建了POD 低維方程,其中模態數取為使得前n階模態能量占比大于99.99%的最小n,并利用低維方程進行數值分岔的計算,得到如圖5 所示的分岔圖,分支的起點均為Re=1000 下的平衡點,其中H 表示Hopf 分岔,黑色為平衡點分支曲線,紅色為周期解分支曲線,實線表示穩定解,虛線表示不穩定解,縱坐標定義為

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

其中bi(i=1,2,···,M) 為模態系數.周期解的norm取為一個周期內norm的時間平均.

對Vr=0.95,分別取Re=2500,4500,6500 下的201 個DNS 得到的速度場和溫度場數據,組合成快照矩陣(603 列),模態數為24.對Vr=1.00,分別取Re=2000,3500,5000 下的201 個DNS 得到的速度場和溫度場數據,組合成快照矩陣,模態數為20.對Vr=1.05,分別取Re=1800,2500,3500 下的201 個DNS 得到的速度場和溫度場數據,組合成快照矩陣,模態數為20.

低維模型的霍普夫分岔點,即為起振的臨界值,與DNS 得到的臨界值對比見表5.

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

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

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

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

2.2.3 對體積比參數外推

體積比是液池熱毛細對流分岔過程的一個重要臨界參數,下面考慮利用POD 低維方程,對體積比參數Vr外推.

用Vr=1.00 時Re=2000,3500,5000 情況下的特征模態,作為Vr=0.95,1.05 的特征模態,對Vr=0.95,1.05 分別構建低維方程,發現在一定的體積比和Re范圍內,低維方程有較好的外推能力,基本反映了流動的主要特征(圖6～圖7).

圖6 體積比為0.95 時不同Re 下由DNS 和低維模型計算出的液層中心點溫度振蕩時間序列、頻譜的對比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 時不同Re 下由DNS 和低維模型計算出的液層中心點溫度振蕩時間序列、頻譜的對比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 結論

作為流動與傳熱耦合的非線性過程,熱毛細對流有復雜的轉捩過程.本文嘗試通過直接數值模擬、POD 降維、數值分岔方法相結合的手段,高效分析熱毛細對流的轉捩.本文構建了貼體坐標系下熱毛細對流的POD-Galerkin 低維修正模型,將低維方程的計算結果與直接數值模擬的結果進行對比,得出以下結論.

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

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

(3) 數值分岔方法應用于低維方程,可以得到與原系統相似的分岔行為,計算出的周期解振蕩頻率與DNS 計算出的振蕩頻率相對誤差在5%以內,說明了在低維模型上用數值分岔方法來分析原系統流動轉捩過程的可行性.

后續可能的改進及進一步工作如下:

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

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

(3) 本文所分析工況的轉捩過程較為簡單,PODGalerkin 方法是否可以還原更復雜的轉捩過程,準確計算轉捩點,得到跨越轉捩點前后流動的變化特性,仍有待探究.