存在相变的自由表面流动问题模拟——三相流

1、第九章 存在相变的自由表面流动问题模拟三相流T.L, MARINDeparment of Mining Engineering , University of Chile,Av. Tupper 2069 , Santiago, ChileE-mail:tmarining.uchile.cl本章介绍了流体自由表面凝固过程的固定网格数值模型，并采用水平集方法来描述自由表面的流动。Voller和Prakash对此方法进行了改进，考虑了液相凝固问题，包括对流和传递过程。在该方法中，液体物性与温度相关，通过设定0到1之间的孔隙率使得体积力依赖于孔隙率和温度的变化，并通过对Navier-Stokes方程的

2、修正来模拟液相或者固相。此外，在传热方程中使用改进的热容表达式来计算融化潜热。1绪论液体在凝固时主要通过对流来传热，由于涉及到移动界面问题，所以通常难于模拟。此外，对流体有效的控制方程对固体不再适用。而且，如果问题中包含自由表面，那么流动将变得更加复杂，需要考虑相变或者凝固问题，这时都需要跟踪相界面的变化。目前有很多种处理液体凝固的计算方法，考虑封闭空间中存在导热和自然对流的情况，Voller和Prakash对此作了总结1。进展之一为采用变形网格来处理液体固体界面的移动问题。也可用固定网格以及焓变随温度的变化来处理该问题，这种方法的特点是较为简便，但在研究零速度封闭容器中的固体时会存在问题。对

3、于给定的计算单元，可以简单的设置速度为零，或者通过设置粘度为潜热容2函数的方法来实现3。在这类情况下，当潜热容接近于零时，粘性会很大，这样才能模拟固体物性。Voller和Prakash1研究了在一定温度范围内发生相变的情况，他们将流体描述为多孔介质，通过在Navier-Stokes方程中引入源项来模拟相变过程，用于研究速度为零的固体情况。存在相变的自由表面流动问题的研究难点在于其实验和计算验证均很困难。Pasandideh-Fard等4采用流体体积元方法(VOF)对锡液滴和水滴热表面上5的凝固进行了实验和计算研究，对自由表面及过程中的焓变进行跟踪。他们的模型引入了液滴接触角以及液滴底部界面传热

4、系数，该系数由实验测试得到并作为模拟的边界条件。2控制方程纯净流体凝固过程由熔点决定。在液体冷却过程中，一旦达到这个温度，在温度继续降低之前液体开始释放相变潜热。但是在多组分体系中，相变存在一个温度范围，从固体开始出现时的液体温度开始，一直到最后一种液体凝固时的温度为止。在这种情况下，融化潜热在温度改变的同时不断释放。相变期间固体组分(Fs)可以表示为温度的函数： (1)其中，T为系统温度，Tm为液相线和固相线的温度平均值，为液相线和固相线温度差值的一半。因此，固相和液相温度为： (2) (3)前面已经提到，总的体系热容H由两部分组成，显热h和潜热H。显热通过以下方程计算： (4)潜热可以表示

5、为温度的函数，根据之前对固体比例的定义： (5)导热和对流传热方程用系统温度形式表示为： (6)在这种情况下，释放潜热带来的影响可以通过对有效热容方程的重新定义来包含到热容项中： (7) (8)这里，表示平滑delta方程，在COMSOL Multiphysics中由fldc2hs函数建立。注意，这里需要对该函数在整个温度范围内积分，但是只有到的温度范围。只有在模拟纯净物质凝固过程的时候，才为零，此时将变为真正的Diracs Delta。通过这种方法，系统总热容可以通过对方程(7)积分得到： (9)速度场和压力场的控制方程即Navier-Stokes方程： (10) (11)体积力F包括重力和

6、水平集方法处理的表面张力，同时F也包含依赖于凝固过程的固体比例源项。体积力分量如下： (12) (13)其中，是液体表面张力，是交界面曲率，是水平集函数，g是重力加速度，Sx和Sy项表述如下。为了将固化过程并入体积力中，可以将液相看作是一种多孔介质，孔隙率依赖于温度。全液相对应于孔隙率为1的状态，全固相对应于孔隙率为0的状态。孔隙率定义如下： (14)源项定义如下： (15) (16)这里A是孔隙率的函数，根据以下方程定义： (17)下面介绍一下这些源项的作用。当温度位于液相线以上时，系统处于全液态，源项取值为零，对Navier-Stokes方程没有任何改动。多孔区域意味着温度处于液相线和固相

7、线之间，A值增大并影响瞬态、对流和扩散项，流动方程近似相当于多孔介质中的Darcy定律。当温度进一步降低，孔隙率接近于0时（固相），该源项决定了所有其它源项，使得速度值趋于零（固相）。方程(17)中的常数C和q是任意选取的，依赖于求解的具体问题，通常C取较大值而q取较小值，避免当变为零时方程被0除。水平集函数定义了两个初态相（例如，气和液）的交界面。但是当固化（或者反问题中融化）时会出现第三相，一般出现在两个初始相的其中之一，所以需要修改源项Sx和Sy，使其严格作用在某一相中。求解该问题非常简单，只需要将这几个源项乘以的Heaviside函数即可，出于我们对相的选择，我们限制了这些源项的作用范

8、围。这样，方程(10)中的体积力最终可以表示为： (18) (19)类似水平集方法，Heaviside函数同样用来定义连续相和非连续相间剧烈的物性变化。这些物性包括密度，粘度和热传导率。3结果与讨论该模拟过程需要三个物理模型：不可压Navier-Stokes流动，对流和扩散传质过程，对流和热传导的传热过程。求解该问题时用到了COMSOL Multiphysics的多物理场功能来同时处理和求解三个应用模式。3.1 下降液滴的固化本节介绍一个液滴在冷表面的流动和固化过程。为简单起见，将其简化为二维模型。整个系统用左下角在x=y=0处的无量纲矩形域来表示，底面边界温度固定，其余边界绝热。液体和连续相

9、（环境）的初始边界条件用以下初始水平集函数来表示： (20)这样就在模型数值域的左下角形成四分之一个圆。注意这种情况下液相值为负，连续相值为正。为简化起见，在该问题中我们不考虑液/气交界面处表面张力的影响，集中处理固化问题。但是就像在水平集方法中提到的，加入这些影响因素也非常简单。打开COMSOL Multiphysics模型导航栏。按照表1中的步骤建立模型。注意Navier-Stokes模型中边界2（底边）为滑移/对称边界条件。虽然正常边界条件是无滑移的，但是该假设对处理边界上水平集“phi”函数的对流非常有帮助，更容易模拟液滴和底面边界的接触过程。同时，由于液滴从底面向上固化，速度自然会计

10、算为零，由于对固相进行求解，所以在这种情况下滑移/对称边界条件并不算太脱离现实。如果使用无滑移边界条件，就需要用到水平集函数的再次初始化。表1 水平集三相固化模型模型导航栏选择2D维数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，点击工具栏上的求解按钮（）。根据方程(20)，初始水平集函数“phi0”用空间坐标“x”和“y”来定义，液相为正（或者在液相区域）、气相为负，如图1所示。Heaviside函数“Hphi”用phi的负值来定义，所以在液滴内部为1，其它任何区域均为0。液相和气相的物性用取决于水平集函数“phi”的Heaviside函数来定义。根据方程(14)，用逻辑表达式定义孔隙率函数，该变量依赖于温度变量“T”和固相线、液相线的温度。然后用该函数定义其它依赖于温度的物

15、性。图1 水平集函数初始条件的表面图该问题的新特性在于使用了映射网格。它提供了棋盘形分布的四方形网格。图2给出了不同时刻的解。表面图只给出了负值范围内（液滴域）的水平集函数，同时等高线给出了液相线(Tm+)、平均融化温度(Tm)和固相线(Tm)温度的等高线，箭头图给出了速度矢量。图2 在0.5，1，1.5和2s时的模型计算结果从图2中可以看出固化模型是如何作用的。气相速度场和液相区域与预想结果相同。在液、固相共存的多孔区域，速度分布仍然存在，但是更趋近于固相区，固相区计算出的速度几乎为零。液滴中的热传递要比气相中的快很多，因为液相的热传导系数更高。3.2 固化速率分析在整个模拟时间段内，固化过

16、程已经发生，但是并不是整个液滴同时固化。将fem结构导出，在COMSOL Script或MATLAB中运行“analysis.m”程序。可以使用“postint”函数对的Heavyside函数在空间域积分，算出不同时刻液滴总面积和完全固化区域的面积。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、面积，它从一开始就应该保持常数0.196m2（质量守恒）；As是固相面积，通过f的Heaviside函数和另一个Heaviside函数（flc2hs）在固相线(Tm)附近相叉得到；Fs是液滴固化的百分比。从图3中可以看出，随着时间的增大，液滴总面积保持不变，所以质量始终守恒，由于假设液相和固相密度相等，所以可以预计固化过程中体积都不会发生变化。同时可以看出，模拟2s以后有16的液滴完全固化。图3 以时间函数形式表示的液滴面积，固体区域面积和固体比例图下面来看相变潜热对固化速率的影响，我们可以通过设定常数L0来忽略该影响。打开“Options”菜单，选择“Constants”，设定L为0，点击完

18、成。再次求解该问题。将fem结构体导出到COMSOL Script或MATLAB中，再次运行“analysis.m”程序，可以得到图4。图4 没有相变潜热时，以时间函数形式表示的液滴面积，固体区域面积和固体比例图经过相同时间后，液滴固化比率达到了45。这与预想结果非常吻合，因为没有考虑相变过程中潜热的释放时间，所以比前一种情况中温度降低的更快。该例题说明在固化问题中考率相变影响非常重要。4总结针对冶金中的一个重要例子，使用水平集方法和固化模型，求解了存在相变的自由表面流动问题。与水平集法求解液体中气泡融合的例子不同，这里连续相设定为气相。即使该模型非常复杂（完全耦合了三个物理现象），COMSO

19、L Multiphysics的灵活性和易用性允许通过各种基本模块来建立模型。例子最后对更复杂和现实问题做了基本分析。致谢感谢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