存在相變的自由表面流動問題模擬——三相流

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

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

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

4、系數(shù)，該系數(shù)由實驗測試得到并作為模擬的邊界條件。2控制方程純凈流體凝固過程由熔點決定。在液體冷卻過程中，一旦達到這個溫度，在溫度繼續(xù)降低之前液體開始釋放相變潛熱。但是在多組分體系中，相變存在一個溫度范圍，從固體開始出現(xiàn)時的液體溫度開始，一直到最后一種液體凝固時的溫度為止。在這種情況下，融化潛熱在溫度改變的同時不斷釋放。相變期間固體組分(Fs)可以表示為溫度的函數(shù)： (1)其中，T為系統(tǒng)溫度，Tm為液相線和固相線的溫度平均值，為液相線和固相線溫度差值的一半。因此，固相和液相溫度為： (2) (3)前面已經提到，總的體系熱容H由兩部分組成，顯熱h和潛熱H。顯熱通過以下方程計算： (4)潛熱可以表示

5、為溫度的函數(shù)，根據(jù)之前對固體比例的定義： (5)導熱和對流傳熱方程用系統(tǒng)溫度形式表示為： (6)在這種情況下，釋放潛熱帶來的影響可以通過對有效熱容方程的重新定義來包含到熱容項中： (7) (8)這里，表示平滑delta方程，在COMSOL Multiphysics中由fldc2hs函數(shù)建立。注意，這里需要對該函數(shù)在整個溫度范圍內積分，但是只有到的溫度范圍。只有在模擬純凈物質凝固過程的時候，才為零，此時將變?yōu)檎嬲腄iracs Delta。通過這種方法，系統(tǒng)總熱容可以通過對方程(7)積分得到： (9)速度場和壓力場的控制方程即Navier-Stokes方程： (10) (11)體積力F包括重力和

6、水平集方法處理的表面張力，同時F也包含依賴于凝固過程的固體比例源項。體積力分量如下： (12) (13)其中，是液體表面張力，是交界面曲率，是水平集函數(shù)，g是重力加速度，Sx和Sy項表述如下。為了將固化過程并入體積力中，可以將液相看作是一種多孔介質，孔隙率依賴于溫度。全液相對應于孔隙率為1的狀態(tài)，全固相對應于孔隙率為0的狀態(tài)。孔隙率定義如下： (14)源項定義如下： (15) (16)這里A是孔隙率的函數(shù)，根據(jù)以下方程定義： (17)下面介紹一下這些源項的作用。當溫度位于液相線以上時，系統(tǒng)處于全液態(tài)，源項取值為零，對Navier-Stokes方程沒有任何改動。多孔區(qū)域意味著溫度處于液相線和固相

7、線之間，A值增大并影響瞬態(tài)、對流和擴散項，流動方程近似相當于多孔介質中的Darcy定律。當溫度進一步降低，孔隙率接近于0時（固相），該源項決定了所有其它源項，使得速度值趨于零（固相）。方程(17)中的常數(shù)C和q是任意選取的，依賴于求解的具體問題，通常C取較大值而q取較小值，避免當變?yōu)榱銜r方程被0除。水平集函數(shù)定義了兩個初態(tài)相（例如，氣和液）的交界面。但是當固化（或者反問題中融化）時會出現(xiàn)第三相，一般出現(xiàn)在兩個初始相的其中之一，所以需要修改源項Sx和Sy，使其嚴格作用在某一相中。求解該問題非常簡單，只需要將這幾個源項乘以的Heaviside函數(shù)即可，出于我們對相的選擇，我們限制了這些源項的作用范

8、圍。這樣，方程(10)中的體積力最終可以表示為： (18) (19)類似水平集方法，Heaviside函數(shù)同樣用來定義連續(xù)相和非連續(xù)相間劇烈的物性變化。這些物性包括密度，粘度和熱傳導率。3結果與討論該模擬過程需要三個物理模型：不可壓Navier-Stokes流動，對流和擴散傳質過程，對流和熱傳導的傳熱過程。求解該問題時用到了COMSOL Multiphysics的多物理場功能來同時處理和求解三個應用模式。3.1 下降液滴的固化本節(jié)介紹一個液滴在冷表面的流動和固化過程。為簡單起見，將其簡化為二維模型。整個系統(tǒng)用左下角在x=y=0處的無量綱矩形域來表示，底面邊界溫度固定，其余邊界絕熱。液體和連續(xù)相

9、（環(huán)境）的初始邊界條件用以下初始水平集函數(shù)來表示： (20)這樣就在模型數(shù)值域的左下角形成四分之一個圓。注意這種情況下液相值為負，連續(xù)相值為正。為簡化起見，在該問題中我們不考慮液/氣交界面處表面張力的影響，集中處理固化問題。但是就像在水平集方法中提到的，加入這些影響因素也非常簡單。打開COMSOL Multiphysics模型導航欄。按照表1中的步驟建立模型。注意Navier-Stokes模型中邊界2（底邊）為滑移/對稱邊界條件。雖然正常邊界條件是無滑移的，但是該假設對處理邊界上水平集“phi”函數(shù)的對流非常有幫助，更容易模擬液滴和底面邊界的接觸過程。同時，由于液滴從底面向上固化，速度自然會計

10、算為零，由于對固相進行求解，所以在這種情況下滑移/對稱邊界條件并不算太脫離現(xiàn)實。如果使用無滑移邊界條件，就需要用到水平集函數(shù)的再次初始化。表1 水平集三相固化模型模型導航欄選擇2D維數(shù)COMSOL Multiphysics| Fluid Flow| Incompressible Navier- Stokes。點擊Multiphysics選項卡，添加。選擇Convection and Diffusion。設定因變量為phi，添加。Heat Transfer| Convection and Conduction。添加，完成。Draw菜單Specify Objects| Square輸入width1

11、，基準角xy0。完成Options菜單：Constantsx0=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=0.5, kg=0.01, L=1, C=1600, q=0.001, Ttm=0Options菜單：Expressionsphi0=sqrt(x-x0)2+(y-y0)2)-rHphi=(1+tanh(-phi/n)/2rho=rhog+(rhol-rhog)*Hphieta=etag+(etal-etag)*Hphik=kg+(kl-kg)*Hphi

12、lambd=(T-Tm+e)/(2*e)*(T=(Tm-e)+(T(Tm+e)Fs=1-lambdA=-C*(1-lambd)2/(lambd3+q)Sx=-A*u, Sy=-A*vcpH=L*fldc2hs(T-Tm,e)cpT=cp+cpHPhysics菜單：Point settingsns模式。確認3點的點約束設定為0壓力情況。Physics菜單：Boundary settingsns模式：邊界1和2：slip/symmetry，邊界3和4：no-slipchcd模式：設定所有邊界條件為Insulation/symmetrycc模式：邊界2設定為溫度，在溫度欄中輸入TC設定邊界1，3和4

13、為絕熱邊界。Physics菜單：Subdomain Settingsns模式：=rho, =eta, Fx=-Sx*Hphi, Fy=-Sy*Hphi+rho*gychcd模式：D(各向同性)=0, R=0, u=u, v=vInit選項卡：phi(t0)=phi0cc模式：k(各向同性)=k, =0, Cp=cpT, Q=0, u=u, v=v。Init選項卡：T(t0)=Ti。完成Mesh菜單：Mapped Mesh選擇邊界1，選中“Constrained edge element distribution”選項，輸入40各邊界基元。選擇邊界2，選中“Constrained edge el

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

15、性。圖1 水平集函數(shù)初始條件的表面圖該問題的新特性在于使用了映射網格。它提供了棋盤形分布的四方形網格。圖2給出了不同時刻的解。表面圖只給出了負值范圍內（液滴域）的水平集函數(shù)，同時等高線給出了液相線(Tm+)、平均融化溫度(Tm)和固相線(Tm)溫度的等高線，箭頭圖給出了速度矢量。圖2 在0.5，1，1.5和2s時的模型計算結果從圖2中可以看出固化模型是如何作用的。氣相速度場和液相區(qū)域與預想結果相同。在液、固相共存的多孔區(qū)域，速度分布仍然存在，但是更趨近于固相區(qū)，固相區(qū)計算出的速度幾乎為零。液滴中的熱傳遞要比氣相中的快很多，因為液相的熱傳導系數(shù)更高。3.2 固化速率分析在整個模擬時間段內，固化過

16、程已經發(fā)生，但是并不是整個液滴同時固化。將fem結構導出，在COMSOL Script或MATLAB中運行“analysis.m”程序?？梢允褂谩皃ostint”函數(shù)對的Heavyside函數(shù)在空間域積分，算出不同時刻液滴總面積和完全固化區(qū)域的面積。t=1:101;A=postint(fem,Hphi,Solnum,t);As=postint(fem,Hphi*(1-flc2hs(T-Tm+e,0.001),Solnum,t);Fs=As./A*100;plot(0.02*(t-1),A,0.02*(t-1),As);figureplot(0.02*(t-1),Fs);在以上代碼中，A是液滴總

17、面積，它從一開始就應該保持常數(shù)0.196m2（質量守恒）；As是固相面積，通過f的Heaviside函數(shù)和另一個Heaviside函數(shù)（flc2hs）在固相線(Tm)附近相叉得到；Fs是液滴固化的百分比。從圖3中可以看出，隨著時間的增大，液滴總面積保持不變，所以質量始終守恒，由于假設液相和固相密度相等，所以可以預計固化過程中體積都不會發(fā)生變化。同時可以看出，模擬2s以后有16的液滴完全固化。圖3 以時間函數(shù)形式表示的液滴面積，固體區(qū)域面積和固體比例圖下面來看相變潛熱對固化速率的影響，我們可以通過設定常數(shù)L0來忽略該影響。打開“Options”菜單，選擇“Constants”，設定L為0，點擊完

18、成。再次求解該問題。將fem結構體導出到COMSOL Script或MATLAB中，再次運行“analysis.m”程序，可以得到圖4。圖4 沒有相變潛熱時，以時間函數(shù)形式表示的液滴面積，固體區(qū)域面積和固體比例圖經過相同時間后，液滴固化比率達到了45。這與預想結果非常吻合，因為沒有考慮相變過程中潛熱的釋放時間，所以比前一種情況中溫度降低的更快。該例題說明在固化問題中考率相變影響非常重要。4總結針對冶金中的一個重要例子，使用水平集方法和固化模型，求解了存在相變的自由表面流動問題。與水平集法求解液體中氣泡融合的例子不同，這里連續(xù)相設定為氣相。即使該模型非常復雜（完全耦合了三個物理現(xiàn)象），COMSO

19、L 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 convection, Comput. Meth. Appl. Eng. 28 (1981) 275.3 D. K. Gartling, Finite elemen