隱式重啟動Arnoldi迭代法在雙液層流動穩定性的應用

莫東鳴，徐 敏

(重慶工業職業技術學院 機械工程學院，重慶 401120)

?

隱式重啟動Arnoldi迭代法在雙液層流動穩定性的應用

莫東鳴，徐 敏

(重慶工業職業技術學院 機械工程學院，重慶 401120)

通過引入微小擾動量，建立了描述兩不相溶混的上部為固壁的環形腔內雙層流體的熱毛細對流的線性擾動方程。針對此復廣義特征值問題，采用隱式重啟動Arnoldi迭代法，在單液層流動穩定性計算中，對程序的正確性進行了驗證, 結果證明該方法可應用于流體流動線性穩定性分析及節能計算。

隱式重啟動Arnoldi迭代法；雙液層；穩定性

0 引 言

流體動力系統的穩定性，是一個古老而現在仍非常活躍的研究領域，穩定性分析涉及到基本定態解何時被破壞，如何被破壞，以及后續的發展等。穩定性分析在工程、氣象、海洋、大氣物理、地球物理學中都有許多應用[1]。不少工程技術如節能技術問題已經得益于這方面的研究成果。

我們可以通過某種方法，比如相似性原理、漸近線方法求得微分方程的解，但這個解所代表的流場、溫度場能否在實際中實現，在對其做穩定性分析之前是無法回答的。如果所得的解在所考慮的參數空間上是不穩定的，那么這個解實際上根本無法實現，因為所求得的解只代表了一種未受擾動情況下的流場、溫度場。而工程實際中不可避免地會出現各種形式的擾動，如果系統能使各種擾動隨時間衰減下去，至少不增大，即系統是漸進穩定的或穩定的，那么系統所特有的流場、溫度場才可能保持下去；否則，系統所特有的狀態——基本定態解，將隨時間變化，從一種運動形式變為另一種運動形式，比如從定常運動變為時相關運動，因而所求得的微分方程基本定態解將毫無工程意義。

在晶體材料生長過程中，最常使用的制備方法為提拉法，為了抑制熔體的蒸發和超過臨界溫度之后出現熱流體波對晶體質量的破壞，可以在熔體上方施加一層不相溶混的液體，此種生長晶體的方法稱為液封提拉法。由于增加了液封層，使得熔體流動穩定性的研究變得復雜。在文章[2]，我們已求得了具有液封環形液池內熱對流的基態解，如果對求得的這個解的穩定性分析的結果能夠表明在很大的控制參數(Marangoni數) 范圍內解是穩定的，那么說明,這個解在很大程度上能夠保持下去，即這種被削弱了的運動形式能夠維持下去，那就說明通過液封削弱熱毛細對流的方法在理論上是可行的；反之，如果穩定性分析的結果是：無論多小的控制參數，解都是不穩定的。由此可見，對所關心的熱物理系統的基本定態解做穩定性分析，顯得非常必要和特別重要。在晶體生長過程中，大量理論分析和實驗研究都重視熔區對流的穩定性對生長晶體的質量將帶來直接的影響。遺憾的是對具有液封的環形雙液層的熱毛細對流的穩定性分析還未見文獻報導。本文對之做穩定性分析無疑可以彌補這方面理論的不足，同時可望對工程實際有直接的指導意義，為有關參數的選擇提供重要依據。

求解穩定性方程最后可以化為求解一個特征值問題，在實際工程中所求特征值通常有些特殊的性質，例如，為了穩定性分析要求k個實部最大的特征值，或在復平面上靠近某個點的特征值。通常，隨著子空間維數m的增大，Ritz對會逼近特征對。但實際計算中，m可能要非常大時，Ritz對才滿足精度要求，可以通過重啟動來解決這一問題。當在一個子空間中用投影類方法求解矩陣的特征值問題時，若該子空間中含有所求特征向量的信息越豐富，則收斂速度越快，這就是選擇“重新啟動”方式的原則之一。當期望對特征向量有更多的了解時，就用特征向量的線性組合更新啟動向量，重新啟動Arnoldi分解。基本思路是當Krylov子空間維數達到一設定值時，如果Ritz對仍不收斂，就重新構造“啟動矢量”。用新的矢量來重新計算Arnoldi分解[1]。

由 1992年Soresen提出的隱式重啟動方法，被公認為是最有效的重開始技術之一。它是一種基于隱式QR位移分解的簡便和穩定的方法，不用顯示計算Arnoldi分解。IRAM算法的優勢在于：把一個大型的稀疏矩陣(帶狀矩陣，階數為n)正交投影為一個遠小于其矩陣階數的小上Hessenberg矩陣(階數為k)，這一過程體現了空間節省的優點；通過求此小矩陣的特征值與特征向量，再經過重啟動迭代技術，收斂得到近似于原矩陣的特定特征值，這一過程體現了空間節省、提高收斂精度、縮短收斂時間的優點[3]。

目前，IRAM因其合適求解大型矩陣特征值問題，且有精度高、計算時間短的優點，在相近領域如電力的動力系統穩定性分析、核科學中子輸運特性等研究中得以應用，但在晶體材料生產過程中的流動穩定性分析研究中還未見報道。

1 物理模型及求解方法

環形腔中盛載兩不相容混的熔體5cSt silicone oil 與 HT-70，物性參數同文獻[4]。雙層液體施以水平溫度梯度，上固壁溫度Tc,下固壁溫度Th(Tc

在模型中引入如下假設：(1)熔體和液封均為不可壓縮的牛頓流體，滿足Boussinesq近似；(2)流速較低，流動為軸對稱二維層流；(3)兩相界面平整無變形，在液液界面和自由表面均考慮熱毛細力的作用，固-液界面滿足無滑移條件；(4)液池頂部和底部邊界均絕熱；(5)表面張力是溫度的線性函數。

圖1 物理模型

對控制方程進行無量綱化，時間、長度、速度、溫度和壓力無量綱尺寸如下： (ro-ri)2/v1,ro-ri,v1/(ro-ri)，ρ1v2/(ro-ri)2，則無量綱方程可寫為：

▽·Vi=0,

(1)

(2)

(3)

式中：ρi，ai,βi和vi分別為i(i=1, 2)層液體的密度,熱導率,熱擴散系數和動力粘度；Vi代表無量綱速度矢量；Θi(Ti-Tc)/(Th-Tc)為無因次溫度;Pi為無量綱壓力；τ為無量綱時間；ez為垂直方向的速度矢量；Pr=v1/a1為熔體的普朗特數。

在Z=H1=h1/(ro-ri)的液液界面邊界條件為：

(4)

式中：R與Z為無量綱坐標；Ma是Marangoni數，其定義式為Ma=γ2-1(Th-Tc)(ro-ri)/(μ1a1)；μ為動力粘度；κ為導熱率；γT=-?γ/?T是界面張力梯度。

徑向與縱向速度用無量綱流函數ψ定義如下：

(5)

(6)

把方程(6) 代入方程 (1-4) 并消去非線性項,則獲得線性擾動方程。于是，方程可化簡為一個廣義特征值問題:

(7)

式中：A、B為具有復元素的復矩陣；n為節點數；x為特征向量，包含所取節點上的未知速度、溫度和壓力值的特征向量。求解方程(7)的特征值，當特征值λ為零時，可得到臨界 Marangoni數Mac。至此，將穩定性分析問題轉化成為一個復廣義特征值問題。

穩定性分析的根本任務是要在參數Pr，m，μ，α，γT等取定后，確定特征值第一次跨過虛軸所對應的Marangoni數，通過一連續增加Ma的過程，實現臨界Marangoni數的確定。即在μ*,κ*,ρ*,ν*,α*,Pr等取定后，取一足夠小的Ma0作為初值，令Ma=Ma0。計算系統實部最大特征值λmax，視其實部Real(λmax) 是否大于0，并增加Ma重新計算；如此記錄特征值，在最大特征值為零的交界處便可以得到臨界Marangoni數。

顯然，我們并不計算全部特征值，而是計算所關心的主導特征值。由于上述矩陣為帶狀矩陣，所以選用IRAM法。ARPACK軟件包是基于隱式重啟動Arnoldi方法的Fortran77子程序的集合。該軟件使用n·O(k)+O(k2)的存儲開銷計算滿足用戶要求的k個特征值，這包括具有最大實部、最大虛部或最大模的特征值。因為低存儲和計算需求，這個技術適合大規模特征問題[5-6]。

2 程序正確性驗證

為了檢驗程序的正確性和網格的收斂性，在環形單液層與環形雙液層系統分別做了相關驗證。在環形單液層系統中，在與文獻[7-8]相同的條件及計算網格下進行了線性穩定性計算，計算結果與文獻的結果非常接近，說明程序是正確的。具體結果分析如下。

對于ro=40 mm和ri=40 mm，在微重力情況下的環形淺液池進行了線性穩定性分析，物性值取與文獻[8]完全相同。Shi等人在二維數值模擬中，取網格(122-542)R×(16-28)Z，在三維數值模擬中，取網格(82-202)R×(123-603)θ×16Z，故在線性穩定性分析計算中取網格300R×60Z。圖2為文中線性穩定性分析臨界Marangoni數結果與三維數值模擬、實驗、其它線性穩定性分析結果的對比。以1 mm厚度為例，Shi等人的三維數值模擬的Mac=1.68×105,mc=27,Emakov等人的線性穩定性分析Mac=1.68×105,mc=28,本文線性穩定性分析計算Mac=1.608×105,mc=28, 相對誤差為4.29%，計算結果與他們的數值模擬與理論結果非常接近，說明程序是正確的。

圖2 單液層微重力情況下發生熱流體波的臨界Marangoni數

圖3給出了1 mm厚度下，Mac=1.61×105，線性穩定性分析確定的自由表面熱流體波形態比較結果，此時熱流體波運動方向與溫度梯度方向間夾角約為58°，如圖3(a)所示。Shi等人的線性穩定性分析結果和三維數值模擬結果分別如圖3(b)和(c)所示，可以看出波數相同，傳播角相近，從而說明程序在計算特征向量時也是可靠的。

圖3

3 結束語

首先通過引入微小擾動量，建立了描述兩不相溶混的上部為固壁的環形腔與有自由表面的環形池內雙層流體的熱毛細對流的線性擾動方程。在此基礎上，用有限差分方法離散線性擾動方程，經離散化處理，得到了離散化方程組及邊界條件，該方程組構成了一個復廣義特征值問題;針對問題的特點采用隱式重啟動Arnoldi迭代法，以保持矩陣的帶型結構，通過編程計算求解；在單液層系統中對程序的正確性及網格的獨立性都進行了驗證，程序證明利用IRAM解決環形池內流體系統穩定性具有適合大型帶狀矩陣特征值問題、計算快速、可求解多重特征值的優點。后續可在單液層穩定性分析的基礎上進行雙液層流動穩定性問題的計算。

[1] P.G.德拉津，W.H.雷德，周祖巍，等.流體動力穩定性[M].宇航出版社，1990.

[2] Li Y R, Zhang W J, Peng L. Thermal convection in an annular two-layer system under combined action of buoyancy and thermocapillary forces [J]. J. Supercond. Nov. Magn., 2010, 23: 1219-1223.

[3] G.W.斯圖爾特. 矩陣計算引論[M].上海科學技術出版社，1980.

[4] 曹志浩. 矩陣特征值問題[M]. 上海科學技術出版社，1980.

[5] Lehoucq RB. Implicitly Restarted Arnoldi Methods and Subspace Iteration[J]. SIAM Journal on Matrix Analysis and Applications, 2001,23(2):551-562.

[6] Lehoucq RB. Sorensen D C, Vu etal P A. ARPACK:Fortran subroutines for solving large scale engenvalue problems[R]. Release 2.1.

[7] Peng L, Li Y R, Shi W Y, .Three-dimensional thermocapillary-buoyancy flow of silicone oil in a differentially heated annular pool [J]. Int.J.Heat Mass Transfer, 2007, 50(5-6):872-880.

[8] Shi W Y， Ermakov M K, Imaishi N.Effect of pool rotation on thermocapillary convection in shallow annular pool of silicone oil[J]. J.Crystal Growth, 2006(294):474-485.

[9] 石萬元，李友榮，M.K.Ermakov,今石宣之.環形液層內熱毛細對流的線性穩定性分析[J].工程熱物理學報,2008,29(7): 1218-1220.

Application of Implicit Restarted Arnoldi Method on Flow Stability of Two-layer Liquid System

MO Dong-ming,XU Min

(Department of Mechanical Engineering, Chongqing Industry Polytechnic College, Chongqing 401120, China)

By introducing a small disturbance quantity, the linear perturbation equation described thermocapillary flow in a two-layer system with upper rigid wall was set up. For the generalized eigenvalue problem, the Arnoldi restart implicit iteration method was adopted. This method was verified in a single layer system in the correctness of the program, the results proved that it can be applied to linear stability analysis and energy saving calculation.

Implicit restarted; Arnoldi method;Two-layer; Lquid system; Flow stability

10.3969/j.issn.1009-3230.2015.03.003

2014-02-10

2015-02-27

重慶市教委科學技術研究項目(自然科學類)(KJ132104)；第二批重慶市高等學校青年骨干教師資助計劃(自然科學類)；重慶工業職業技術學院科研項目(GZY201313)

莫東鳴(1982-)，女，廣西南寧人，博士，講師，重慶大學動力工程與工程熱物理專業，從事晶體材料制備中熱物理問題的研究。

TP137.7

B

1009-3230(2015)03-0013-05