存在相變的自由表面流動(dòng)問(wèn)題模擬——三相流

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

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

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

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

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

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

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

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

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

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

11、，基準(zhǔn)角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模式。確認(rèn)3點(diǎn)的點(diǎn)約束設(shè)定為0壓力情況。Physics菜單：Boundary settingsns模式：邊界1和2：slip/symmetry，邊界3和4：no-slipchcd模式：設(shè)定所有邊界條件為Insulation/symmetrycc模式：邊界2設(shè)定為溫度，在溫度欄中輸入TC設(shè)定邊界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選項(xiàng)卡：phi(t0)=phi0cc模式：k(各向同性)=k, =0, Cp=cpT, Q=0, u=u, v=v。Init選項(xiàng)卡：T(t0)=Ti。完成Mesh菜單：Mapped Mesh選擇邊界1，選中“Constrained edge element distribution”選項(xiàng)，輸入40各邊界基元。選擇邊界2，選中“Constrained edge el

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

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

16、程已經(jīng)發(fā)生，但是并不是整個(gè)液滴同時(shí)固化。將fem結(jié)構(gòu)導(dǎo)出，在COMSOL Script或MATLAB中運(yùn)行“analysis.m”程序。可以使用“postint”函數(shù)對(duì)的Heavyside函數(shù)在空間域積分，算出不同時(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、面積，它從一開始就應(yīng)該保持常數(shù)0.196m2（質(zhì)量守恒）；As是固相面積，通過(guò)f的Heaviside函數(shù)和另一個(gè)Heaviside函數(shù)（flc2hs）在固相線(Tm)附近相叉得到；Fs是液滴固化的百分比。從圖3中可以看出，隨著時(shí)間的增大，液滴總面積保持不變，所以質(zhì)量始終守恒，由于假設(shè)液相和固相密度相等，所以可以預(yù)計(jì)固化過(guò)程中體積都不會(huì)發(fā)生變化。同時(shí)可以看出，模擬2s以后有16的液滴完全固化。圖3 以時(shí)間函數(shù)形式表示的液滴面積，固體區(qū)域面積和固體比例圖下面來(lái)看相變潛熱對(duì)固化速率的影響，我們可以通過(guò)設(shè)定常數(shù)L0來(lái)忽略該影響。打開“Options”菜單，選擇“Constants”，設(shè)定L為0，點(diǎn)擊完

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

19、L Multiphysics的靈活性和易用性允許通過(guò)各種基本模塊來(lái)建立模型。例子最后對(duì)更復(fù)雜和現(xiàn)實(shí)問(wèn)題做了基本分析。致謝感謝Toronto大學(xué)Torstein Utigard教授的建議和指導(dǎo)。參考文獻(xiàn)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