莫東鳴，徐 敏

(重慶工業(yè)職業(yè)技術(shù)學(xué)院 機(jī)械工程學(xué)院，重慶 401120)

?

隱式重啟動(dòng)Arnoldi迭代法在雙液層流動(dòng)穩(wěn)定性的應(yīng)用

莫東鳴，徐 敏

(重慶工業(yè)職業(yè)技術(shù)學(xué)院 機(jī)械工程學(xué)院，重慶 401120)

通過(guò)引入微小擾動(dòng)量，建立了描述兩不相溶混的上部為固壁的環(huán)形腔內(nèi)雙層流體的熱毛細(xì)對(duì)流的線性擾動(dòng)方程。針對(duì)此復(fù)廣義特征值問(wèn)題，采用隱式重啟動(dòng)Arnoldi迭代法，在單液層流動(dòng)穩(wěn)定性計(jì)算中，對(duì)程序的正確性進(jìn)行了驗(yàn)證, 結(jié)果證明該方法可應(yīng)用于流體流動(dòng)線性穩(wěn)定性分析及節(jié)能計(jì)算。

隱式重啟動(dòng)Arnoldi迭代法；雙液層；穩(wěn)定性

0 引 言

流體動(dòng)力系統(tǒng)的穩(wěn)定性，是一個(gè)古老而現(xiàn)在仍非?；钴S的研究領(lǐng)域，穩(wěn)定性分析涉及到基本定態(tài)解何時(shí)被破壞，如何被破壞，以及后續(xù)的發(fā)展等。穩(wěn)定性分析在工程、氣象、海洋、大氣物理、地球物理學(xué)中都有許多應(yīng)用[1]。不少工程技術(shù)如節(jié)能技術(shù)問(wèn)題已經(jīng)得益于這方面的研究成果。

我們可以通過(guò)某種方法，比如相似性原理、漸近線方法求得微分方程的解，但這個(gè)解所代表的流場(chǎng)、溫度場(chǎng)能否在實(shí)際中實(shí)現(xiàn)，在對(duì)其做穩(wěn)定性分析之前是無(wú)法回答的。如果所得的解在所考慮的參數(shù)空間上是不穩(wěn)定的，那么這個(gè)解實(shí)際上根本無(wú)法實(shí)現(xiàn)，因?yàn)樗蟮玫慕庵淮砹艘环N未受擾動(dòng)情況下的流場(chǎng)、溫度場(chǎng)。而工程實(shí)際中不可避免地會(huì)出現(xiàn)各種形式的擾動(dòng)，如果系統(tǒng)能使各種擾動(dòng)隨時(shí)間衰減下去，至少不增大，即系統(tǒng)是漸進(jìn)穩(wěn)定的或穩(wěn)定的，那么系統(tǒng)所特有的流場(chǎng)、溫度場(chǎng)才可能保持下去；否則，系統(tǒng)所特有的狀態(tài)——基本定態(tài)解，將隨時(shí)間變化，從一種運(yùn)動(dòng)形式變?yōu)榱硪环N運(yùn)動(dòng)形式，比如從定常運(yùn)動(dòng)變?yōu)闀r(shí)相關(guān)運(yùn)動(dòng)，因而所求得的微分方程基本定態(tài)解將毫無(wú)工程意義。

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

求解穩(wěn)定性方程最后可以化為求解一個(gè)特征值問(wèn)題，在實(shí)際工程中所求特征值通常有些特殊的性質(zhì)，例如，為了穩(wěn)定性分析要求k個(gè)實(shí)部最大的特征值，或在復(fù)平面上靠近某個(gè)點(diǎn)的特征值。通常，隨著子空間維數(shù)m的增大，Ritz對(duì)會(huì)逼近特征對(duì)。但實(shí)際計(jì)算中，m可能要非常大時(shí)，Ritz對(duì)才滿足精度要求，可以通過(guò)重啟動(dòng)來(lái)解決這一問(wèn)題。當(dāng)在一個(gè)子空間中用投影類方法求解矩陣的特征值問(wèn)題時(shí)，若該子空間中含有所求特征向量的信息越豐富，則收斂速度越快，這就是選擇“重新啟動(dòng)”方式的原則之一。當(dāng)期望對(duì)特征向量有更多的了解時(shí)，就用特征向量的線性組合更新啟動(dòng)向量，重新啟動(dòng)Arnoldi分解?；舅悸肥钱?dāng)Krylov子空間維數(shù)達(dá)到一設(shè)定值時(shí)，如果Ritz對(duì)仍不收斂，就重新構(gòu)造“啟動(dòng)矢量”。用新的矢量來(lái)重新計(jì)算Arnoldi分解[1]。

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

目前，IRAM因其合適求解大型矩陣特征值問(wèn)題，且有精度高、計(jì)算時(shí)間短的優(yōu)點(diǎn)，在相近領(lǐng)域如電力的動(dòng)力系統(tǒng)穩(wěn)定性分析、核科學(xué)中子輸運(yùn)特性等研究中得以應(yīng)用，但在晶體材料生產(chǎn)過(guò)程中的流動(dòng)穩(wěn)定性分析研究中還未見(jiàn)報(bào)道。

1 物理模型及求解方法

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

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

圖1 物理模型

對(duì)控制方程進(jìn)行無(wú)量綱化，時(shí)間、長(zhǎng)度、速度、溫度和壓力無(wú)量綱尺寸如下： (ro-ri)2/v1,ro-ri,v1/(ro-ri)，ρ1v2/(ro-ri)2，則無(wú)量綱方程可寫為：

▽·Vi=0,

(1)

(2)

(3)

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

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

(4)

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

徑向與縱向速度用無(wú)量綱流函數(shù)ψ定義如下：

(5)

(6)

把方程(6) 代入方程 (1-4) 并消去非線性項(xiàng),則獲得線性擾動(dòng)方程。于是，方程可化簡(jiǎn)為一個(gè)廣義特征值問(wèn)題:

(7)

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

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

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

2 程序正確性驗(yàn)證

為了檢驗(yàn)程序的正確性和網(wǎng)格的收斂性，在環(huán)形單液層與環(huán)形雙液層系統(tǒng)分別做了相關(guān)驗(yàn)證。在環(huán)形單液層系統(tǒng)中，在與文獻(xiàn)[7-8]相同的條件及計(jì)算網(wǎng)格下進(jìn)行了線性穩(wěn)定性計(jì)算，計(jì)算結(jié)果與文獻(xiàn)的結(jié)果非常接近，說(shuō)明程序是正確的。具體結(jié)果分析如下。

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

圖2 單液層微重力情況下發(fā)生熱流體波的臨界Marangoni數(shù)

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

圖3

3 結(jié)束語(yǔ)

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

[1] P.G.德拉津，W.H.雷德，周祖巍，等.流體動(dòng)力穩(wěn)定性[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.斯圖爾特. 矩陣計(jì)算引論[M].上海科學(xué)技術(shù)出版社，1980.

[4] 曹志浩. 矩陣特征值問(wèn)題[M]. 上海科學(xué)技術(shù)出版社，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] 石萬(wàn)元，李友榮，M.K.Ermakov,今石宣之.環(huán)形液層內(nèi)熱毛細(xì)對(duì)流的線性穩(wěn)定性分析[J].工程熱物理學(xué)報(bào),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

重慶市教委科學(xué)技術(shù)研究項(xiàng)目(自然科學(xué)類)(KJ132104)；第二批重慶市高等學(xué)校青年骨干教師資助計(jì)劃(自然科學(xué)類)；重慶工業(yè)職業(yè)技術(shù)學(xué)院科研項(xiàng)目(GZY201313)

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

TP137.7

B

1009-3230(2015)03-0013-05