9存在相變的自由表面流體問題模擬

1、第九章存在相變的自由表面流動問題模擬一一三相流T.L, MARINDeparme nt of Mining Engin eeri ng , Uni versity of Chile,Av. Tupper 2069 , San tiago, ChileE-mail:tmari ni ng.uchile.cl本章介紹了流體自由表面凝固過程的固定網(wǎng)格數(shù)值模型，并采用水平集方法來描述自由表面的流動。Voller和Prakash對此方法進行了改進，考慮了液相凝 固問題，包括對流和傳遞過程。在該方法中，液體物性與溫度相關，通過設定0到1之間的孔隙率使得體積力依賴于孔隙率和溫度的變化，并通過對 Navier

2、-Stokes方程的修正來模擬液相或者固相。此外，在傳熱方程中使用改進的 熱容表達式來計算融化潛熱。1 .緒論液體在凝固時主要通過對流來傳熱，由于涉及到移動界面問題，所以通常難于模擬。此外，對流體有效的控制方程對固體不再適用。而且，如果問題中包含 自由表面，那么流動將變得更加復雜，需要考慮相變或者凝固問題，這時都需要 跟蹤相界面的變化。目前有很多種處理液體凝固的計算方法，考慮封閉空間中存在導熱和自然對 流的情況，Voller和Prakash對此作了總結1。進展之一為采用變形網(wǎng)格來處理 液體一固體界面的移動問題。也可用固定網(wǎng)格以及焓變隨溫度的變化來處理該問 題，這種方法的特點是較為簡便，但在研究

3、零速度封閉容器中的固體時會存在問 題。對于給定的計算單元，可以簡單的設置速度為零，或者通過設置粘度為潛熱 容2函數(shù)的方法來實現(xiàn)3。在這類情況下，當潛熱容接近于零時，粘性會很大， 這樣才能模擬固體物性。Voller和Prakash1研究了在一定溫度范圍內發(fā)生相變 的情況，他們將流體描述為多孔介質，通過在Navier-Stokes方程中引入源項來模擬相變過程，用于研究速度為零的固體情況。存在相變的自由表面流動問題的研究難點在于其實驗和計算驗證均很困難。Pasandideh-Fard等4采用流體體積元方法(VOF)對錫液滴和水滴熱表面上5的 凝固進行了實驗和計算研究，對自由表面及過程中的焓變進行跟蹤

4、。他們的模型引入了液滴接觸角以及液滴底部界面?zhèn)鳠嵯禂?shù)，該系數(shù)由實驗測試得到并作為模擬的邊界條件。2.控制方程純凈流體凝固過程由熔點決定。在液體冷卻過程中，一旦達到這個溫度，在 溫度繼續(xù)降低之前液體開始釋放相變潛熱。但是在多組分體系中，相變存在一個溫度范圍，從固體開始出現(xiàn)時的液體溫度開始，一直到最后一種液體凝固時的溫 度為止。在這種情況下，融化潛熱在溫度改變的同時不斷釋放。相變期間固體組分(Fs)可以表示為溫度的函數(shù)：0TA(TmZ)T +名一TFs T 二Tm ； T _ © _ ；(1)21T ： ： Tm - ；&為液相線和固相線其中，T為系統(tǒng)溫度，Tm為液相線和固相線的

5、溫度平均值, 溫度差值的一半。因此，固相和液相溫度為：Ts 二 Tl =Tm ；前面已經(jīng)提到，總的體系熱容 H由兩部分組成，顯熱h和潛熱H。顯熱通過以 下方程計算：h 二 Cp T(4)潛熱可以表示為溫度的函數(shù)，根據(jù)之前對固體比例的定義：LT 一 Tm ；H T = L FsTm ； T _ Tm - ；oT<(Tm-E)導熱和對流傳熱方程用系統(tǒng)溫度形式表示為：? Cp ' - -k T CpTu 二 0(6):t在這種情況下，釋放潛熱帶來的影響可以通過對有效熱容方程的重新定義來包含 到熱容項中：CpT 二 Cp - CpH(7)CpH 二L、T : L fldc2hsT-Tm,

6、；(8)這里，”(T)表示平滑delta方程，在 COMSOL Multiphysics中由fldc2hs函數(shù)建 立。注意，這里需要對該函數(shù)在整個溫度范圍內積分，但是只有T =Tm-；到T =Tm ；的溫度范圍。只有在模擬純凈物質凝固過程的時候，；才為零，此時(T)將變?yōu)檎嬲腄irac 'Delta。通過這種方法，系統(tǒng)總熱容可以通過對方程 (7)積分得到：T?T?T2 _.H =CpT dT =Cp dT L T dT(9)T1T111速度場和壓力場的控制方程即 Navier-Stokes方程：utu u i 亠門 u u 、p = F(10)-tu=0(11)體積力F包括重力和水平

7、集方法處理的表面張力，同時F也包含依賴于凝固過程的固體比例源項。體積力分量如下：FxSx(12)ex(13)FygSy&y其中，匚是液體表面張力，I是交界面曲率， '是水平集函數(shù)，g是重力加速度， Sx和Sy項表述如下。為了將固化過程并入體積力中，可以將液相看作是一種多孔介質， 孔隙率' 依賴于溫度。全液相對應于孔隙率為1的狀態(tài)，全固相對應于孔隙率為0的狀態(tài)。 孔隙率定義如下：1T 啟(Tm + E)T T =1-Fs T2 - Tm ； T 一 Tm-；(14)0T ： ： Tm - ；源項定義如下：Sx _ -Au(15)(16)Sy - - Av2這里A是孔隙率的

8、函數(shù)，根據(jù)以下方程定義:(17)下面介紹一下這些源項的作用。當溫度位于液相線以上時，系統(tǒng)處于全液態(tài)，源 項取值為零，對Navier-Stokes方程沒有任何改動。多孔區(qū)域意味著溫度處于液 相線和固相線之間，A值增大并影響瞬態(tài)、對流和擴散項，流動方程近似相當于 多孔介質中的Darcy定律。當溫度進一步降低，孔隙率接近于 0時(固相)，該 源項決定了所有其它源項，使得速度值趨于零(固相)。方程(17)中的常數(shù)C和q 是任意選取的，依賴于求解的具體問題，通常C取較大值而q取較小值，避免當變?yōu)榱銜r方程被0除。水平集函數(shù) '定義了兩個初態(tài)相(例如，氣和液)的交界面。但是當固化(或 者反問題中融化

9、)時會出現(xiàn)第三相，一般出現(xiàn)在兩個初始相的其中之一， 所以需 要修改源項Sx和Sy，使其嚴格作用在某一相中。求解該問題非常簡單，只需要將這幾個源項乘以 '的Heaviside函數(shù)即可，出于我們對相的選擇，我們限制了(18)這些源項的作用范圍。這樣，方程(10)中的體積力最終可以表示為:FxSxHexFygSyHcy類似水平集方法，Heaviside函數(shù)同樣用來定義連續(xù)相和非連續(xù)相間劇烈的物性 變化。這些物性包括密度，粘度和熱傳導率。(19)3.結果與討論該模擬過程需要三個物理模型：不可壓 Navier-Stokes流動，對流和擴散傳 質過程，對流和熱傳導的傳熱過程。求解該問題時用到了 C

10、OMSOL Multiphysics 的多物理場功能來同時處理和求解三個應用模式。3.1下降液滴的固化本節(jié)介紹一個液滴在冷表面的流動和固化過程。為簡單起見，將其簡化為二維模型。整個系統(tǒng)用左下角在 x=y=0處的無量綱矩形域來表示，底面邊界溫度 固定，其余邊界絕熱。液體和連續(xù)相（環(huán)境）的初始邊界條件用以下初始水平集 函數(shù)來表示：t = 0 二.x2 y2 一0.5（20）這樣就在模型數(shù)值域的左下角形成四分之一個圓。注意這種情況下液相值為負，連續(xù)相 '值為正。為簡化起見，在該問題中我們不考慮液/氣交界面處表面張力的影響，集中 處理固化問題。但是就像在水平集方法中提到的，加入這些影響因素也非

11、常簡單。打開COMSOL Multiphysics模型導航欄。按照表1中的步驟建立模型。注意 Navier-Stokes模型中邊界2 （底邊）為滑移/對稱邊界條件。雖然正常邊界條件是 無滑移的，但是該假設對處理邊界上水平集“phi”函數(shù)的對流非常有幫助，更容易模擬液滴和底面邊界的接觸過程。同時，由于液滴從底面向上固化，速度自然會計算為零，由于對固相進行求解，所以在這種情況下滑移/對稱邊界條件并不算太脫離現(xiàn)實。如果使用無滑移邊界條件，就需要用到水平集函數(shù)的再次初始 化。表1水平集三相固化模型模型導航欄選擇2D維數(shù)COMSOL Multiphysics| Fluid Flow| Incompres

12、sible Navier- Stokes點擊Multiphysics選項卡，添加。選擇 Convection and Diffusion。設疋因變量為phi，添加。Heat Transfer| Conv ecti on and Con duction 添加，完成。Draw菜單Specify Objects| Square輸入width = 1，基準角x=y= 0。完成Opti ons 菜單：Con sta ntsX0=y0=0, r=0.5, rhog=1, rhol=10, etag=etal=1, gy=-10, n=0.02, Ti=1, TC=-0.2, e=0.1, cp=1, kl

13、=0.5, kg=0.01, L=1, C=1600, q=0.001, tm=0Opti ons 菜單： Expressi onsphi0=sqrt(x-x0)A2+(y-yj)A2)-r Hphi=(1+ta nh(-phi/n)/2 rho=rhog+(rhol-rhog)*Hphi eta=etag+(etal-etag)*Hphi k=kg+ (kl-kg)*Hphilambd=(T-Tm+e)/(2* e)*(T<= (Tm+e)*(T>= (Tm-e)+(T>(Tm+e) Fs=1-lambdA=-C*(1-lambd)A2/(lambdA3+ q)Sx=-A*

14、u, Sy=-A*vCpH=L*fldc2hs(T-Tm,e)CpT=Cp+CpHPhysics 菜單： Point sett ingsns模式。確認3點的點約束設定為0壓力情況。Physics 菜單： Boun dary sett ingsns模式：邊界 1 和 2： slip/symmetry，邊界 3 和 4： no-slip chcd模式： 設疋所有邊界條件為 In sulatio n/symmetry cc模式：邊界2設定為溫度，在溫度欄中輸入 TC設定邊界1，3和4為絕熱邊界。Physics 菜單： Subdoma in Sett ingsns模式：尸rho, r=eta, Fx=

15、-Sx*H phi, Fy=-Sy*H phi+rho*gy chcd 模式：D(各向同性)=0, R=0, u=u, v=vInit 選項卡:phi(t0)=phi0cc模式：k(各向同性)=k, p=0, Cp=CpT, Q=0, u=u, v=v。 Init選項卡:T(t0)=Ti。完成Mesh菜單： Map ped Mesh選擇邊界 1，選中“ Constrained edge element distribution' 選 項，輸入40各邊界基元。選擇邊界 2，選中“ Constrained edge element distribution' 選 項，輸入40各邊界

16、基元。重繪網(wǎng)格，完成。Solve菜單選擇時間依賴求解器。設定輸出時間0:0.02:2,點擊工具欄上 的求解按鈕(=)。根據(jù)方程(20),初始水平集函數(shù)“ phiO”用空間坐標“ x”和“y”來定義， 液相為正(或者在液相區(qū)域)、氣相為負，如圖1所示。Heaviside函數(shù)“ Hphi” 用phi的負值來定義，所以在液滴內部為1，其它任何區(qū)域均為0。液相和氣相 的物性用取決于水平集函數(shù)“ phi ”的Heaviside函數(shù)來定義。根據(jù)方程(14)，用 邏輯表達式定義孔隙率函數(shù)，該變量依賴于溫度變量“T”和固相線、液相線的溫度。然后用該函數(shù)定義其它依賴于溫度的物性。ime-U XFect Ccre

17、rtnaticn, chCrnomn 3!圖1水平集函數(shù)初始條件的表面圖該問題的新特性在于使用了映射網(wǎng)格。它提供了棋盤形分布的四方形網(wǎng)格。圖2給出了不同時刻的解。表面圖只給出了負值范圍內（液滴域）的水平集 函數(shù)，同時等高線給出了液相線（Tm+Q、平均融化溫度（Tm）和固相線（Tm 9溫度 的等高線，箭頭圖給出了速度矢量。 i« 4 V*r1-.M .ill圖2在0.5 , 1, 1.5和2s時的模型計算結果從圖2中可以看出固化模型是如何作用的。氣相速度場和液相區(qū)域與預想結 果相同。在液、固相共存的多孔區(qū)域，速度分布仍然存在，但是更趨近于固相區(qū), 固相區(qū)計算出的速度幾乎為零。液滴中的熱

18、傳遞要比氣相中的快很多，因為液相的熱傳導系數(shù)更高。3.2固化速率分析在整個模擬時間段內，固化過程已經(jīng)發(fā)生，但是并不是整個液滴同時固化。 將fem結構導出，在 COMSOL Script或MATLAB中運行“ analysis.m”程序?？梢允褂谩?postint”函數(shù)對 '的Heavyside函數(shù)在空間域積分，算出不同時刻液 滴總面積和完全固化區(qū)域的面積t=1:101;A=posti nt(fem,' Hphi ' , ' Sol num ' ,t);As=posti nt(fem,' Hphfi*(2hs(T- Tm+e,0.001)'

19、 , ' Sol num ' ,t);Fs=As./A*100;plot(0.02*(t-1),A,0.02*(t-1),As);figureplot(0.02*(t-1),Fs);在以上代碼中，A是液滴總面積，它從一開始就應該保持常數(shù)0.196m2（質量守恒）；As是固相面積，通過f的Heaviside函數(shù)和另一個Heaviside函數(shù)（flc2hs） 在固相線（Tm 9附近相叉得到；Fs是液滴固化的百分比。從圖3中可以看出，隨著時間的增大，液滴總面積保持不變，所以質量始終 守恒，由于假設液相和固相密度相等，所以可以預計固化過程中體積都不會發(fā)生 變化。同時可以看出，模擬2s以

20、后有16%的液滴完全固化。，固體區(qū)域面積和固體比例圖圖3以時間函數(shù)形式表示的液滴面積下面來看相變潛熱對固化速率的影響，我們可以通過設定常數(shù)L= 0來忽略該影響。打開“ Options”菜單，選擇“ Constants”，設定L為0，點擊完成。 再次求解該問題。將fem結構體導出到COMSOL Script或MATLAB中，再次運行“analysis.m” 程序，可以得到圖4。圖4沒有相變潛熱時，以時間函數(shù)形式表示的液滴面積，固體區(qū)域面積和 固體比例圖經(jīng)過相同時間后，液滴固化比率達到了 45%。這與預想結果非常吻合，因 為沒有考慮相變過程中潛熱的釋放時間，所以比前一種情況中溫度降低的更快。 該例

21、題說明在固化問題中考率相變影響非常重要。4 .總結針對冶金中的一個重要例子，使用水平集方法和固化模型，求解了存在相變 的自由表面流動問題。與水平集法求解液體中氣泡融合的例子不同， 這里連續(xù)相 設定為氣相。即使該模型非常復雜（完全耦合了三個物理現(xiàn)象） ， COMSOL Multiphysics 的靈活性和易用性允許通過各種基本模塊來建立模型。例子最后對 更復雜和現(xiàn)實問題做了基本分析。致謝感謝Toronto大學Torstein Utigard教授的建議和指導。參考文獻1 V. R. Voller and C. Prakash, A fixed grid numerical modeling methodology for convection-diffusion mushy region phase-change problems, Int. J. Heat Mass Transfer 30 (1987) 1709.2 K. Morgan, A numerical analysis of freezing and melting with convect