第九章 存在相变的自由表面流动问题模拟——三相流T.L, MARINDeparment of Mining Engineering , University of Chile,Av. Tupper 2069 , Santiago, ChileE-mail:tmarin@ing.uchile.cl本章介绍了流体自由表面凝固过程的固定网格数值模型，并采用水平集方法来描述自由表面的流动Voller和Prakash对此方法进行了改进，考虑了液相凝固问题，包括对流和传递过程在该方法中，液体物性与温度相关，通过设定0到1之间的孔隙率使得体积力依赖于孔隙率和温度的变化，并通过对Navier-Stokes方程的修正来模拟液相或者固相此外，在传热方程中使用改进的热容表达式来计算融化潜热1．绪论液体在凝固时主要通过对流来传热，由于涉及到移动界面问题，所以通常难于模拟此外，对流体有效的控制方程对固体不再适用而且，如果问题中包含自由表面，那么流动将变得更加复杂，需要考虑相变或者凝固问题，这时都需要跟踪相界面的变化目前有很多种处理液体凝固的计算方法，考虑封闭空间中存在导热和自然对流的情况，Voller和Prakash对此作了总结[1]。
进展之一为采用变形网格来处理液体－固体界面的移动问题也可用固定网格以及焓变随温度的变化来处理该问题，这种方法的特点是较为简便，但在研究零速度封闭容器中的固体时会存在问题对于给定的计算单元，可以简单的设置速度为零，或者通过设置粘度为潜热容[2]函数的方法来实现[3]在这类情况下，当潜热容接近于零时，粘性会很大，这样才能模拟固体物性Voller和Prakash[1]研究了在一定温度范围内发生相变的情况，他们将流体描述为多孔介质，通过在Navier-Stokes方程中引入源项来模拟相变过程，用于研究速度为零的固体情况存在相变的自由表面流动问题的研究难点在于其实验和计算验证均很困难Pasandideh-Fard等[4]采用流体体积元方法(VOF)对锡液滴和水滴热表面上[5]的凝固进行了实验和计算研究，对自由表面及过程中的焓变进行跟踪他们的模型引入了液滴接触角以及液滴底部界面传热系数，该系数由实验测试得到并作为模拟的边界条件2．控制方程纯净流体凝固过程由熔点决定在液体冷却过程中，一旦达到这个温度，在温度继续降低之前液体开始释放相变潜热但是在多组分体系中，相变存在一个温度范围，从固体开始出现时的液体温度开始，一直到最后一种液体凝固时的温度为止。
在这种情况下，融化潜热在温度改变的同时不断释放相变期间固体组分(Fs)可以表示为温度的函数： (1)其中，T为系统温度，Tm为液相线和固相线的温度平均值，ε为液相线和固相线温度差值的一半因此，固相和液相温度为： (2) (3)前面已经提到，总的体系热容H由两部分组成，显热h和潜热△H显热通过以下方程计算： (4)潜热可以表示为温度的函数，根据之前对固体比例的定义： (5)导热和对流传热方程用系统温度形式表示为： (6)在这种情况下，释放潜热带来的影响可以通过对有效热容方程的重新定义来包含到热容项中： (7) (8)这里，表示平滑delta方程，在COMSOL Multiphysics中由fldc2hs函数建立注意，这里需要对该函数在整个温度范围内积分，但是只有到的温度范围只有在模拟纯净物质凝固过程的时候，才为零，此时将变为真正的Dirac’s Delta通过这种方法，系统总热容可以通过对方程(7)积分得到： (9)速度场和压力场的控制方程即Navier-Stokes方程： (10) (11)体积力F包括重力和水平集方法处理的表面张力，同时F也包含依赖于凝固过程的固体比例源项。
体积力分量如下： (12) (13)其中，是液体表面张力，是交界面曲率，是水平集函数，g是重力加速度，Sx和Sy项表述如下为了将固化过程并入体积力中，可以将液相看作是一种多孔介质，孔隙率依赖于温度全液相对应于孔隙率为1的状态，全固相对应于孔隙率为0的状态孔隙率定义如下： (14)源项定义如下： (15) (16)这里A是孔隙率的函数，根据以下方程定义： (17)下面介绍一下这些源项的作用当温度位于液相线以上时，系统处于全液态，源项取值为零，对Navier-Stokes方程没有任何改动多孔区域意味着温度处于液相线和固相线之间，A值增大并影响瞬态、对流和扩散项，流动方程近似相当于多孔介质中的Darcy定律当温度进一步降低，孔隙率接近于0时（固相），该源项决定了所有其它源项，使得速度值趋于零（固相）方程(17)中的常数C和q是任意选取的，依赖于求解的具体问题，通常C取较大值而q取较小值，避免当变为零时方程被0除水平集函数定义了两个初态相（例如，气和液）的交界面但是当固化（或者反问题中融化）时会出现第三相，一般出现在两个初始相的其中之一，所以需要修改源项Sx和Sy，使其严格作用在某一相中。
求解该问题非常简单，只需要将这几个源项乘以的Heaviside函数即可，出于我们对相的选择，我们限制了这些源项的作用范围这样，方程(10)中的体积力最终可以表示为： (18) (19)类似水平集方法，Heaviside函数同样用来定义连续相和非连续相间剧烈的物性变化这些物性包括密度，粘度和热传导率3．结果与讨论该模拟过程需要三个物理模型：不可压Navier-Stokes流动，对流和扩散传质过程，对流和热传导的传热过程求解该问题时用到了COMSOL Multiphysics的多物理场功能来同时处理和求解三个应用模式3.1 下降液滴的固化本节介绍一个液滴在冷表面的流动和固化过程为简单起见，将其简化为二维模型整个系统用左下角在x=y=0处的无量纲矩形域来表示，底面边界温度固定，其余边界绝热液体和连续相（环境）的初始边界条件用以下初始水平集函数来表示： (20)这样就在模型数值域的左下角形成四分之一个圆注意这种情况下液相值为负，连续相值为正为简化起见，在该问题中我们不考虑液/气交界面处表面张力的影响，集中处理固化问题但是就像在水平集方法中提到的，加入这些影响因素也非常简单。
打开COMSOL Multiphysics模型导航栏按照表1中的步骤建立模型注意Navier-Stokes模型中边界2（底边）为滑移/对称边界条件虽然正常边界条件是无滑移的，但是该假设对处理边界上水平集“phi”函数的对流非常有帮助，更容易模拟液滴和底面边界的接触过程同时，由于液滴从底面向上固化，速度自然会计算为零，由于对固相进行求解，所以在这种情况下滑移/对称边界条件并不算太脱离现实如果使用无滑移边界条件，就需要用到水平集函数的再次初始化表1 水平集三相固化模型模型导航栏选择2D维数COMSOL Multiphysics| Fluid Flow| Incompressible Navier- Stokes点击Multiphysics选项卡，添加选择Convection and Diffusion设定因变量为phi，添加Heat Transfer| Convection and Conduction添加，完成Draw菜单Specify Objects| Square输入width＝1，基准角x＝y＝0完成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)*Hphilambd=(T-Tm+e)/(2*e)*((T<=(Tm+e))*(T>=(Tm-e)))+(T>(Tm+e))Fs=1-lambdA=-C*(1-lambd)^2/(lambd^3+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为绝热边界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=vInit选项卡：T(t0)=Ti完成Mesh菜单：Mapped Mesh选择边界1，选中“Constrained edge element distribution”选项，输入40各边界基元选择边界2，选中“Constrained edge element distribution”选项，输入40各边界基元重绘网格，完成Solve菜单选择时间依赖求解器。
设定输出时间0:0.02:2，点击工具栏上的求解按钮（＝）根据方程(20)，初始水平集函数“phi0”用空间坐标“x”和“y”来定义，液相为正（或者在液相区域）、气相为负，如图1所示Heaviside函数“Hphi”用phi的负值来定义，所以在液滴内部为1，其它任何区域均为0液相和气相的物性用取决于水平集函数“phi”的Heaviside函数来定义根据方程(14)，用逻辑表达式定义孔隙率函数，该变量依赖于温度变量“T”和固相线、液相线的温度然后用该函数定义其它依赖于温度的物性图1 水平集函数初始条件的表面图该问题的新特性在于使用了映射网格它提供了棋盘形分布的四方形网格图2给出了不同时刻的解表面图只给出了负值范围内（液滴域）的水平集函数，同时等高线给出了液相线(Tm+ε)、平均融化温度(Tm)和固相线(Tm－ε)温度的等高线，箭头图给出了速度矢量图2 在0.5，1，1.5和2s时的模型计算结果从图2中可以看出固化模型是如何作用的气相速度场和液相区域与预想结果相同在液、固相共存的多孔区域，速度分布仍然存在，但是更趋近于固相区，固相区计算出的速度几乎为零液滴中的热传递要比气相中的快很多，因为液相的热传导系数更高。
3.2 固化速率分析在整个模拟时间段内，固化过程已经发生，但是并不是整个液滴同时固化将fem结构导出，在COMSOL Script或MATLAB中运行“analysis.m”程。